text:
"""
    AbstractTrees

This package is intended to provide an abstract interface for working
with tree structures.
Though the package itself is not particularly sophisticated, it defines
the interface that can be used by other packages to talk about trees.
"""
module AbstractTrees

using Base: HasLength, SizeUnknown, HasEltype, EltypeUnknown

include("traits.jl")
include("base.jl")
include("indexing.jl")
include("cursors.jl")
include("iteration.jl")
include("builtins.jl")
include("printing.jl")

#interface
export ParentLinks, StoredParents, ImplicitParents
export SiblingLinks, StoredSiblings, ImplicitSiblings
export ChildIndexing, IndexedChildren, NonIndexedChildren
export NodeType, HasNodeType, NodeTypeUnknown
export nodetype, nodevalue, nodevalues, children, parentlinks, siblinglinks, childindexing, childtype, childrentype
#extended interface
export nextsibling, prevsibling

# properties
export ischild, isroot, isroot, intree, isdescendant, treesize, treebreadth, treeheight, descendleft, getroot

# cursors
export TreeCursor, TrivialCursor, ImplicitCursor, SiblingCursor

#indexing
export getdescendant, childindices, rootindex, parentindex, nextsiblingindex, prevsiblingindex, IndexNode

# iteration
export TreeIterator, PreOrderDFS, PostOrderDFS, Siblings, Leaves, StatelessBFS, MapNode
export treemap

# printing
export print_tree, repr_tree

end # module


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

#! format: off

"""
    mutable struct RenewableDispatch <: RenewableGen
        name::String
        available::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        prime_mover::PrimeMovers
        reactive_power_limits::Union{Nothing, Min_Max}
        power_factor::Float64
        operation_cost::TwoPartCost
        base_power::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `active_power::Float64`
- `reactive_power::Float64`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `reactive_power_limits::Union{Nothing, Min_Max}`
- `power_factor::Float64`, validation range: `(0, 1)`, action if invalid: `error`
- `operation_cost::TwoPartCost`: Operation Cost of Generation [`TwoPartCost`](@ref)
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RenewableDispatch <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    reactive_power_limits::Union{Nothing, Min_Max}
    power_factor::Float64
    "Operation Cost of Generation [`TwoPartCost`](@ref)"
    operation_cost::TwoPartCost
    "Base power of the unit in MVA"
    base_power::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RenewableDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, reactive_power_limits, power_factor, operation_cost, base_power, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    RenewableDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, reactive_power_limits, power_factor, operation_cost, base_power, services, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function RenewableDispatch(; name, available, bus, active_power, reactive_power, rating, prime_mover, reactive_power_limits, power_factor, operation_cost, base_power, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    RenewableDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, reactive_power_limits, power_factor, operation_cost, base_power, services, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function RenewableDispatch(::Nothing)
    RenewableDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        rating=0.0,
        prime_mover=PrimeMovers.OT,
        reactive_power_limits=nothing,
        power_factor=1.0,
        operation_cost=TwoPartCost(nothing),
        base_power=0.0,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`RenewableDispatch`](@ref) `name`."""
get_name(value::RenewableDispatch) = value.name
"""Get [`RenewableDispatch`](@ref) `available`."""
get_available(value::RenewableDispatch) = value.available
"""Get [`RenewableDispatch`](@ref) `bus`."""
get_bus(value::RenewableDispatch) = value.bus
"""Get [`RenewableDispatch`](@ref) `active_power`."""
get_active_power(value::RenewableDispatch) = get_value(value, value.active_power)
"""Get [`RenewableDispatch`](@ref) `reactive_power`."""
get_reactive_power(value::RenewableDispatch) = get_value(value, value.reactive_power)
"""Get [`RenewableDispatch`](@ref) `rating`."""
get_rating(value::RenewableDispatch) = get_value(value, value.rating)
"""Get [`RenewableDispatch`](@ref) `prime_mover`."""
get_prime_mover(value::RenewableDispatch) = value.prime_mover
"""Get [`RenewableDispatch`](@ref) `reactive_power_limits`."""
get_reactive_power_limits(value::RenewableDispatch) = get_value(value, value.reactive_power_limits)
"""Get [`RenewableDispatch`](@ref) `power_factor`."""
get_power_factor(value::RenewableDispatch) = value.power_factor
"""Get [`RenewableDispatch`](@ref) `operation_cost`."""
get_operation_cost(value::RenewableDispatch) = value.operation_cost
"""Get [`RenewableDispatch`](@ref) `base_power`."""
get_base_power(value::RenewableDispatch) = value.base_power
"""Get [`RenewableDispatch`](@ref) `services`."""
get_services(value::RenewableDispatch) = value.services
"""Get [`RenewableDispatch`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::RenewableDispatch) = value.dynamic_injector
"""Get [`RenewableDispatch`](@ref) `ext`."""
get_ext(value::RenewableDispatch) = value.ext
"""Get [`RenewableDispatch`](@ref) `time_series_container`."""
get_time_series_container(value::RenewableDispatch) = value.time_series_container
"""Get [`RenewableDispatch`](@ref) `internal`."""
get_internal(value::RenewableDispatch) = value.internal

"""Set [`RenewableDispatch`](@ref) `available`."""
set_available!(value::RenewableDispatch, val) = value.available = val
"""Set [`RenewableDispatch`](@ref) `bus`."""
set_bus!(value::RenewableDispatch, val) = value.bus = val
"""Set [`RenewableDispatch`](@ref) `active_power`."""
set_active_power!(value::RenewableDispatch, val) = value.active_power = set_value(value, val)
"""Set [`RenewableDispatch`](@ref) `reactive_power`."""
set_reactive_power!(value::RenewableDispatch, val) = value.reactive_power = set_value(value, val)
"""Set [`RenewableDispatch`](@ref) `rating`."""
set_rating!(value::RenewableDispatch, val) = value.rating = set_value(value, val)
"""Set [`RenewableDispatch`](@ref) `prime_mover`."""
set_prime_mover!(value::RenewableDispatch, val) = value.prime_mover = val
"""Set [`RenewableDispatch`](@ref) `reactive_power_limits`."""
set_reactive_power_limits!(value::RenewableDispatch, val) = value.reactive_power_limits = set_value(value, val)
"""Set [`RenewableDispatch`](@ref) `power_factor`."""
set_power_factor!(value::RenewableDispatch, val) = value.power_factor = val
"""Set [`RenewableDispatch`](@ref) `operation_cost`."""
set_operation_cost!(value::RenewableDispatch, val) = value.operation_cost = val
"""Set [`RenewableDispatch`](@ref) `base_power`."""
set_base_power!(value::RenewableDispatch, val) = value.base_power = val
"""Set [`RenewableDispatch`](@ref) `services`."""
set_services!(value::RenewableDispatch, val) = value.services = val
"""Set [`RenewableDispatch`](@ref) `ext`."""
set_ext!(value::RenewableDispatch, val) = value.ext = val
"""Set [`RenewableDispatch`](@ref) `time_series_container`."""
set_time_series_container!(value::RenewableDispatch, val) = value.time_series_container = val

ext: jl


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra
using Quaternions

# ## Simulation
s = get_simulation("centroidal_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
	# joinpath(module_dir(), "src/dynamics/centroidal_quadruped/gaits/inplace_trot_v4.jld2"),
	joinpath(@__DIR__, "reference/inplace_bounce.jld2"),
    # joinpath(module_dir(), "src/dynamics/centroidal_quadruped/gaits/stand_euler_v0.jld2"),
    load_type = :split_traj_alt));

H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 3000
κ_mpc = 2.0e-4

v0 = -0.3
obj = TrackingVelocityObjective(model, env, H_mpc,
    v = [Diagonal(1e-3 * [[1,1,1]; 1e+3*[1,1,1]; fill([1,1,1], 4)...]) for t = 1:H_mpc],
	q = [relative_state_cost(1e-0*[1e-5,1e-5,1], 3e-1*[1,1,1], 1e-0*[0.2,0.2,1]) for t = 1:H_mpc],
	u = [Diagonal(3e-3 * vcat(fill([1,1,1], 4)...)) for t = 1:H_mpc],
	v_target = [1/ref_traj.h * [v0;0;0; 0;0;0; v0;0;0; v0;0;0; v0;0;0; v0;0;0] for t = 1:H_mpc],)

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					κ_tol = κ_mpc,
					r_tol = 1.0e-4, # TODO maybe relax this parameter
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1e5),
    n_opts = NewtonOptions(
        r_tol = 3e-5,
		solver=:ldl_solver,
        max_iter = 5),
    mpc_opts = CIMPCOptions(
		# live_plotting=true
		));

# ## Disturbances
idx_d1 = 20
idx_d2 = idx_d1 + 350
idx_d3 = idx_d2 + 350
idx_d4 = idx_d3 + 350
idx_d5 = idx_d4 + 350
idx_d6 = idx_d5 + 350
idx = [idx_d1, idx_d2, idx_d3, idx_d4, idx_d5, idx_d6]
impulses = [[0,0,10.0], [0,0,-10.0], [0,5,0.0], [0,-5,0.0], [5,0,0.0], [-5,0,0.0]]
# impulses = [[0,10,0.0], [0,0,0.0], [0,0,0.0], [0,0,0.0], [0,0,0.0], [0,0,0.0]]
d = impulse_disturbances(impulses, idx)

w = [[0.05,0.0,0.0] for i=1:H_sim/N_sample]
w = [[0,0.07,0.0] for i=1:H_sim/N_sample]
w = [[0,0.0,0.0] for i=1:H_sim/N_sample]
d = open_loop_disturbances(w, N_sample)

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator

# u_ref = deepcopy(ref_traj.u[1])
# p = open_loop_policy([u_ref for i=1:1000], N_sample=N_sample)
# p = open_loop_policy(deepcopy(ref_traj.u), N_sample=N_sample)
# sim = simulator(s, H_sim, h=h_sim, dist=d);
sim = simulator(s, H_sim, h=h_sim, policy=p, dist=d);

using BenchmarkTools
# ## Simulate
q1_sim0 = deepcopy(q1_sim)
# q1_sim0[4:6] += [0.02, 0.04, 0.02]
Δx = 0.0*[0.1, 0.2, 0.1]
q1_sim0[1:3] += Δx
q1_sim0[7:9] += Δx
q1_sim0[10:12] += Δx
q1_sim0[13:15] += Δx
q1_sim0[16:18] += Δx
RoboDojo.simulate!(sim, q1_sim0, v1_sim)

# # ## Visualizer
# vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.open(vis)

# ## Visualize
set_light!(vis)
set_floor!(vis, grid=true)
set_background!(vis)
anim = visualize!(vis, model, sim.traj.q; Δt=h_sim)

# # ## Timing result
# # Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio
plot(sim.stats.policy_time)

# convert_frames_to_video_and_gif("centroidal_inplace_bounce_forward")


text:
function MStep_mfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)
    #update = MFAModelM(model.pivec,model.mu,model.B,model.D,model.g,model.q,model.tau,model.logL,model.BIC)
    Mquantities.n_i = sum(model.tau, dims = 1)

    #CM-Step 1: Updating pi and mu_i
    model.pivec = Mquantities.n_i ./ Squantities.n
    #mu_new = Array{Float64,2}(undef,p,g)
    @inbounds for i in 1:model.g
        @views model.mu[:,i] = sum(Squantities.Y .* model.tau[:,i], dims = 1) ./ sum(model.tau[:,i])
    end

    #CM-Step 2: Updating q
    #Vtilde = @MArray zeros(Float64,p,p,model.g)
    #Lambda = @MMatrix zeros(Float64,q_max,model.g)
    #Lambda_new = Array{Float64,2}(undef,q_max,model.g)
    #eigenvecs =  Array{Any,1}(undef,model.g)
    @inbounds for i in 1:model.g
        @views Mquantities.Ymu = Squantities.Y .- model.mu[:,i]'
        @views Mquantities.Vi = (Array(Mquantities.Ymu') * (Mquantities.Ymu .* model.tau[:,i]))/Mquantities.n_i[i]
        @views Mquantities.D_is = Diagonal(1 ./ sqrt.(model.D[:,:,i]))
        Mquantities.Vtilde[:,:,i] = Mquantities.D_is*Mquantities.Vi*Mquantities.D_is
        @views Mquantities.V_eig_temp = eigvals(Array(Mquantities.Vtilde[:,:,i]))
        #if any(imag.(V_eig_temp) .!= 0)
        #    print("uh oh")
        #    print(imag.(V_eig_temp))
        #    print("converting to real number >:(")
        #end

        if any(imag(Mquantities.V_eig_temp) .> 0)
            println("Complex eigenvalues appeared!!")
            model.status = "failed"
            return
        end

       # Mquantities.V_eig = real.(Mquantities.V_eig_temp)
        Mquantities.eig_order = sortperm(Mquantities.V_eig_temp, rev = true)
        @views Mquantities.Lambda[:,i] = Mquantities.V_eig_temp[Mquantities.eig_order][1:Squantities.qmax]
        @views Mquantities.eigenvecs[i] = real.(eigvecs(Array(Mquantities.Vtilde[:,:,i]))[:,Mquantities.eig_order])
        if any(Mquantities.Lambda[:,i] .< 0)
            println("negative eigenvalues appeared!!")
            Mquantities.Lambda[(Mquantities.Lambda[:,i] .< 0),i] .= 1e-20
        end
    end
    #CM-Step 2.5: Optionally updating q
    if Squantities.update_q #Fix this so that we don't accidentally change the value of q when we shouldnt.
        Mquantities.qup.Lambda_cs = cumsum(log.(Mquantities.Lambda) - Mquantities.Lambda + ones(size(Mquantities.Lambda)), dims = 1)
        Mquantities.qup.n_param = Int.((model.g - 1) .+ 2*model.g*Squantities.p .+ model.g * (Squantities.p .* (1:Squantities.qmax) .- (1:Squantities.qmax) .* ((1:Squantities.qmax) .- 1)/2))
        Mquantities.qup.q_scores = sum(Mquantities.qup.Lambda_cs .* Mquantities.n_i, dims = 2) .+ log(Squantities.n)*Mquantities.qup.n_param
        Mquantities.qup.q_up = argmin(vec(Mquantities.qup.q_scores))
        Mquantities.qup.q_change = (Mquantities.qup.q_up != model.q)
        model.q = Mquantities.qup.q_up
        if Mquantities.qup.q_change
            Mquantities.qup.B_update = zeros(ComplexF64,Squantities.p,model.q,model.g)
        end
    end

    #CM-Step 3: Updating B_i, then D_i in same loop.
    Mquantities.qi_break = false
    @inbounds for i in 1:model.g
        #Update B
        @views Mquantities.qi = sum(Mquantities.Lambda[:,i] .> 1)
        if Mquantities.qi > model.q
            Mquantities.qi = model.q
        end

        if Mquantities.qi == 0
            Mquantities.qi_break = true
           break
        end
        @views Mquantities.Uqi = Mquantities.eigenvecs[i][:,1:Mquantities.qi]
        @views Mquantities.Ri = I(model.q)[1:Mquantities.qi,:]
        if Squantities.update_q
            if Mquantities.qup.q_change
                @views Mquantities.qup.B_update[:,:,i] = sqrt.(model.D[:,:,i]) * Mquantities.Uqi *
                Diagonal(sqrt.(Mquantities.Lambda[1:Mquantities.qi,i] - vec(ones(1,Mquantities.qi)))) * Mquantities.Ri
            else
                @views model.B[:,:,i] = sqrt.(model.D[:,:,i]) * Mquantities.Uqi *
            Diagonal(sqrt.(Mquantities.Lambda[1:Mquantities.qi,i] - vec(ones(1,Mquantities.qi)))) * Mquantities.Ri
            end
        else
            @views model.B[:,:,i] = sqrt.(model.D[:,:,i]) * Mquantities.Uqi *
            Diagonal(sqrt.(Mquantities.Lambda[1:Mquantities.qi,i] - vec(ones(1,Mquantities.qi)))) * Mquantities.Ri
        end

        #Update D
        Mquantities.C = Array(1.0*I(Squantities.p)); Mquantities.Us = Mquantities.Uqi*(Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi,i])-I(Mquantities.qi)); Mquantities.V = Array(Mquantities.Uqi')
        Mquantities.bvec = zeros(Squantities.p,1)
        #Mquantities.psitilde = Array{Float64,2}(undef,p,1)
        @inbounds for r in 1:Squantities.p
            @views Mquantities.c = Mquantities.C[1:r,r]
            @views Mquantities.v = Mquantities.V[:,1:r]*Mquantities.c
            @views Mquantities.bvec[r:Squantities.p] = Mquantities.Us[r:Squantities.p,:]*Mquantities.v
            Mquantities.b = Mquantities.bvec[r] + 1
            @views Mquantities.a = Array(Mquantities.v')*Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi, i] .- Mquantities.Lambda[1:Mquantities.qi,i])*Mquantities.v + Array(Mquantities.c')*Mquantities.Vtilde[1:r,1:r,i] * Mquantities.c
            Mquantities.psitilde[r] = max((((Mquantities.a .- Mquantities.b)/Mquantities.b^2 .+ 1)*model.D[r,r,i])[1],Squantities.eta)
            Mquantities.omega = Mquantities.psitilde[r]/model.D[r,r,i] - 1
            if r<Squantities.p && Mquantities.omega!=0
                Mquantities.ratio = Mquantities.omega/(1 + Mquantities.omega*Mquantities.b)
                @views Mquantities.C[1:r, (r+1):Squantities.p] = Mquantities.C[1:r, (r+1):Squantities.p] .- Mquantities.ratio .* Mquantities.c*Array(Mquantities.bvec[(r+1):Squantities.p]')
            end
        end
        @views model.D[:,:,i] = Diagonal(Mquantities.psitilde[:,1])
    end
    if Squantities.update_q
        if Mquantities.qup.q_change
            model.B = Mquantities.qup.B_update
        end
    end #problem is here...
    #if Mquantities.i_break
    #    return "failedModel"
    #else
    #    return model
    #end
end

function MStep_msmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)
    #update = MFAModelM(model.pivec,model.mu,model.B,model.D,model.g,model.q,model.tau,model.logL,model.BIC)
    #p = size(Y,2); n = size(Y,1)
    #if p==2
    #    q_max = 1
    #else
    #    q_max = Int(floor(p + (1 - sqrt(1 + 8*p))/2))
    #end
    Mquantities.n_i = sum(model.tau, dims = 1)

    #CM-Step 1: Updating pi and mu_i
    model.pivec = Mquantities.n_i ./ Squantities.n
    #mu_new = Array{Float64,2}(undef,p,g)
    #Mquantities.tau_xi_ij = zeros(n, model.g)
    @inbounds for i in 1:model.g
        Mquantities.tau_xi_ij[:,i] = model.tau[:,i] .* model.w_esteps[:,1,i]
        @views model.mu[:,i] = sum(Squantities.Y .* Mquantities.tau_xi_ij[:,i], dims = 1) ./ sum(Mquantities.tau_xi_ij[:,i])
    end

    #CM-Step 2: Updating q
    #Vtilde = @MArray zeros(Float64,p,p,model.g)
    #Lambda = @MMatrix zeros(Float64,q_max,model.g)
    #Lambda_new = Array{Float64,2}(undef,q_max,model.g)
    #eigenvecs =  Array{Any,1}(undef,model.g)
    ecm_failed = false
    @inbounds for i in 1:model.g
        @views Mquantities.Ymu = Squantities.Y .- model.mu[:,i]'
        @views Mquantities.Vi = (Array(Mquantities.Ymu') * (Mquantities.Ymu .* Mquantities.tau_xi_ij[:,i]))/Mquantities.n_i[i]
        @views Mquantities.D_is = Diagonal(1 ./ sqrt.(model.D[:,:,i]))
        Mquantities.Vtilde[:,:,i] = Mquantities.D_is*Mquantities.Vi*Mquantities.D_is
        @views Mquantities.V_eig_temp = eigvals(Array(Mquantities.Vtilde[:,:,i]))
        #if any(imag.(V_eig_temp) .!= 0)
        #    print("uh oh")
        #    print(imag.(V_eig_temp))
        #    print("converting to real number >:(")
        #end
        #V_eig = real.(V_eig_temp)
         #end

        if all( Mquantities.V_eig_temp .< 1 )
            ecm_failed = true
            #println("Complex updating step used")
            #println(Mquantities.V_eig_temp)

        #else
            #println("...")
        end

        if any(imag(Mquantities.V_eig_temp) .> 0)
            println("Complex eigenvalues appeared!!")
            model.status = "failed"
            return
        end

        Mquantities.eig_order = sortperm(Mquantities.V_eig_temp, rev = true)
        @views Mquantities.Lambda[:,i] = Mquantities.V_eig_temp[Mquantities.eig_order][1:Squantities.qmax]
        @views Mquantities.eigenvecs[i] = real.(eigvecs(Array(Mquantities.Vtilde[:,:,i]))[:,Mquantities.eig_order])
        if any(real.(Mquantities.Lambda[:,i]) .< 0)
            Mquantities.Lambda[(real.(Mquantities.Lambda[:,i]) .< 0),i] .= 1e-20
            println("Negative eigenvalues appeared!")
        end
    end

    if (!ecm_failed | model.name != "MGHFA" )

        #CM-Step 3: Updating B_i, then D_i in same loop.
        Mquantities.qi_break = false
        @inbounds for i in 1:model.g
            #Update B
            @views Mquantities.qi = sum(Mquantities.Lambda[:,i] .> 1)
            if Mquantities.qi > model.q
                Mquantities.qi = model.q
            end

            if Mquantities.qi == 0
                #qi_break = true
                Mquantities.qi = argmin(cumsum(log.(Mquantities.Lambda[:,i]) .- Mquantities.Lambda[:,i] .+ 1))
                #println("No eigenvalues greater than 1")
            break
            end
            @views Mquantities.Uqi = Mquantities.eigenvecs[i][:,1:Mquantities.qi]
            @views Mquantities.Ri = I(model.q)[1:Mquantities.qi,:]
            @views model.B[:,:,i] = sqrt.(model.D[:,:,i]) * Mquantities.Uqi * Diagonal(sqrt.(Complex.(Mquantities.Lambda[1:Mquantities.qi,i] - ones(1,Mquantities.qi)))) * Mquantities.Ri

            #Update D
            Mquantities.C = Array(1.0*I(Squantities.p)); Mquantities.Us = Mquantities.Uqi*(Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi,i])-I(Mquantities.qi)); Mquantities.V = Array(Mquantities.Uqi')
            Mquantities.bvec = zeros(Squantities.p,1)
            Mquantities.psitilde = Array{Float64,2}(undef,Squantities.p,1)
            @inbounds for r in 1:Squantities.p
                @views Mquantities.c = Mquantities.C[1:r,r]
                @views Mquantities.v = Mquantities.V[:,1:r]*Mquantities.c
                @views Mquantities.bvec[r:Squantities.p] = Mquantities.Us[r:Squantities.p,:]*Mquantities.v
                Mquantities.b = Mquantities.bvec[r] + 1
                @views Mquantities.a = Array(Mquantities.v')*Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi, i] .- Mquantities.Lambda[1:Mquantities.qi,i])*Mquantities.v + Array(Mquantities.c')*Mquantities.Vtilde[1:r,1:r,i] * Mquantities.c
                Mquantities.psitilde[r] = max((((Mquantities.a .- Mquantities.b)/Mquantities.b^2 .+ 1)*model.D[r,r,i])[1],Squantities.eta)
                Mquantities.omega = Mquantities.psitilde[r]/model.D[r,r,i] - 1
                if r<Squantities.p && Mquantities.omega!=0
                    Mquantities.ratio = Mquantities.omega/(1 + Mquantities.omega*Mquantities.b)
                    @views Mquantities.C[1:r, (r+1):Squantities.p] = Mquantities.C[1:r, (r+1):Squantities.p] .- Mquantities.ratio .* Mquantities.c*Array(Mquantities.bvec[(r+1):Squantities.p]')
                end
            end
            @views model.D[:,:,i] = Diagonal(Mquantities.psitilde[:,1])

            if any(isnan.(model.D[:,:,i]))
                println("Not good...")
            end

        end
    else
        logL_mmvmnfa(model, Mquantities, Squantities)
        eu3 = zeros(Squantities.p,Squantities.p,model.g)
        @inbounds for i in 1:model.g
            Ymu = Squantities.Y .- model.mu[:,i]'
            bi = model.B[:,:,i]' * inv(model.B[:,:,i]*model.B[:,:,i]' + model.D[:,:,i])
            eu1 = bi * ( ( model.w_esteps[:, 1, i] .* Ymu) .- model.beta[:,i]' )
            eu2 = bi * ( model.w_esteps[:, 2, i] .* Ymu  .- model.beta[:,i]')
            @inbounds for j in 1:Squantities.n
                eu3[:,:,j] = model.w_esteps[j, 2, i]*( I(model.q) - bi*model.B[:,:,i] + bi*Ymu[j,:]'*Ymu[j,:]*bi') - bi*(Ymu[j,:]'*model.beta[:,:,i] + model.beta[:,:,i]' * Ymu[j,:])*bi' + model.w_esteps[j, 1, i] * bi * model.beta[:,:,i]' * model.beta[:,:,i] * bi'
            end

            #V = (Array(Ymu') * (Ymu .* model.tau[:,i]))/Mquantities.n_i[i]
            #gamma = inv(real.(model.B[:,:,i]*model.B[:,:,i]' + model.D[:,:,i]))*model.B[:,:,i]
            #omega = I(model.q) - gamma' * model.B[:,:,i]
            #model.B[:,:,i] =  V*gamma*inv(gamma'*V*gamma + omega)
            #model.D[:,:,i] = diagm(diag(V - V*gamma*model.B[:,:,i]'))
        end

    end
    #CM Step 4: Updating Psi
    @inbounds for i in 1:model.g
            @views model.Ψ[:,i] = Squantities.msteps(model.tau[:,i], model.w_esteps[:,:,i], Mquantities)
    end

end

function MStep_mmvmnfa(model::MMVMNFAModel, Mquantities::MMVMNFAMutableQuantities, Squantities::MMVMNFAStaticQuantities)
    #update = MFAModelM(model.pivec,model.mu,model.B,model.D,model.g,model.q,model.tau,model.logL,model.BIC)
    #p = size(Y,2); n = size(Y,1)
    #if p==2
    #    q_max = 1
    #else
    #    q_max = Int(floor(p + (1 - sqrt(1 + 8*p))/2))
    #end
    Mquantities.n_i = sum(model.tau, dims = 1)

    #CM-Step 1: Updating pi and mu_i
    model.pivec = Mquantities.n_i ./ Squantities.n
    #mu_new = Array{Float64,2}(undef,p,g)
    #tau_rho_xi_ij = zeros(n, model.g); tau_rho_ij = zeros(n, model.g)
    #tr_ij = zeros(n, model.g); tx_ij = zeros(n, model.g)
    @inbounds for i in 1:model.g
        @views Mquantities.xi_bar_i = sum(model.tau[:,i] .* model.w_esteps[:,1,i])/sum(model.tau[:,i])
        @views Mquantities.rho_bar_i = sum(model.tau[:,i] .* model.w_esteps[:,2,i])/sum(model.tau[:,i])
        @views Mquantities.tau_rho_xi_ij[:,i] = model.tau[:,i] .* (Mquantities.xi_bar_i .* model.w_esteps[:,2,i] .- 1)
        @views Mquantities.tau_rho_ij[:,i] = model.tau[:,i] .* (Mquantities.rho_bar_i .- model.w_esteps[:,2,i])
        @views Mquantities.tr_ij[:,i] = model.tau[:,i] .* model.w_esteps[:,2,i]
        @views Mquantities.tx_ij[:,i] = model.tau[:,i] .* model.w_esteps[:,1,i]
        @views model.mu[:,i] = sum(Squantities.Y .* Mquantities.tau_rho_xi_ij[:,i], dims = 1) ./ sum(Mquantities.tau_rho_xi_ij[:,i])
        @views model.beta[:,i] = sum(Squantities.Y .* Mquantities.tau_rho_ij[:,i], dims = 1) ./ sum(Mquantities.tau_rho_xi_ij[:,i])
        #if any(sum(Y .== model.mu[:,i]', dims = 2) == p)
        #    println("Data point has dominated and become the mean of a component.")
        #    break
        #end
    end

    #CM-Step 2: Updating q
    #Vtilde = @MArray zeros(Float64,p,p,model.g)
    #Lambda = @MMatrix zeros(Float64,q_max,model.g)
    #Lambda_new = Array{Float64,2}(undef,q_max,model.g)
    #eigenvecs =  Array{Any,1}(undef,model.g)
    ecm_failed = false
        @inbounds for i in 1:model.g
            Mquantities.beta_matrix = reshape(repeat(model.beta[:,i],inner=1,outer=Squantities.n), Squantities.p, Squantities.n) |> permutedims
            @views Mquantities.Ymu = Squantities.Y .- model.mu[:,i]'
            @views Mquantities.Vi = (Array(Mquantities.Ymu') * (Mquantities.Ymu .* Mquantities.tr_ij[:,i]) -
                        Array(Mquantities.Ymu') * (Mquantities.beta_matrix .* model.tau[:,i]) -
                        Array(Mquantities.beta_matrix') * (Mquantities.Ymu .* model.tau[:,i]) +
                        Array(Mquantities.beta_matrix') * (Mquantities.beta_matrix .* Mquantities.tx_ij[:,i]))/Mquantities.n_i[i]
            @views Mquantities.D_is = Diagonal(1 ./ sqrt.(model.D[:,:,i]))
            Mquantities.Vtilde[:,:,i] = Mquantities.D_is*Mquantities.Vi*Mquantities.D_is
            @views Mquantities.V_eig_temp = eigvals(Array(Mquantities.Vtilde[:,:,i]))
            #if any(imag.(V_eig_temp) .!= 0)
            #    print("uh oh")
            #    print(imag.(V_eig_temp))
            #    print("converting to real number >:(")
            #end
            if all( Mquantities.V_eig_temp .< 1 )
                ecm_failed = true
                #println("Complex updating step used")
                #println(Mquantities.V_eig_temp)

            #else
                #println("...")
            end

            if any(imag(Mquantities.V_eig_temp) .> 0)
                println("Complex eigenvalues appeared!!")
                model.status = "failed"
                return
            end
            #V_eig = real.(V_eig_temp)

            Mquantities.eig_order = sortperm(Mquantities.V_eig_temp, rev = true)

            @views Mquantities.Lambda[:,i] = Mquantities.V_eig_temp[Mquantities.eig_order][1:Squantities.qmax]
            @views Mquantities.eigenvecs[i] = eigvecs(Array(Mquantities.Vtilde[:,:,i]))[:,Mquantities.eig_order]
            if any(Mquantities.Lambda[:,i] .< 0)
                println("Negative eigenvalues appeared!")
                Mquantities.Lambda[(Mquantities.Lambda[:,i] .< 0),i] .= 1e-20
            end
        end

        if true #(!ecm_failed | model.name != "MGHFA")
        #println("ECM update used")
        #CM-Step 3: Updating B_i, then D_i in same loop.
        Mquantities.qi_break = false
        @inbounds for i in 1:model.g
            #Update B
            @views Mquantities.qi = sum(Mquantities.Lambda[:,i] .> 1)
            if Mquantities.qi > model.q
                Mquantities.qi = model.q
            end

            if Mquantities.qi == 0
                Mquantities.qi = argmin(cumsum(log.(Mquantities.Lambda[:,i]) .- Mquantities.Lambda[:,i] .+ 1))
                #qi_break = true
                #println("No eigenvalues greater than 1")
            #break
            end
            @views Mquantities.Uqi = Mquantities.eigenvecs[i][:,1:Mquantities.qi]
            @views Mquantities.Ri = I(model.q)[1:Mquantities.qi,:]
            @views model.B[:,:,i] = sqrt.(model.D[:,:,i]) * Mquantities.Uqi * Diagonal(sqrt.(Complex.(Mquantities.Lambda[1:Mquantities.qi,i:i]' - ones(1,Mquantities.qi)))) * Mquantities.Ri

            #Update D
            Mquantities.C = Array(1.0*I(Squantities.p)); Mquantities.Us = Mquantities.Uqi*(Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi,i])-I(Mquantities.qi)); Mquantities.V = Array(Mquantities.Uqi')
            Mquantities.bvec = zeros(Squantities.p,1)
            Mquantities.psitilde = Array{Float64,2}(undef,Squantities.p,1)
            #Mquantities.bvec = zeros(p,1)
        # Mquantities.psitilde = Array{Float64,2}(undef,p,1) #Debugging up to here ...
            @inbounds for r in 1:Squantities.p
                @views Mquantities.c = Mquantities.C[1:r,r]
                @views Mquantities.v = Mquantities.V[:,1:r]*Mquantities.c
                @views Mquantities.bvec[r:Squantities.p] = Mquantities.Us[r:Squantities.p,:]*Mquantities.v
                Mquantities.b = Mquantities.bvec[r] + 1
                @views Mquantities.a = Array(Mquantities.v')*Diagonal(1 ./ Mquantities.Lambda[1:Mquantities.qi, i] .- Mquantities.Lambda[1:Mquantities.qi,i])*Mquantities.v + Array(Mquantities.c')*Mquantities.Vtilde[1:r,1:r,i] * Mquantities.c
                Mquantities.psitilde[r] = max((((Mquantities.a .- Mquantities.b)/Mquantities.b^2 .+ 1)*model.D[r,r,i])[1],Squantities.eta)
                Mquantities.omega = Mquantities.psitilde[r]/model.D[r,r,i] - 1
                if r<Squantities.p && Mquantities.omega!=0
                    Mquantities.ratio = Mquantities.omega/(1 + Mquantities.omega*Mquantities.b)
                    @views Mquantities.C[1:r, (r+1):Squantities.p] = Mquantities.C[1:r, (r+1):Squantities.p] .- Mquantities.ratio .* Mquantities.c*Array(Mquantities.bvec[(r+1):Squantities.p]')
                end
            end
            @views model.D[:,:,i] = Diagonal(Mquantities.psitilde[:,1])
        end
    else
        println("AECM update used")
        logL_mmvmnfa(model, Mquantities, Squantities)
        println(model.logL)
        eu3 = zeros(model.q,model.q,Squantities.n)
        @inbounds for i in 1:model.g
            Ymu = Squantities.Y .- model.mu[:,i]'
            bi = model.B[:,:,i]' * inv(model.B[:,:,i]*model.B[:,:,i]' + model.D[:,:,i])

            eu1 = (bi * ( Ymu .- model.w_esteps[:, 1, i].*model.beta[:,i]' )')'
            eu2 = (bi * ( model.w_esteps[:, 2, i] .* Ymu  .- model.beta[:,i]')')'
            @inbounds for j in 1:Squantities.n
                eu3[:,:,j] = model.w_esteps[j, 2, i]*( I(model.q) - bi*model.B[:,:,i] + bi*Ymu[j:j,:]'  *Ymu[j:j,:]*bi') - bi*(Ymu[j:j,:]'*model.beta[:,i:i]' + model.beta[:,i:i] * Ymu[j:j,:])*bi' + model.w_esteps[j, 1, i] * bi * model.beta[:,i:i] * model.beta[:,i:i]' * bi'
            end

            #*inv(sum([model.tau[j,i] .* eu3[:,:,j] for j in 1:Squantities.n],dims = 3))

            model.B[:,:,i] = sum([model.tau[j,i] * (Ymu[j:j,:]'*eu2[j:j,:] .- model.beta[:,i:i]*eu1[j:j,:]) for j in 1:Squantities.n]  ) * inv(sum([model.tau[j,i] .* eu3[:,:,j] for j in 1:Squantities.n], dims = 1)[1])

            logL_mmvmnfa(model, Mquantities, Squantities)

            #update for B not quite right??

            #println(model.logL)
            model.D[:,:,i] = 1/(sum(model.tau[:,i]))*diagm(diag(sum(
            [ model.tau[j,i] * (
            model.w_esteps[j,2,i]*Ymu[j:j,:]'*Ymu[j:j,:] - 2*model.beta[:,i:i]*Ymu[j:j,:] + model.w_esteps[j,1,i]*model.beta[:,i:i]*model.beta[:,i:i]' - 2*Ymu[j:j,:]'*eu2[j:j,:]*model.B[:,:,i]' + 2*model.beta[:,i:i]*eu1[j:j,:]*model.B[:,:,i]' + model.B[:,:,i]*eu3[:,:,j]*model.B[:,:,i]'
            )
            for j in 1:Squantities.n]
            ,dims=1)[1]))
            #logL_mmvmnfa(model, Mquantities, Squantities)
            #println(model.logL)
            #eu3[:,:,j] for j in 1,dims = 1)
            #V = (Array(Ymu') * (Ymu .* model.tau[:,i]))/Mquantities.n_i[i]
            #gamma = inv(real.(model.B[:,:,i]*model.B[:,:,i]' + model.D[:,:,i]))*model.B[:,:,i]
            #omega = I(model.q) - gamma' * model.B[:,:,i]
            #model.B[:,:,i] =  V*gamma*inv(gamma'*V*gamma + omega)
            #model.D[:,:,i] = diagm(diag(V - V*gamma*model.B[:,:,i]'))
        end

    end

    #CM Step 4: Updating Psi
    @inbounds for i in 1:model.g
        @views model.Ψ[:,i] = Squantities.msteps(model.tau[:,i], model.w_esteps[:,:,i], Mquantities)
    end

end

max_stars_repo_name: john-c-davey/FactorMixtures

max_issues_repo_name: john-c-davey/FactorMixtures

max_forks_repo_name: john-c-davey/FactorMixtures


text:
function JuliaAppGUIPrivate()
    println("It works!!")
end

function JuliaAppGUI()
    println(JuliaAppGUIPrivate())
end


text:
using BayesianStatisticsTheFunWay
using Test

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


text:
for i = 1:5
    print("\rRunning argparse_test$i")
    try
        include("argparse_test$i.jl")
    catch err
        println()
        rethrow(err)
    end
end
println()


text:
using ImplicitGlobalGrid, CUDA, Plots
#(...)

@views function diffusion3D()
    # Physics
    #(...)

    # Numerics
    #(...)
    me, dims   = init_global_grid(nx, ny, nz);              # Initialize the implicit global grid
    #(...)

    # Array initializations
    #(...)

    # Initial conditions (heat capacity and temperature with two Gaussian anomalies each)
    #(...)

    # Preparation of visualisation
    gr()
    ENV["GKSwstype"]="nul"
    anim = Animation();
    nx_v = (nx-2)*dims[1];
    ny_v = (ny-2)*dims[2];
    nz_v = (nz-2)*dims[3];
    T_v  = zeros(nx_v, ny_v, nz_v);
    T_nohalo = zeros(nx-2, ny-2, nz-2);

    # Time loop
    #(...)
    for it = 1:nt
        if mod(it, 1000) == 1                                                                 # Visualize only every 1000th time step
            T_nohalo .= T[2:end-1,2:end-1,2:end-1];                                           # Copy data to CPU removing the halo.
            gather!(T_nohalo, T_v)                                                            # Gather data on process 0 (could be interpolated/sampled first)
            if (me==0) heatmap(transpose(T_v[:,ny_v÷2,:]), aspect_ratio=1); frame(anim); end  # Visualize it on process 0.
        end
        #(...)
    end

    # Postprocessing
    if (me==0) gif(anim, "diffusion3D.gif", fps = 15) end                                     # Create a gif movie on process 0.
    if (me==0) mp4(anim, "diffusion3D.mp4", fps = 15) end                                     # Create a mp4 movie on process 0.
    finalize_global_grid();                                                                   # Finalize the implicit global grid
end

diffusion3D()


text:
module TestSlicedArray

include("preamble.jl")

const TEST_ALONGS = [
    (True(), ),
    (False(), ),

    (True(), True()),
    (True(), False()),
    (False(), True()),
    (False(), False()),
]

function makedata(V::Type, alongs::NTuple{L,TypedBool}) where {L}
    psize = testdims(L)
    insize = tuple_getindex(psize, alongs)
    outsize = tuple_getindex(psize, invert(alongs))
    M = length(insize)
    N = length(outsize)
    alongs_int = tuple_getindex(ntuple(identity, L), alongs)

    parent = rand!(Array{V,L}(undef, psize...))

    nested = Array{Array{V,M},N}(undef, outsize...)
    i = 0
    Base.mapslices(parent, dims = alongs_int) do el
        i += 1 # TODO rand!
        nested[i] = zeros(V, insize...)
        nested[i] .= el
        el
    end

    return (
        alongs_int = alongs_int,
        alongs_glob = map(a -> a === True() ? Colon() : *, alongs),
        insize = insize,
        outsize = outsize,
        M = M,
        N = N,
        nested = nested,
        parent = parent,
    )
end

showalongs(alongs) = "($(join(map(a -> a === True() ? ':' : '*', alongs), ", ")))"

@testset "alongs = $(showalongs(alongs)), V = $V" for alongs in TEST_ALONGS, V in (Float64,)
    @testset "constructors" begin
        @unpack alongs_int, alongs_glob, parent, M, N = makedata(V, alongs)
        Expected = SlicedArray{<:AbstractArray{V,M},N,M,Array{V,M+N},typeof(alongs)}

        @test slice(parent, alongs) isa Expected
        @test_inferred slice(parent, alongs)
        @test slice(parent, alongs...) isa Expected
        @test_inferred slice(parent, alongs...)

        @test slice(parent, alongs_glob) isa Expected
        @test_inferred slice(parent, alongs_glob)
        @test slice(parent, alongs_glob...) isa Expected
        @test_inferred slice(parent, alongs_glob...)

        @test slice(parent, alongs_int) isa Expected
        @test slice(parent, alongs_int...) isa Expected
    end

    @testset "IndexStyle" begin
        S = slice(makedata(V, alongs).parent, alongs)
        if ndims(S) == 1
            @test IndexStyle(typeof(S)) === IndexLinear()
        else
            @test IndexStyle(typeof(S)) === IndexCartesian()
        end
    end

    let V = V, alongs = alongs
        test_array_AB() do
            data = makedata(V, alongs)
            A = SlicedArray(data.parent, alongs)
            B = data.nested
            return A, B
        end
    end

    @testset "mapslices ($name)" for (name, f) in (
        ("identity", identity),
        ("sum", el -> sum(el)),
        ("reshape(..., reverse(dims))", el -> el isa AbstractArray ? reshape(el, reverse(size(el))) : el),
        ("reshape(..., Val(1))", el -> el isa AbstractArray ? reshape(el, Val(1)) : el),
    )
        data = makedata(V, alongs)
        B1 = mapslices(f, data.parent, dims = data.alongs_int)
        B2 = SpecialArrays.mapslices(f, data.parent, dims = data.alongs_int)
        @test B1 == B2
        @test_inferred SpecialArrays.mapslices(f, data.parent, dims = alongs)
    end

    @testset "align" begin
        data = makedata(V, alongs)
        @test SpecialArrays.align(slice(data.parent, alongs), alongs) === data.parent
        @test SpecialArrays.align(data.nested, alongs) == data.parent
    end

    @testset "flatten" begin
        @unpack parent, nested = makedata(V, alongs)
        S = slice(parent, alongs)
        if SpecialArrays.iscontiguous(alongs)
            @test flatten(S) === S.parent === parent
        else
            flat = reshape(reduce(hcat, nested), (inner_size(nested)..., size(nested)...))
            @test flatten(S) == flat
        end
    end

    @testset "inner_*" begin
        data = makedata(V, alongs)
        S = SlicedArray(data.parent, alongs)

        @test inner_size(S) == size(first(S))
        @test_inferred inner_size(S)
        @test_noalloc inner_size($S)

        @test inner_axes(S) == axes(first(S))
        @test_inferred inner_axes(S)
        @test_noalloc inner_axes($S)
    end
end

@testset "ContiguousSlicedVector ($(showalongs(alongs))" for alongs in (
    (False(),),
    (True(), False()),
    (True(), True(), False()),
)
    V = Float64

    # TODO we define resize!/sizehint! for all SlicedArrays, not just ContiguousSlicedVector,
    # but ElasticArray only supports resizing along the last dimension, so that's all we
    # (currently) test for.
    @testset "resize!/empty!" begin
        data = makedata(V, alongs)
        S = slice(ElasticArray(data.parent), alongs)
        len = length(S)
        @test empty!(S) === S
        @test length(S) == 0
        @test resize!(S, len + 1) === S
        @test length(S) == len + 1
    end
    @testset "sizehint!" begin
        data = makedata(V, alongs)
        S = slice(ElasticArray(data.parent), alongs)
        @test sizehint!(S, length(S) + 1) === S
    end

    @testset "append!" begin
        data = makedata(V, alongs)
        S = empty!(slice(ElasticArray(data.parent), alongs))
        src1 = slice(rand!(ElasticArray(copy(data.parent))), alongs)
        src2 = slice(rand!(ElasticArray(copy(data.parent))), alongs)
        @test append!(S, src1) === S
        @test append!(S, src2) === S
        l = length(src1)
        @test S[1:l] == src1
        @test S[l+1:end] == src2
    end
    @testset "prepend!" begin
        data = makedata(V, alongs)
        S = empty!(slice(ElasticArray(data.parent), alongs))
        src1 = slice(rand!(ElasticArray(copy(data.parent))), alongs)
        src2 = slice(rand!(ElasticArray(copy(data.parent))), alongs)
        @test prepend!(S, src1) === S
        @test prepend!(S, src2) === S
        l = length(src1)
        @test S[1:l] == src2
        @test S[l+1:end] == src1
    end
    @testset "pop!" begin
        data = makedata(V, alongs)
        S = slice(ElasticArray(data.parent), alongs)
        xs = Array{V,data.M}[]
        @test all(1:length(S)) do _
            x = pop!(S)
            pushfirst!(xs, x)
            !(x isa SubArray)
        end
        @test xs == data.nested
    end
    @testset "push!" begin
        data = makedata(V, alongs)
        S = slice(ElasticArray(data.parent), alongs)
        empty!(S)
        @test all(data.nested) do x
            push!(S, x) === S
        end
        @test S == data.nested
    end
    @testset "pushfirst!" begin
        data = makedata(V, alongs)
        data.parent[:] .= 1:length(data.parent)
        S = slice(ElasticArray(data.parent), alongs)
        empty!(S)
        @test all(data.nested) do x
            pushfirst!(S, x) === S
        end
        @test S == reverse(data.nested)
    end
end

@testset "parent/dataids" begin
    S = slice(rand(2, 3, 4), 1, 3)
    @test parent(S) === S.parent
    @test Base.dataids(S) === Base.dataids(S.parent)
end

@testset "UnsafeArrays" begin
    S = slice(rand(2, 3, 4), 1, 3)
    Sv = uview(S)
    @test parent(Sv) isa UnsafeArray{eltype(S.parent),ndims(S.parent)}
    @test S == Sv
end

@testset "slice(::SlicedArray, ...)" begin
    A = rand(2, 3, 4)
    S1 = slice(A, *, :, :)
    S2 = slice(slice(A, *, :, *), *, :)
    @test S1 == S2
    @test S1.parent === S2.parent
end

# TODO Adapt.adapt_storage

end # module


text:
module CLITest

using ErlPort.CLI

function run()::Nothing
    return nothing
end # run()

end # module


text:
struct SerialVariableTargetAdjacency{T}
    """
    Start indices of respective column data in `colentries`
    """
    colstart::Vector{T}
end

"""
    $(SIGNATURES)

    Comparison of two adjacencies
"""
function Base.:(==)(a::SerialVariableTargetAdjacency{Ta}, b::SerialVariableTargetAdjacency{Tb}) where {Ta,Tb}
    Ta==Tb && a.colstart==b.colstart
end

"""
$(TYPEDSIGNATURES)
Create an empty SerialVariableTargetAdjacency
"""
SerialVariableTargetAdjacency(t::Type{T}) where T=SerialVariableTargetAdjacency{T}([one(T)])

"""
$(TYPEDSIGNATURES)
Create an empty SerialVariableTargetAdjacency with default type
"""
SerialVariableTargetAdjacency()=SerialVariableTargetAdjacency(Int64)

"""
$(TYPEDSIGNATURES)
Show adjacency (in trasposed form; preliminary)
"""
function Base.show(io::IO,adj::SerialVariableTargetAdjacency)
    for isource=1:num_sources(adj)
        for itarget=1:num_targets(adj,isource)
            print(adj[itarget,isource], " ")
        end
        println()
    end
end

"""
$(TYPEDSIGNATURES)
Access adjacency as if it is a 2D Array
"""
Base.getindex(adj::SerialVariableTargetAdjacency,i,isource)=adj.colstart[isource]+i-1
Base.getindex(adj::SerialVariableTargetAdjacency,::Colon,isource)=adj.colstart[isource]:adj.colstart[isource+1]-1

"""
$(TYPEDSIGNATURES)
Number of targets for given source
"""
ExtendableGrids.num_targets(adj::SerialVariableTargetAdjacency,isource)=adj.colstart[isource+1]-adj.colstart[isource]

"""
$(TYPEDSIGNATURES)
Number of sources in adjacency
"""
ExtendableGrids.num_sources(adj::SerialVariableTargetAdjacency)=length(adj.colstart)-1

"""
$(TYPEDSIGNATURES)
Maximum number of targets per source
"""
ExtendableGrids.max_num_targets_per_source(adj::SerialVariableTargetAdjacency)=maximum(adj.colstart[2:end].-adj.colstart[1:end-1])

"""
$(TYPEDSIGNATURES)
Append a column to adjacency.
"""
function Base.append!(adj::SerialVariableTargetAdjacency,len)
    push!(adj.colstart,adj.colstart[end]+len)
end


text:
# 2D Grid examples
# ===============
#

# ## Rectangle
function rectangle()
    X=collect(0:0.05:1)
    Y=collect(0:0.05:1)
    simplexgrid(X,X)
end
# ![](rectangle.svg)
#
# ## Rectangle with local refinement
#
function rectangle_localref()
    hmin=0.01
    hmax=0.1
    XLeft=geomspace(0.0,0.5,hmax,hmin)
    XRight=geomspace(0.5,1.0,hmin,hmax)
    X=glue(XLeft, XRight)
    simplexgrid(X,X)
end
# ![](rectangle_localref.svg)

#
# ## Rectangle with multiple regions
#
function rectangle_multiregion()
    X=collect(0:0.05:1)
    Y=collect(0:0.05:1)
    grid=simplexgrid(X,Y)
    cellmask!(grid,[0.0,0.0],[1.0,0.5],3)
    bfacemask!(grid,[0.0,0.0],[0.0,0.5],5)
    bfacemask!(grid,[1.0,0.0],[1.0,0.5],6)
    bfacemask!(grid,[0.0,0.5],[1.0,0.5],7)
end
# ![](rectangle_multiregion.svg)

#
# ## Subgrid from rectangle
#
function rectangle_subgrid()
    X=collect(0:0.05:1)
    Y=collect(0:0.05:1)
    grid=simplexgrid(X,Y)
    cellmask!(grid,[0.25,0.25],[0.75,0.75],2)
    bfacemask!(grid,[0.25,0.25],[0.25,0.75],5)
    bfacemask!(grid,[0.25,0.25],[0.75,0.25],5)
    bfacemask!(grid,[0.25,0.75],[0.75,0.75],5)
    bfacemask!(grid,[0.75,0.25],[0.75,0.75],5)
    subgrid(grid,[1])
end
# ![](rectangle_subgrid.svg)


text:
r = PH.ScalarPenaltyParameter(25.0)
atol = 1e-8
rtol = 1e-12
max_iter = 500
obj_val = 178.3537498401004
var_vals = Dict([
    "x[1]_{0,1,2,3}" => 7.5625,
    "x[2]_{0,1,2,3}" => 0.0,
    "x[3]_{0,1,2,3}" => 1.0,
    "y_{0,1}" => 1.75,
    "y_{2,3}" => 0.0,
    "z[1]_{0}" => -0.65625,
    "z[2]_{0}" => -0.65625,
    "z[1]_{1}" => 2.84375,
    "z[2]_{1}" => 2.84375,
    "z[1]_{2}" => -1.78125,
    "z[2]_{2}" => -1.78125,
    "z[1]_{3}" => 4.71875,
    "z[2]_{3}" => 4.71875,
])

# Setup distributed stuff
nw = 2
diff = nw - nworkers()
diff > 0 && Distributed.addprocs(nw)
@assert Distributed.nworkers() == nw
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)
@everywhere using ProgressiveHedging
@everywhere const PH = ProgressiveHedging
@everywhere using Ipopt
@everywhere using JuMP
include("common.jl")

@testset "Error Handling" begin
    ch_size = 10
    worker_inf = PH._launch_workers(ch_size, ch_size)

    @everywhere struct MakeWorkerError <: PH.Message end

    PH._send_message(worker_inf, first(workers()), MakeWorkerError())

    my_task = @async begin
        PH._wait_for_shutdown(worker_inf)
    end
    if timeout_wait(my_task)
        @test istaskfailed(my_task)
        @test typeof(my_task.exception) <: RemoteException
    else
        error("Timed out on test")
    end

    my_task = @async begin
        PH._wait_for_shutdown(worker_inf)
        return PH._is_running(worker_inf)
    end
    if timeout_wait(my_task)
        @test !fetch(my_task)
    else
        error("Timed out on test")
    end
end

@testset "Solve" begin

    # Solve the problem
    (n, err, rerr, obj, soln, phd) = PH.solve(build_scen_tree(),
                                              create_model,
                                              r,
                                              atol=atol,
                                              rtol=rtol,
                                              opt=Ipopt.Optimizer,
                                              opt_args=(print_level=0,tol=1e-12),
                                              max_iter=max_iter,
                                              report=0,
                                              timing=false,
                                              warm_start=true)

    @test err < atol
    @test isapprox(obj, obj_val)
    @test n < max_iter

    for row in eachrow(soln)
        var = row[:variable] * "_{" * row[:scenarios] * "}"
        @test isapprox(row[:value], var_vals[var], atol=1e-7)
    end

end

@testset "User Specified Workers" begin

    st = build_scen_tree()
    worker_assignments = Dict(k=>Set{PH.ScenarioID}() for k in 2:3)
    for s in PH.scenarios(st)
        if PH.value(s) < 3
            push!(worker_assignments[2], s)
        else
            push!(worker_assignments[3], s)
        end
    end

    # Solve the problem
    (n, err, rerr, obj, soln, phd) = PH.solve(build_scen_tree(),
                                              create_model,
                                              r;
                                              atol=atol,
                                              rtol=rtol,
                                              opt=Ipopt.Optimizer,
                                              opt_args=(print_level=0,tol=1e-12),
                                              max_iter=max_iter,
                                              report=0,
                                              worker_assignments=worker_assignments,
                                              timing=false,
                                              warm_start=true)

    @test err < atol
    @test isapprox(obj, obj_val)
    @test n < max_iter

    for row in eachrow(soln)
        var = row[:variable] * "_{" * row[:scenarios] * "}"
        @test isapprox(row[:value], var_vals[var], atol=1e-7)
    end

    for s in PH.scenarios(phd)
        if PH.value(s) < 3
            @test phd.scenario_map[s].pid == 2
        else
            @test phd.scenario_map[s].pid == 3
        end
    end

end

@testset "Lower Bound Algorithm" begin
    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ScalarPenaltyParameter(2.0),
                                              atol=atol,
                                              rtol=rtol,
                                              max_iter=max_iter,
                                              report=0,
                                              lower_bound=5,
                                              timing=false,
                                              warm_start=false
                                              )

    @test err < atol
    @test isapprox(obj, 8.25, atol=1e-6)
    @test n < max_iter

    lb_df = PH.lower_bounds(phd)
    for row in eachrow(lb_df)
        @test row[:bound] <= 8.25
    end

    @test size(lb_df,1) == 12
    @test isapprox(lb_df[size(lb_df,1), "absolute gap"], 0.0, atol=1e-7)
    @test isapprox(lb_df[size(lb_df,1), "relative gap"], 0.0, atol=(1e-7/8.25))
end

# Tear down distrbuted stuff
Distributed.rmprocs(Distributed.workers())
@assert Distributed.nprocs() == 1


text:
#==============================================================================#
# Pricing.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/pricing-2017-10-15.normal.json
#==============================================================================#

__precompile__()

module Pricing

using AWSCore

"""
    using AWSSDK.Pricing.describe_services
    describe_services([::AWSConfig], arguments::Dict)
    describe_services([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.pricing
    pricing([::AWSConfig], "DescribeServices", arguments::Dict)
    pricing([::AWSConfig], "DescribeServices", <keyword arguments>)

# DescribeServices Operation

Returns the metadata for one service or a list of the metadata for all services. Use this without a service code to get the service codes for all services. Use it with a service code, such as `AmazonEC2`, to get information specific to that service, such as the attribute names available for that service. For example, some of the attribute names available for EC2 are `volumeType`, `maxIopsVolume`, `operation`, `locationType`, and `instanceCapacity10xlarge`.

# Arguments

## `ServiceCode = ::String`
The code for the service whose information you want to retrieve, such as `AmazonEC2`. You can use the `ServiceCode` to filter the results in a `GetProducts` call. To retrieve a list of all services, leave this blank.

## `FormatVersion = ::String`
The format version that you want the response to be in.

Valid values are: `aws_v1`

## `NextToken = ::String`
The pagination token that indicates the next set of results that you want to retrieve.

## `MaxResults = ::Int`
The maximum number of results that you want returned in the response.

# Returns

`DescribeServicesResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

# Example: To retrieve a list of services and service codes

Input:
```
[
    "FormatVersion" => "aws_v1",
    "MaxResults" => 1,
    "ServiceCode" => "AmazonEC2"
]
```

Output:
```
Dict(
    "FormatVersion" => "aws_v1",
    "NextToken" => "abcdefg123",
    "Services" => [
        Dict(
            "AttributeNames" => [
                "volumeType",
                "maxIopsvolume",
                "instanceCapacity10xlarge",
                "locationType",
                "operation"
            ],
            "ServiceCode" => "AmazonEC2"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/pricing-2017-10-15/DescribeServices)
"""
@inline describe_services(aws::AWSConfig=default_aws_config(); args...) = describe_services(aws, args)

@inline describe_services(aws::AWSConfig, args) = AWSCore.Services.pricing(aws, "DescribeServices", args)

@inline describe_services(args) = describe_services(default_aws_config(), args)

"""
    using AWSSDK.Pricing.get_attribute_values
    get_attribute_values([::AWSConfig], arguments::Dict)
    get_attribute_values([::AWSConfig]; ServiceCode=, AttributeName=, <keyword arguments>)

    using AWSCore.Services.pricing
    pricing([::AWSConfig], "GetAttributeValues", arguments::Dict)
    pricing([::AWSConfig], "GetAttributeValues", ServiceCode=, AttributeName=, <keyword arguments>)

# GetAttributeValues Operation

Returns a list of attribute values. Attibutes are similar to the details in a Price List API offer file. For a list of available attributes, see [Offer File Definitions](http://docs.aws.amazon.com/awsaccountbilling/latest/aboutv2/reading-an-offer.html#pps-defs) in the [AWS Billing and Cost Management User Guide](http://docs.aws.amazon.com/awsaccountbilling/latest/aboutv2/billing-what-is.html).

# Arguments

## `ServiceCode = ::String` -- *Required*
The service code for the service whose attributes you want to retrieve. For example, if you want the retrieve an EC2 attribute, use `AmazonEC2`.

## `AttributeName = ::String` -- *Required*
The name of the attribute that you want to retrieve the values for, such as `volumeType`.

## `NextToken = ::String`
The pagination token that indicates the next set of results that you want to retrieve.

## `MaxResults = ::Int`
The maximum number of results to return in response.

# Returns

`GetAttributeValuesResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

# Example: To retrieve a list of attribute values

This operation returns a list of values available for the given attribute.

Input:
```
[
    "AttributeName" => "volumeType",
    "MaxResults" => 2,
    "ServiceCode" => "AmazonEC2"
]
```

Output:
```
Dict(
    "AttributeValues" => [
        Dict(
            "Value" => "Throughput Optimized HDD"
        ),
        Dict(
            "Value" => "Provisioned IOPS"
        )
    ],
    "NextToken" => "GpgauEXAMPLEezucl5LV0w==:7GzYJ0nw0DBTJ2J66EoTIIynE6O1uXwQtTRqioJzQadBnDVgHPzI1en4BUQnPCLpzeBk9RQQAWaFieA4+DapFAGLgk+Z/9/cTw9GldnPOHN98+FdmJP7wKU3QQpQ8MQr5KOeBkIsAqvAQYdL0DkL7tHwPtE5iCEByAmg9gcC/yBU1vAOsf7R3VaNN4M5jMDv3woSWqASSIlBVB6tgW78YL22KhssoItM/jWW+aP6Jqtq4mldxp/ct6DWAl+xLFwHU/CbketimPPXyqHF3/UXDw=="
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/pricing-2017-10-15/GetAttributeValues)
"""
@inline get_attribute_values(aws::AWSConfig=default_aws_config(); args...) = get_attribute_values(aws, args)

@inline get_attribute_values(aws::AWSConfig, args) = AWSCore.Services.pricing(aws, "GetAttributeValues", args)

@inline get_attribute_values(args) = get_attribute_values(default_aws_config(), args)

"""
    using AWSSDK.Pricing.get_products
    get_products([::AWSConfig], arguments::Dict)
    get_products([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.pricing
    pricing([::AWSConfig], "GetProducts", arguments::Dict)
    pricing([::AWSConfig], "GetProducts", <keyword arguments>)

# GetProducts Operation

Returns a list of all products that match the filter criteria.

# Arguments

## `ServiceCode = ::String`
The code for the service whose products you want to retrieve.

## `Filters = [[ ... ], ...]`
The list of filters that limit the returned products. only products that match all filters are returned.
```
 Filters = [[
        "Type" => <required> "TERM_MATCH",
        "Field" => <required> ::String,
        "Value" => <required> ::String
    ], ...]
```

## `FormatVersion = ::String`
The format version that you want the response to be in.

Valid values are: `aws_v1`

## `NextToken = ::String`
The pagination token that indicates the next set of results that you want to retrieve.

## `MaxResults = ::Int`
The maximum number of results to return in the response.

# Returns

`GetProductsResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

# Example: To retrieve a list of products

This operation returns a list of products that match the given criteria.

Input:
```
[
    "Filters" => [
        [
            "Field" => "ServiceCode",
            "Type" => "TERM_MATCH",
            "Value" => "AmazonEC2"
        ],
        [
            "Field" => "volumeType",
            "Type" => "TERM_MATCH",
            "Value" => "Provisioned IOPS"
        ]
    ],
    "FormatVersion" => "aws_v1",
    "MaxResults" => 1
]
```

Output:
```
Dict(
    "FormatVersion" => "aws_v1",
    "NextToken" => "57r3EXAMPLEujbzWfHF7Ciw==:ywSmZsD3mtpQmQLQ5XfOsIMkYybSj+vAT+kGmwMFq+K9DGmIoJkz7lunVeamiOPgthdWSO2a7YKojCO+zY4dJmuNl2QvbNhXs+AJ2Ufn7xGmJncNI2TsEuAsVCUfTAvAQNcwwamtk6XuZ4YdNnooV62FjkV3ZAn40d9+wAxV7+FImvhUHi/+f8afgZdGh2zPUlH8jlV9uUtj0oHp8+DhPUuHXh+WBII1E/aoKpPSm3c=",
    "PriceList" => [
        "{"product":{"productFamily":"Storage","attributes":{"storageMedia":"SSD-backed","maxThroughputvolume":"320 MB/sec","volumeType":"Provisioned IOPS","maxIopsvolume":"20000","servicecode":"AmazonEC2","usagetype":"CAN1-EBS:VolumeUsage.piops","locationType":"AWS Region","location":"Canada (Central)","servicename":"Amazon Elastic Compute Cloud","maxVolumeSize":"16 TiB","operation":""},"sku":"WQGC34PB2AWS8R4U"},"serviceCode":"AmazonEC2","terms":{"OnDemand":{"WQGC34PB2AWS8R4U.JRTCKXETXF":{"priceDimensions":{"WQGC34PB2AWS8R4U.JRTCKXETXF.6YS6EN2CT7":{"unit":"GB-Mo","endRange":"Inf","description":"\$0.138 per GB-month of Provisioned IOPS SSD (io1)  provisioned storage - Canada (Central)","appliesTo":[],"rateCode":"WQGC34PB2AWS8R4U.JRTCKXETXF.6YS6EN2CT7","beginRange":"0","pricePerUnit":{"USD":"0.1380000000"}}},"sku":"WQGC34PB2AWS8R4U","effectiveDate":"2017-08-01T00:00:00Z","offerTermCode":"JRTCKXETXF","termAttributes":{}}}},"version":"20170901182201","publicationDate":"2017-09-01T18:22:01Z"}"
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/pricing-2017-10-15/GetProducts)
"""
@inline get_products(aws::AWSConfig=default_aws_config(); args...) = get_products(aws, args)

@inline get_products(aws::AWSConfig, args) = AWSCore.Services.pricing(aws, "GetProducts", args)

@inline get_products(args) = get_products(default_aws_config(), args)

end # module Pricing

#==============================================================================#
# End of file
#==============================================================================#


text:
"""
    DGMulti(approximation_type::AbstractDerivativeOperator;
            element_type::AbstractElemShape,
            surface_flux=flux_central,
            surface_integral=SurfaceIntegralWeakForm(surface_flux),
            volume_integral=VolumeIntegralWeakForm(),
            kwargs...)

Create a summation by parts (SBP) discretization on the given `element_type`
using a tensor product structure based on the 1D SBP derivative operator
passed as `approximation_type`.

For more info, see the documentations of
[StartUpDG.jl](https://jlchan.github.io/StartUpDG.jl/dev/)
and
[SummationByPartsOperators.jl](https://ranocha.de/SummationByPartsOperators.jl/stable/).
"""
function DGMulti(approximation_type::AbstractDerivativeOperator;
                 element_type::AbstractElemShape,
                 surface_flux=flux_central,
                 surface_integral=SurfaceIntegralWeakForm(surface_flux),
                 volume_integral=VolumeIntegralWeakForm(),
                 kwargs...)

  rd = RefElemData(element_type, approximation_type; kwargs...)
  return DG(rd, nothing #= mortar =#, surface_integral, volume_integral)
end

function DGMulti(element_type::AbstractElemShape,
                 approximation_type::AbstractDerivativeOperator,
                 volume_integral,
                 surface_integral;
                 kwargs...)

  DGMulti(approximation_type, element_type=element_type,
          surface_integral=surface_integral, volume_integral=volume_integral)
end

function construct_1d_operators(D::AbstractDerivativeOperator, tol)
  nodes_1d = collect(grid(D))
  M = SummationByPartsOperators.mass_matrix(D)
  if M isa UniformScaling
    weights_1d = M * ones(Bool, length(nodes_1d))
  else
    weights_1d = diag(M)
  end

  # StartUpDG assumes nodes from -1 to +1. Thus, we need to re-scale everything.
  # We can adjust the grid spacing as follows.
  xmin = SummationByPartsOperators.xmin(D)
  xmax = SummationByPartsOperators.xmax(D)
  factor = 2 / (xmax - xmin)
  @. nodes_1d = factor * (nodes_1d - xmin) - 1
  @. weights_1d = factor * weights_1d

  D_1d = droptol!(inv(factor) * sparse(D), tol)
  I_1d = Diagonal(ones(Bool, length(nodes_1d)))

  return nodes_1d, weights_1d, D_1d, I_1d
end

function StartUpDG.RefElemData(element_type::Line,
                               D::AbstractDerivativeOperator;
                               tol = 100*eps())

  approximation_type = D
  N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree

  # 1D operators
  nodes_1d, weights_1d, D_1d = construct_1d_operators(D, tol)

  # volume
  rq = r = nodes_1d
  wq = weights_1d
  Dr = D_1d
  M = Diagonal(wq)
  Pq = LinearAlgebra.I
  Vq = LinearAlgebra.I

  VDM = nothing # unused generalized Vandermonde matrix

  rst = (r,)
  rstq = (rq,)
  Drst = (Dr,)

  # face
  face_vertices = StartUpDG.face_vertices(element_type)
  face_mask = [1, length(nodes_1d)]

  rf = [-1.0; 1.0]
  nrJ = [-1.0; 1.0]
  wf = [1.0; 1.0]
  if D isa AbstractPeriodicDerivativeOperator
    # we do not need any face stuff for periodic operators
    Vf = spzeros(length(wf), length(wq))
  else
    Vf = sparse([1, 2], [1, length(nodes_1d)], [1.0, 1.0])
  end
  LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))

  rstf = (rf,)
  nrstJ = (nrJ,)

  # low order interpolation nodes
  r1 = StartUpDG.nodes(element_type, 1)
  V1 = StartUpDG.vandermonde(element_type, 1, r) / StartUpDG.vandermonde(element_type, 1, r1)

  return RefElemData(
    element_type, approximation_type, N,
    face_vertices, V1,
    rst, VDM, face_mask,
    N, rst, LinearAlgebra.I, # plotting
    rstq, wq, Vq, # quadrature
    rstf, wf, Vf, nrstJ, # faces
    M, Pq, Drst, LIFT)
end

function StartUpDG.RefElemData(element_type::Quad,
                               D::AbstractDerivativeOperator;
                               tol = 100*eps())

  approximation_type = D
  N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree

  # 1D operators
  nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol)

  # volume
  s, r = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d)) # this is to match
                                                          # ordering of nrstJ
  rq = r; sq = s
  wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d))
  wq = wr .* ws
  Dr = kron(I_1d, D_1d)
  Ds = kron(D_1d, I_1d)
  M = Diagonal(wq)
  Pq = LinearAlgebra.I
  Vq = LinearAlgebra.I

  VDM = nothing # unused generalized Vandermonde matrix

  rst = (r, s)
  rstq = (rq, sq)
  Drst = (Dr, Ds)

  # face
  face_vertices = StartUpDG.face_vertices(element_type)
  face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s)...)

  rf, sf, wf, nrJ, nsJ = StartUpDG.init_face_data(element_type,
    quad_rule_face=(nodes_1d, weights_1d))
  if D isa AbstractPeriodicDerivativeOperator
    # we do not need any face stuff for periodic operators
    Vf = spzeros(length(wf), length(wq))
  else
    Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask)))
  end
  LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))

  rstf = (rf, sf)
  nrstJ = (nrJ, nsJ)

  # low order interpolation nodes
  r1, s1 = StartUpDG.nodes(element_type, 1)
  V1 = StartUpDG.vandermonde(element_type, 1, r, s) / StartUpDG.vandermonde(element_type, 1, r1, s1)

  return RefElemData(
    element_type, approximation_type, N,
    face_vertices, V1,
    rst, VDM, face_mask,
    N, rst, LinearAlgebra.I, # plotting
    rstq, wq, Vq, # quadrature
    rstf, wf, Vf, nrstJ, # faces
    M, Pq, Drst, LIFT)
end

function StartUpDG.RefElemData(element_type::Hex,
                               D::AbstractDerivativeOperator;
                               tol = 100*eps())

  approximation_type = D
  N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree

  # 1D operators
  nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol)

  # volume
  # to match ordering of nrstJ
  s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d))
  rq = r; sq = s; tq = t
  wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d))
  wq = wr .* ws .* wt
  Dr = kron(I_1d, I_1d, D_1d)
  Ds = kron(I_1d, D_1d, I_1d)
  Dt = kron(D_1d, I_1d, I_1d)
  M = Diagonal(wq)
  Pq = LinearAlgebra.I
  Vq = LinearAlgebra.I

  VDM = nothing # unused generalized Vandermonde matrix

  rst = (r, s, t)
  rstq = (rq, sq, tq)
  Drst = (Dr, Ds, Dt)

  # face
  face_vertices = StartUpDG.face_vertices(element_type)
  face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s, t)...)

  rf, sf, tf, wf, nrJ, nsJ, ntJ = let
    rf, sf = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d))
    wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d))
    wf = wr .* ws
    StartUpDG.init_face_data(element_type, quad_rule_face=(rf, sf, wf))
  end
  if D isa AbstractPeriodicDerivativeOperator
    # we do not need any face stuff for periodic operators
    Vf = spzeros(length(wf), length(wq))
  else
    Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask)))
  end
  LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))

  rstf = (rf, sf, tf)
  nrstJ = (nrJ, nsJ, ntJ)

  # low order interpolation nodes
  r1, s1, t1 = StartUpDG.nodes(element_type, 1)
  V1 = StartUpDG.vandermonde(element_type, 1, r, s, t) / StartUpDG.vandermonde(element_type, 1, r1, s1, t1)

  return RefElemData(
    element_type, approximation_type, N,
    face_vertices, V1,
    rst, VDM, face_mask,
    N, rst, LinearAlgebra.I, # plotting
    rstq, wq, Vq, # quadrature
    rstf, wf, Vf, nrstJ, # faces
    M, Pq, Drst, LIFT)
end

function Base.show(io::IO, mime::MIME"text/plain", rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {NDIMS, ElementType<:StartUpDG.AbstractElemShape, ApproximationType<:AbstractDerivativeOperator}
  @nospecialize rd
  print(io, "RefElemData for an approximation using an ")
  show(IOContext(io, :compact => true), rd.approximationType)
  print(io, " on $(rd.elementType) element")
end

function Base.show(io::IO, rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {NDIMS, ElementType<:StartUpDG.AbstractElemShape, ApproximationType<:AbstractDerivativeOperator}
  @nospecialize rd
  print(io, "RefElemData{", summary(rd.approximationType), ", ", rd.elementType, "}")
end

function StartUpDG.inverse_trace_constant(rd::RefElemData{NDIMS, ElementType, ApproximationType})  where {NDIMS, ElementType<:Union{Line, Quad, Hex}, ApproximationType<:AbstractDerivativeOperator}
  D = rd.approximationType

  # the inverse trace constant is the maximum eigenvalue corresponding to
  #       M_f * v = λ * M * v
  # where M_f is the face mass matrix and M is the volume mass matrix.
  # Since M is diagonal and since M_f is just the boundary "mask" matrix
  # (which extracts the first and last entries of a vector), the maximum
  # eigenvalue is the inverse of the first or last mass matrix diagonal.
  left_weight = SummationByPartsOperators.left_boundary_weight(D)
  right_weight = SummationByPartsOperators.right_boundary_weight(D)
  max_eigenvalue = max(inv(left_weight), inv(right_weight))

  # For tensor product elements, the trace constant for higher dimensional
  # elements is the one-dimensional trace constant multiplied by `NDIMS`. See
  #     "GPU-accelerated discontinuous Galerkin methods on hybrid meshes."
  #     Chan, Jesse, et al (2016), https://doi.org/10.1016/j.jcp.2016.04.003
  # for more details (specifically, Appendix A.1, Theorem A.4).
  return NDIMS * max_eigenvalue
end

# type alias for specializing on a periodic SBP operator
const DGMultiPeriodicFDSBP{NDIMS, ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} where {NDIMS, ElemType, ApproxType<:SummationByPartsOperators.AbstractPeriodicDerivativeOperator, SurfaceIntegral, VolumeIntegral}

"""
    DGMultiMesh(dg::DGMulti)

Constructs a single-element [`DGMultiMesh`](@ref) for a single periodic element given
a DGMulti with `approximation_type` set to a periodic (finite difference) SBP operator from
SummationByPartsOperators.jl.
"""
function DGMultiMesh(dg::DGMultiPeriodicFDSBP{NDIMS};
                     coordinates_min=ntuple(_ -> -one(real(dg)), NDIMS),
                     coordinates_max=ntuple(_ -> one(real(dg)), NDIMS)) where {NDIMS}

  rd = dg.basis

  e = ones(size(rd.r))
  z = zero(e)

  VXYZ = ntuple(_ -> [], NDIMS)
  EToV = NaN # StartUpDG.jl uses size(EToV, 1) for the number of elements, this lets us reuse that.
  FToF = []

  # We need to scale the domain from `[-1, 1]^NDIMS` (default in StartUpDG.jl)
  # to the given `coordinates_min, coordinates_max`
  xyz = xyzq = map(copy, rd.rst)
  for dim in 1:NDIMS
    factor = (coordinates_max[dim] - coordinates_min[dim]) / 2
    @. xyz[dim] = factor * (xyz[dim] + 1) + coordinates_min[dim]
  end
  xyzf = ntuple(_ -> [], NDIMS)
  wJq = diag(rd.M)

  # arrays of connectivity indices between face nodes
  mapM = mapP = mapB = []

  # volume geofacs Gij = dx_i/dxhat_j
  coord_diffs = coordinates_max .- coordinates_min

  J_scalar = prod(coord_diffs) / 2^NDIMS
  J = e * J_scalar

  if NDIMS == 1
    rxJ = J_scalar * 2 / coord_diffs[1]
    rstxyzJ = @SMatrix [rxJ * e]
  elseif NDIMS == 2
    rxJ = J_scalar * 2 / coord_diffs[1]
    syJ = J_scalar * 2 / coord_diffs[2]
    rstxyzJ = @SMatrix [rxJ * e z; z syJ * e]
  elseif NDIMS == 3
    rxJ = J_scalar * 2 / coord_diffs[1]
    syJ = J_scalar * 2 / coord_diffs[2]
    tzJ = J_scalar * 2 / coord_diffs[3]
    rstxyzJ = @SMatrix [rxJ * e z z; z syJ * e z; z z tzJ * e]
  end

  # surface geofacs
  nxyzJ = ntuple(_ -> [], NDIMS)
  Jf = []

  periodicity = ntuple(_ -> true, NDIMS)

  md = MeshData(VXYZ, EToV, FToF, xyz, xyzf, xyzq, wJq,
                mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf,
                periodicity)

  boundary_faces = []
  return DGMultiMesh{NDIMS, rd.elementType, typeof(md), typeof(boundary_faces)}(md, boundary_faces)
end

# This is used in `estimate_dt`. `estimate_h` uses that `Jf / J = O(h^{NDIMS-1}) / O(h^{NDIMS}) = O(1/h)`.
# However, since we do not initialize `Jf` for periodic FDSBP operators, we specialize `estimate_h`
# based on the reference grid provided by SummationByPartsOperators.jl and information about the domain size
# provided by `md::MeshData``.
function StartUpDG.estimate_h(e, rd::RefElemData{NDIMS, ElementType, ApproximationType}, md::MeshData)  where {NDIMS, ElementType<:StartUpDG.AbstractElemShape, ApproximationType<:SummationByPartsOperators.AbstractPeriodicDerivativeOperator}
  D = rd.approximationType
  x = grid(D)

  # we assume all SummationByPartsOperators.jl reference grids are rescaled to [-1, 1]
  xmin = SummationByPartsOperators.xmin(D)
  xmax = SummationByPartsOperators.xmax(D)
  factor = 2 / (xmax - xmin)

  # If the domain has size L^NDIMS, then `minimum(md.J)^(1 / NDIMS) = L`.
  # WARNING: this is not a good estimate on anisotropic grids.
  return minimum(diff(x)) * factor * minimum(md.J)^(1 / NDIMS)
end

# specialized for DGMultiPeriodicFDSBP since there are no face nodes
# and thus no inverse trace constant for periodic domains.
function estimate_dt(mesh::DGMultiMesh, dg::DGMultiPeriodicFDSBP)
  rd = dg.basis # RefElemData
  return StartUpDG.estimate_h(rd, mesh.md)
end

# do nothing for interface terms if using a periodic operator
function prolong2interfaces!(cache, u, mesh::DGMultiMesh, equations,
                             surface_integral, dg::DGMultiPeriodicFDSBP)
  @assert nelements(mesh, dg, cache) == 1
  nothing
end

function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
                              mesh::DGMultiMesh,
                              have_nonconservative_terms::Val{false}, equations,
                              dg::DGMultiPeriodicFDSBP)
  @assert nelements(mesh, dg, cache) == 1
  nothing
end

function calc_surface_integral!(du, u, surface_integral::SurfaceIntegralWeakForm,
                                mesh::DGMultiMesh, equations,
                                dg::DGMultiPeriodicFDSBP, cache)
  @assert nelements(mesh, dg, cache) == 1
  nothing
end