text:
"""
    MenuAdventures.ExtraDirections

A sub-module with a bunch of extra directions

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> using MenuAdventures.ExtraDirections

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "ExtraDirections")) do
            check_choices() do interface
                you = Person(
                    "Brandon",
                    grammatical_person = second_person,
                    indefinite_article = "",
                )
                center_room = Room("central room")
                universe = Universe(you, interface = interface)
                universe[center_room, you] = Containing()
                universe[center_room, Room("room")] = West()
                universe[center_room, Room("room")] = NorthWest()
                universe[center_room, Room("room")] = North()
                universe[center_room, Room("room")] = NorthEast()
                universe[center_room, Room("room")] = East()
                universe[center_room, Room("room")] = SouthEast()
                universe[center_room, Room("room")] = South()
                universe[center_room, Room("room")] = SouthWest()
                universe[center_room, Room("room")] = Up()
                universe[center_room, Room("room")] = Down()
                universe[center_room, Room("room")] = Inside()
                universe[center_room, Room("room")] = Outside()
                universe
            end
        end
true
```
"""
module ExtraDirections

import Base: show
using MenuAdventures: Direction
import MenuAdventures: opposite

"""
    NorthWest()

A [`Direction`](@ref).
"""
struct NorthWest <: Direction end

export NorthWest

function show(io::IO, ::NorthWest)
    print(io, "north-west")
end

opposite(::NorthWest) = SouthEast()

"""
    SouthWest()

A [`Direction`](@ref).
"""
struct SouthWest <: Direction end

export SouthWest

function show(io::IO, ::SouthWest)
    print(io, "south-west")
end

opposite(::SouthWest) = NorthEast()

"""
    SouthEast()

A [`Direction`](@ref).
"""
struct SouthEast <: Direction end

export SouthEast

function show(io::IO, ::SouthEast)
    print(io, "south-east")
end

opposite(::SouthEast) = NorthWest()

"""
    NorthEast()

A [`Direction`](@ref).
"""
struct NorthEast <: Direction end

export NorthEast

function show(io::IO, ::NorthEast)
    print(io, "north-east")
end

opposite(::NorthEast) = SouthWest()

"""
    Up()

A [`Direction`](@ref).
"""
struct Up <: Direction end

export Up

function show(io::IO, ::Up)
    print(io, "up")
end

opposite(::Up) = Down()

"""
    Down()

A [`Direction`](@ref).
"""
struct Down <: Direction end

export Down

function show(io::IO, ::Down)
    print(io, "down")
end

opposite(::Down) = Up()

"""
    Inside()

A [`Direction`](@ref).
"""
struct Inside <: Direction end

export Inside

function show(io::IO, ::Inside)
    print(io, "inside")
end

opposite(::Inside) = Outside()

"""
    Outside()

A [`Direction`](@ref).
"""
struct Outside <: Direction end

export Outside

function show(io::IO, ::Outside)
    print(io, "outside")
end

opposite(::Outside) = Inside()

end


text:
using NLsolve

"""
	StaticMatch(mtypes, ftypes, mdist, fdist, surplus)

Construct a Choo & Siow (2006) marriage market model and solve for the equilibrium.
"""
struct StaticMatch

	# m/f types
	"vector of vectors of values for each of N male traits"
	mtypes::Array{Array{T, 1}, 1} where T <: Real
	"vector of vectors of values for each of L female traits"
	ftypes::Array{Array{T, 1}, 1} where T <: Real

	# m/f masses
	"male masses: I_1 x ... x I_N"
	mdist::Array
	"female masses: J_1 x ... x J_L"
	fdist::Array

	"surplus array: (I_1 x ... x I_N) x (J_1 x ... x J_L)"
	surplus::Array

	# equilibrium
	"masses of single males: I_1 x ... x I_N"
	msingle::Array
	"masses of single females: J_1 x ... x J_L"
	fsingle::Array
	"masses of married couples: (I_1 x ... x I_N) x (J_1 x ... x J_L)"
	matches::Array
	"wife's share of marital surplus: (I_1 x ... x I_N) x (J_1 x ... x J_L)"
	wifeshare::Array

	"Inner constructor solves equilibrium and performs sanity checks"
	function StaticMatch(mtypes::Array{Array{T, 1}, 1} where T <: Real,
		                 ftypes::Array{Array{T, 1}, 1} where T <: Real,
		                 mdist::Array, fdist::Array, surplus::Array)

		# CHECK: masses of m/f must be proper probability distro
		if minimum(mdist) < 0 || minimum(fdist) < 0
			error("invalid type distribution")
		end

		# CHECK: masses of m/f match the number of types
		if ndims(mdist) != length(mtypes) || ndims(fdist) != length(ftypes)
			error("type distributions inconsistent with type vectors")
		end

		# compute equilibrium
		msingle, fsingle, matches = equilibrium(surplus, mdist, fdist)
		wifeshare = surplusdiv(matches, msingle, fsingle) ./ (2 .* surplus)

		# TEST: masses of every match outcome must be strictly positive
		negtol = 0
		if minimum(msingle) < -negtol || minimum(fsingle) < -negtol
			@warn "Non-positive mass of singles."
		end
		if minimum(matches) < -negtol
			@warn "Non-positive match mass."
		end

		# create instance
		new(mtypes, ftypes, mdist, fdist, surplus, msingle, fsingle, matches, wifeshare)
	end # constructor

end # struct

"Outer constructor that takes the production function to build the surplus array"
function StaticMatch(men::Array{Array{T, 1}, 1} where T <: Real, wom::Array{Array{T, 1}, 1} where T <: Real,
	                 mmass::Array, fmass::Array, prodfn::Function)
	# Note: prodfn(man::Vector, woman::Vector)
	surp = generate_surplus(men, wom, mmass, fmass, prodfn)
	return StaticMatch(men, wom, mmass, fmass, surp)
end

"Outer constructor for one dimensional case"
function StaticMatch(men::Vector, wom::Vector,
	                 mmass::Array, fmass::Array, prodfn::Function)
	# augment production function
	g(x::Vector, y::Vector) = prodfn(x[1], y[1])
	return StaticMatch([men,], [wom,], mmass, fmass, g)
end

"Outer constructor for one dimensional males case"
function StaticMatch(men::Vector, wom::Array{Array{T, 1}, 1} where T <: Real,
	                 mmass::Array, fmass::Array, prodfn::Function)
	return StaticMatch([men,], wom, mmass, fmass, prodfn)
end

"Outer constructor for one dimensional females case"
function StaticMatch(men::Array{Array{T, 1}, 1} where T <: Real, wom::Vector,
	                 mmass::Array, fmass::Array, prodfn::Function)
	return StaticMatch(men, [wom,], mmass, fmass, prodfn)
end

"""
Quasi-demand functions

Assuming that surplus is generated from equal gains on both sides.
This is w.l.o.g. as gains and transfers cannot be separately identified.
"""
function demand(surp::Array, transfers::Array, mmass::Array, fmass::Array)

	# un-normalized raw demands (Note: surp is defined in per-person terms)
	m_raw_dmd = exp.(surp .- transfers)
	f_raw_dmd = exp.(surp .+ transfers)

	# total raw demands per sex, by summing over opposite sex
	m_raw_dmd_tot = [sum(m_raw_dmd[i,j] for j in CartesianIndices(fmass)) for i in CartesianIndices(mmass)]
	f_raw_dmd_tot = [sum(f_raw_dmd[i,j] for i in CartesianIndices(mmass)) for j in CartesianIndices(fmass)]

	# solve for singles: normalization factors
	m_sng_dmd = mmass ./ (1 .+ m_raw_dmd_tot)
	f_sng_dmd = fmass ./ (1 .+ f_raw_dmd_tot)

	# normalized demands
	m_dmd = [m_sng_dmd[i] * m_raw_dmd[i,j] for i in CartesianIndices(mmass), j in CartesianIndices(fmass)]
	f_dmd = [f_sng_dmd[j] * f_raw_dmd[i,j] for i in CartesianIndices(mmass), j in CartesianIndices(fmass)]

	return m_sng_dmd, m_dmd, f_sng_dmd, f_dmd
end # demand

"""
Compute equilibrium shares of singles and marriages.

Solves the matching equilibrium by searching for the transfers
(analogous to prices) which clear the market.
"""
function equilibrium(surpl::Array, mmass::Array, fmass::Array)

	"Function to pass to solver to find zero."
	function eqnsys!(res::Array, tx::Array)
		msd, mdmd, fsd, fdmd = demand(surpl, tx, mmass, fmass)
		res .= mdmd .- fdmd
		return res
	end

	#initial guess of transfers/prices
	guess = zero(surpl)

	# NLsolve
	result = nlsolve(eqnsys!, guess, ftol = 1e-16, autodiff = :forward, method = :newton)

	if !converged(result)
		@warn "NLsolve failed to converge to equilibrium. Solver result:"
		@show result
	end

	eqm_transfers = reshape(result.zero, size(surpl))
	singmen, mmatches, singwom, fmatches = demand(surpl, eqm_transfers, mmass, fmass)
	matches = (mmatches .+ fmatches) ./ 2 # average out any approximation error

	return singmen, singwom, matches
end # equilibrium

"Construct production array from function."
function generate_surplus(mtypes::Array{Array{T, 1}, 1} where T <: Real,
                          ftypes::Array{Array{T, 1}, 1} where T <: Real,
                          mmass::Array, fmass::Array, prodfn::Function)
	# get dimensions
	Dm = [length(v) for v in mtypes]
	Df = [length(v) for v in ftypes]

	# initialize arrays
	surp = Array{Float64}(undef, Dm..., Df...)
	gent = Vector{Float64}(undef, length(mtypes)) # one man's vector of traits
	lady = Vector{Float64}(undef, length(ftypes))

	for coord in CartesianIndices(surp)
		for trt in 1:length(mtypes) # loop through traits in coord
			gent[trt] = mtypes[trt][coord[trt]]
		end
		for trt in 1:length(ftypes) # loop through traits in coord
			lady[trt] = ftypes[trt][coord[trt+length(mtypes)]]
		end
		surp[coord] = prodfn(gent, lady)
	end

	return surp
end # generate_surplus

"Wife's consumption out of surplus (aggregate share)."
function surplusdiv(matches::Array, sm::Array, sw::Array)
	return [log(matches[i,j]) - log(sw[j]) for i in CartesianIndices(sm), j in CartesianIndices(sw)]
end # surplusdiv

"Estimate the marital surplus from observed matches and singles."
function estimate_static_surplus(matches::Array, msingle::Array, fsingle::Array)
	return [log(matches[i,j]) - 0.5 * (log(msingle[i]) + log(fsingle[j]))
	        for i in CartesianIndices(msingle), j in CartesianIndices(fsingle)]
end # estimate_static_surplus

ext: jl


text:
abstract type AbstractDeviceFormulation end

function _validate_device_formulation(::Type{D}) where D<:Union{AbstractDeviceFormulation,
                                                                PSY.Device}

    if !isconcretetype(D)
        throw(ArgumentError("The device model must contain only concrete types, $(D) is an Abstract Type"))
    end

end

mutable struct DeviceModel{D<:PSY.Device,
                           B<:AbstractDeviceFormulation}
    device_type::Type{D}
    formulation::Type{B}
    feedforward::Union{Nothing, AbstractFeedForwardAffect}

    function DeviceModel(::Type{D},
                    ::Type{B},
                    feedforward::Union{Nothing, AbstractFeedForwardAffect}) where {D<:PSY.Device,
                                                            B<:AbstractDeviceFormulation}

    _validate_device_formulation(D)
    _validate_device_formulation(B)

    new{D, B}(D, B, feedforward)

    end

end

"""
    DeviceModel(::Type{D}, ::Type{B}) where {D<:PSY.Device,
                                       B<:AbstractDeviceFormulation}
This validates the device formulation for the Power System Device and the
abstract device formulation and returns  Power System Device and the
abstract device formulation if the power system device is a concrete type.

# Arguments
-`::Type{D}`: Power System Device
-`::Type{B}`: Abstract Device Formulation

# Outputs
`DeviceModel(D, B, nothing)`: D::PSY.Device, B::AbstractDeviceFormulation

# Example
```julia
branches = Dict{Symbol, DeviceModel}
    (:L => DeviceModel(PSY.Line, PSI.StaticLine),
    :T => DeviceModel(PSY.Transformer2W, PSI.StaticTransformer),
    :TT => DeviceModel(PSY.TapTransformer , PSI.StaticTransformer),
    :dc_line => DeviceModel(PSY.HVDCLine, PSI.HVDCDispatch))
```
"""
function DeviceModel(::Type{D},
                     ::Type{B}) where {D<:PSY.Device,
                                       B<:AbstractDeviceFormulation}

                    _validate_device_formulation(D)
                    _validate_device_formulation(B)

    return DeviceModel(D, B, nothing)

end


text:
import Lumberjack
dbglog(x::Any) = Lumberjack.debug(string(x))
dbglog(x::Any...) = Lumberjack.debug(string(x))
infolog(x::Any) = Lumberjack.info(string(x))
infolog(x::Any...) = Lumberjack.info(string(x))

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# function plots trajectory
###############

function plotTraj(xp,up,N,ego,L,nOb,vOb,lOb,disp_title,plotNumb)

	# obcenter1 = [(ob1[1]+ob1[3])/2-ob1[3];
	#             	 (ob1[2]+ob1[4])/2-ob1[4]]
	#
	# obcenter2 = [(ob2[1]+ob2[3])/2-ob2[3];
	#              (ob2[2]+ob2[4])/2-ob2[4]]
	#
	# obcenter3 = [(ob3[1]+ob3[3])/2-ob3[3];
	#              (ob3[2]+ob3[4])/2-ob3[4]]

    W_ev = ego[2]+ego[4]
	L_ev = ego[1]+ego[3]

	# W_tv1 = ob1[2]+ob1[4]
	# L_tv1 = ob1[1]+ob1[3]
	#
	# W_tv2 = ob2[2]+ob2[4]
	# L_tv2 = ob2[1]+ob2[3]
	#
	# W_tv3 = ob3[2]+ob3[4]
	# L_tv3 = ob3[1]+ob3[3]

	up = [up ; zeros(1,2)]	# final position no input

	w = W_ev/2;
	offset = L_ev/2 - ego[3]

	# initial state
	x0_s = xp[1,:]
	Rot0 = [cos(x0_s[3]) -sin(x0_s[3]); sin(x0_s[3]) cos(x0_s[3])]
	x0 = [x0_s[1]; x0_s[2]]
    centerCar0 = x0 + Rot0*[offset;0]

	# end state
	xF_s = xp[end,:]
	RotF = [cos(xF_s[3]) -sin(xF_s[3]); sin(xF_s[3]) cos(xF_s[3])]
	xF = [xF_s[1]; xF_s[2]]
    centerCarF = xF + RotF*[offset;0]

	for i = 1:1:N+1

		figure(plotNumb)
		plot(xp[1:i,1],xp[1:i,2],"b") 	# plot trajectory so far
		title(disp_title)
		hold(1)

		# plot trajectory
		for j = 1 : nOb
			for k = 1 : vOb[j]
				plot([lOb[j][k][1],lOb[j][k+1][1]] , [lOb[j][k][2],lOb[j][k+1][2]] ,"k")
			end
		end

	    Rot = [cos(xp[i,3]) -sin(xp[i,3]);sin(xp[i,3]) cos(xp[i,3])]

		x_cur = [xp[i,1];
	          xp[i,2]]

	    centerCar = x_cur + Rot*[offset;0]

	    carBox(centerCar,xp[i,3],W_ev/2,L_ev/2)
	    carBox(x_cur + (Rot*[L;w-0.15]), xp[i,3] + up[i,1],0.15,0.3)
	    carBox(x_cur + (Rot*[L;-w+0.15]),xp[i,3] + up[i,1],0.15,0.3)
	    carBox(x_cur + (Rot*[0; w-0.15]) ,xp[i,3],0.15,0.3)
	    carBox(x_cur + (Rot*[0;-w+0.15]) ,xp[i,3],0.15,0.3)

		# plot start position
		plot(x0[1],x0[2],"ob")
		carBox(centerCar0,x0_s[3],W_ev/2,L_ev/2)
		carBox(x0 + (Rot0*[L;w-0.15])  ,x0_s[3],0.15,0.3)
		carBox(x0 + (Rot0*[L;-w+0.15]) ,x0_s[3],0.15,0.3)
	    carBox(x0 + (Rot0*[0; w-0.15]) ,x0_s[3], 0.15,0.3)
		carBox(x0 + (Rot0*[0;-w+0.15]) ,x0_s[3], 0.15,0.3)

		# plot end position
		carBox_dashed(centerCarF,xF_s[3],W_ev/2,L_ev/2)
		carBox_dashed(xF + (RotF*[L;w-0.15])  ,xF_s[3],0.15,0.3)
		carBox_dashed(xF + (RotF*[L;-w+0.15]) ,xF_s[3],0.15,0.3)
	    carBox_dashed(xF + (RotF*[0; w-0.15]) ,xF_s[3], 0.15,0.3)
		carBox_dashed(xF + (RotF*[0;-w+0.15]) ,xF_s[3], 0.15,0.3)
		if i == N+1
			plot(xF[1],xF[2],"ob")
		end

	    axis("equal")

	    hold(0)

	    sleep(0.05)
	end
end

# plot cars
function carBox(x0,phi,w,l)
    car1 = x0[1:2] + [cos(phi)*l;sin(phi)*l] + [sin(phi)*w;-cos(phi)*w];
    car2 = x0[1:2] + [cos(phi)*l;sin(phi)*l] - [sin(phi)*w;-cos(phi)*w];
    car3 = x0[1:2] - [cos(phi)*l;sin(phi)*l] + [sin(phi)*w;-cos(phi)*w];
    car4 = x0[1:2] - [cos(phi)*l;sin(phi)*l] - [sin(phi)*w;-cos(phi)*w];
    plot([car1[1],car2[1],car4[1],car3[1],car1[1]],[car1[2],car2[2],car4[2],car3[2],car1[2]],"k")
end

# plot cars
function carBox_dashed(x0,phi,w,l)
    car1 = x0[1:2] + [cos(phi)*l;sin(phi)*l] + [sin(phi)*w;-cos(phi)*w];
    car2 = x0[1:2] + [cos(phi)*l;sin(phi)*l] - [sin(phi)*w;-cos(phi)*w];
    car3 = x0[1:2] - [cos(phi)*l;sin(phi)*l] + [sin(phi)*w;-cos(phi)*w];
    car4 = x0[1:2] - [cos(phi)*l;sin(phi)*l] - [sin(phi)*w;-cos(phi)*w];
    plot([car1[1],car2[1],car4[1],car3[1],car1[1]],[car1[2],car2[2],car4[2],car3[2],car1[2]],":k")
end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
@testset "wrong theta_init" begin
    method = :CICO_ONE_PASS
    @test_throws ArgumentError get_endpoint(
        [-10., 2., 2.1],
        1,
        f_3p_1im_dep,
        method;
        loss_crit = 9.
        )
end

@testset "wrong scan_bound" begin
    method = :CICO_ONE_PASS
    @test_throws ArgumentError get_endpoint(
        [3., 2., 2.1],
        1,
        f_3p_1im_dep,
        method;
        loss_crit = 9.,
        scan_bound = 6.,
        theta_bounds = [(-Inf,5.),(-Inf,Inf),(-Inf,Inf)]
        )
end

@testset "wrong scan_bound and theta_init" begin
    method = :CICO_ONE_PASS
    @test_throws ArgumentError get_endpoint(
        [3., 2., 2.1],
        1,
        f_3p_1im_dep,
        method;
        loss_crit = 9.,
        scan_bound = 2.
        )
end

@testset "wrong theta_init in log scale" begin
    method = :CICO_ONE_PASS
    @test_throws ArgumentError get_endpoint(
        [3., 2., -1],
        1,
        f_3p_1im_dep,
        method;
        loss_crit = 9.,
        scale = [:direct,:direct,:log]
        )
end

@testset "wrong theta_bounds in log scale" begin
    method = :CICO_ONE_PASS
    @test_throws ArgumentError get_endpoint(
        [3., 2., 2.],
        3,
        f_3p_1im_dep,
        method;
        loss_crit = 9.,
        scale = [:direct,:direct,:log],
        theta_bounds = [(-Inf, Inf), (-Inf, Inf), (-5., Inf)]
        )
end


text:
module WallBooster

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Wall Booster (Right)" => Ahorn.EntityPlacement(
        Maple.WallBooster,
        "rectangle",
        Dict{String, Any}(
            "left" => true
        )
    ),
    "Wall Booster (Left)" => Ahorn.EntityPlacement(
        Maple.WallBooster,
        "rectangle",
        Dict{String, Any}(
            "left" => false
        )
    )
)

Ahorn.minimumSize(entity::Maple.WallBooster) = 0, 8
Ahorn.resizable(entity::Maple.WallBooster) = false, true

function Ahorn.selection(entity::Maple.WallBooster)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.WallBooster, room::Maple.Room)
    left = get(entity.data, "left", false)

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    tileHeight = div(height, 8)

    if left
        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "objects/wallBooster/fireMid00", 0, (i - 1) * 8)
        end

        Ahorn.drawImage(ctx, "objects/wallBooster/fireTop00", 0, 0)
        Ahorn.drawImage(ctx, "objects/wallBooster/fireBottom00", 0, (tileHeight - 1) * 8)

    else
        Ahorn.Cairo.save(ctx)
        Ahorn.scale(ctx, -1, 1)

        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "objects/wallBooster/fireMid00", -8, (i - 1) * 8)
        end

        Ahorn.drawImage(ctx, "objects/wallBooster/fireTop00", -8, 0)
        Ahorn.drawImage(ctx, "objects/wallBooster/fireBottom00", -8, (tileHeight - 1) * 8)

        Ahorn.restore(ctx)
    end
end

# Offset X position so it flips in place
function Ahorn.flipped(entity::Maple.WallBooster, horizontal::Bool)
    if horizontal
        entity.left = !entity.left
        entity.x += entity.left ? 8 : -8

        return entity
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
#this structure holds the partitioned portfolio
struct TrisectedPortfolio
  P::Portfolio
  partitions::Vector{Portfolio}
  assignments::Vector{Int}

  wₛ::Vector{Float64} #sums of the weights of the partitions
  S::Matrix{Float64} #covariance matrix of the portfolio
end

#constructor
function TrisectedPortfolio(P::Portfolio, start1=1)::TrisectedPortfolio
  targetPartitionSize::Int = P.U.M ÷ 3
  assignments::Vector{Int} = fill(3, P.U.M)

  #do the initial assignments
  assignments[1:targetPartitionSize] .= 1
  assignments[(targetPartitionSize+1):(2*targetPartitionSize)] .= 2

  #create 3 new portfolios, one for each partition
  #the portfolios are a placeholder until the refresh command is called as part of rotatetrisected
  partitions::Vector{Portfolio} = (i::Int->Portfolio(P, P.wₜ)).(1:3)

  #make the trisected portfolio
  TP::TrisectedPortfolio = TrisectedPortfolio(
    P,
    partitions,
    assignments,
    Vector{Float64}(undef,3),
    Matrix{Float64}(undef,3,3)
  )

  #rotate the portfolio, which will also compute wₛ and S
  rotatetrisected!(TP, start1 - 1)

  return TP
end

#computes the covariance of the partitions and the test portfolio
function Statistics.cov(TG::TrisectedPortfolio, P::Portfolio)::Vector{Float64}
  local SGP::Vector{Float64} = Vector{Float64}(undef, 3)

  for i::Int ∈ 1:3
    SGP[i] = cov(TG.partitions[i], P)
  end

  return SGP
end

#rotate the assignments by a set number
function rotatetrisected!(TP::TrisectedPortfolio, by::Int = 1)::TrisectedPortfolio

  #rotate the assignments
  for i::Int ∈ 1:TP.P.U.M
    TP.assignments[i] = TP.assignments[1 + ((i + by - 1) % (TP.P.U.M))]
  end

  #refresh the remaining contents given the changes
  refreshtrisected!(TP)
  return TP
end

function refreshtrisected!(TP::TrisectedPortfolio)::TrisectedPortfolio
  #get the sum of the weights for each parittion
  TP.wₛ .= (i::Int->sum(TP.P.wₜ .* (i .== TP.assignments))).(1:3)

  #update the portfolio partitions
  (i::Int->Portfolio!(TP.partitions[i], TP.P.wₜ .* (i .== TP.assignments), normalizeto=TP.wₛ[i])).(1:3)

  #update the covariance matrix
  cov!(TP.partitions, TP.S)

  return TP
end

#refreshes the trisected portfolio given a new set of target weights
function refreshtrisected!(TP, wₜ)::TrisectedPortfolio
  Portfolio!(TP.P, wₜ)

  return refreshtrisected!(TP)
end

function getω₁R(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64})
  @inbounds ω₁::Float64 = (
    (TG.S[2,2] * (SGP[3] * TG.wₛ[1] + (-1) * SGP[1] * TG.wₛ[3])^2 + SGP[2]^2 * (TG.S[3,3]
    * TG.wₛ[1]^2 + TG.wₛ[3] * ((-2) * TG.S[1,3] * TG.wₛ[1] + TG.S[1,1] * TG.wₛ[3])) +
    (-2) * SGP[2] * (SGP[1] * TG.S[3,3] * TG.wₛ[1] * TG.wₛ[2] + (-1) * TG.S[1,3] * SGP[3]
    * TG.wₛ[1] * TG.wₛ[2] + (-1) * TG.S[1,2] * SGP[3] * TG.wₛ[1] * TG.wₛ[3] + (-1) *
    TG.S[1,3] * SGP[1] * TG.wₛ[2] * TG.wₛ[3] + TG.S[1,1] * SGP[3] * TG.wₛ[2] * TG.wₛ[3]
    + TG.S[1,2] * SGP[1] * TG.wₛ[3]^2 + TG.S[2,3] * TG.wₛ[1] * (SGP[3] * TG.wₛ[1] + (-1)
    * SGP[1] * TG.wₛ[3]) ) + TG.wₛ[2] * (SGP[3]^2 * ((-2) * TG.S[1,2] * TG.wₛ[1] + TG.S[1,1]
    * TG.wₛ[2]) + 2 * SGP[1] * SGP[3] * (TG.S[2,3] * TG.wₛ[1] + ( -1) * TG.S[1,3] * TG.wₛ[2]
    + TG.S[1,2] * TG.wₛ[3]) + SGP[1]^2 * (TG.S[3,3] * TG.wₛ[2] + (-2) * TG.S[2,3] * TG.wₛ[3])))^(-1)
    * (SGP[2]^2 * TG.S[3,3] * TG.wₛ[1] + (-2) * TG.S[2,3] * SGP[2] * SGP[3] * TG.wₛ[1]
    + TG.S[2,2] * SGP[3]^2 * TG.wₛ[1] + (-1) * SGP[1] * SGP[2] * TG.S[3,3] * TG.wₛ[2]
    + SGP[1] * TG.S[2,3] * SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.S[3,3] * Θ.SGP * TG.wₛ[1]
    * TG.wₛ[2] + TG.S[2,3] * SGP[3] * Θ.SGP * TG.wₛ[1] * TG.wₛ[2] + SGP[1] * TG.S[3,3]
    * Θ.SGP * TG.wₛ[2]^2 + SGP[1] * TG.S[2,3] * SGP[2] * TG.wₛ[3] + ( -1) * SGP[1] *
    TG.S[2,2] * SGP[3] * TG.wₛ[3] + TG.S[2,3] * SGP[2] * Θ.SGP * TG.wₛ[1] * TG.wₛ[3]
    + (-1) * TG.S[2,2] * SGP[3] * Θ.SGP * TG.wₛ[1] * TG.wₛ[3] + (-2) * SGP[1] * TG.S[2,3]
    * Θ.SGP * TG.wₛ[2] * TG.wₛ[3] + SGP[1] * TG.S[2,2] * Θ.SGP * TG.wₛ[3]^2 + TG.S[1,3]
    * SGP[2] * (SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.wₛ[3]) + (-1) * TG.S[1,2] * SGP[3]
    * (SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.wₛ[3]) + (-1) * TG.S[1,3] * Θ.SGP * TG.wₛ[2]
    * (SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.wₛ[3]) + TG.S[1,2] * Θ.SGP * TG.wₛ[3] *
    (SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.wₛ[3]) + (-1) * (SGP[3] * TG.wₛ[2] + (-1)
    * SGP[2] * TG.wₛ[3])^2 * ((SGP[3] * TG.wₛ[2] + (-1) * SGP[2] * TG.wₛ[3])^(-2) * ((-1)
    * TG.S[1,1] * SGP[2]^2 * TG.S[3,3] + 2 * TG.S[1,1] * TG.S[2,3] * SGP[2] * SGP[3]
    + TG.S[1,2]^2 * SGP[3]^2 + (-1) * TG.S[1,1] * TG.S[2,2] * SGP[3]^2 + (-2) * TG.S[1,2]
    * SGP[2] * TG.S[3,3] * Θ.SGP * TG.wₛ[1] + 2 * TG.S[1,2] * TG.S[2,3] * SGP[3] * Θ.SGP
    * TG.wₛ[1] + SGP[2]^2 * TG.S[3,3] * Θ.SG * TG.wₛ[1]^2 + (-2) * TG.S[2,3] * SGP[2]
    * SGP[3] * Θ.SG * TG.wₛ[1]^2 + TG.S[2,2] * SGP[3]^2 * Θ.SG * TG.wₛ[1]^2 + TG.S[2,3]^2
    * Θ.SGP^2 * TG.wₛ[1]^2 + (-1) * TG.S[2,2] * TG.S[3,3] * Θ.SGP^2 * TG.wₛ[1]^2 + 2
    * TG.S[1,1] * SGP[2] * TG.S[3,3] * Θ.SGP * TG.wₛ[2] + (-2) * TG.S[1,1] * TG.S[2,3]
    * SGP[3] * Θ.SGP * TG.wₛ[2] + (-2) * TG.S[1,2] * SGP[3]^2 * Θ.SG * TG.wₛ[1] * TG.wₛ[2]
    + 2 * TG.S[1,2] * TG.S[3,3] * Θ.SGP^2 * TG.wₛ[1] * TG.wₛ[2] + TG.S[1,1] * SGP[3]^2
    * Θ.SG * TG.wₛ[2]^2 + (-1) * TG.S[1,1] * TG.S[3,3] * Θ.SGP^2 * TG.wₛ[2]^2 + TG.S[1,3]^2
    * (SGP[2] + (-1) * Θ.SGP * TG.wₛ[2])^2 + (-2) * TG.S[1,1] * TG.S[2,3] * SGP[2] *
    Θ.SGP * TG.wₛ[3] + (-2) * TG.S[1,2]^2 * SGP[3] * Θ.SGP * TG.wₛ[3] + 2 * TG.S[1,1]
    * TG.S[2,2] * SGP[3] * Θ.SGP * TG.wₛ[3] + 2 * TG.S[1,2] * SGP[2] * SGP[3] * Θ.SG
    * TG.wₛ[1] * TG.wₛ[3] + (-2) * TG.S[1,2] * TG.S[2,3] * Θ.SGP^2 * TG.wₛ[1] * TG.wₛ[3]
    + (-2) * TG.S[1,1] * SGP[2] * SGP[3] * Θ.SG * TG.wₛ[2] * TG.wₛ[3] + 2 * TG.S[1,1]
    * TG.S[2,3] * Θ.SGP^2 * TG.wₛ[2] * TG.wₛ[3] + TG.S[1,1] * SGP[2]^2 * Θ.SG * TG.wₛ[3]^2
    + TG.S[1,2]^2 * Θ.SGP^2 * TG.wₛ[3]^2 + ( -1) * TG.S[1,1] * TG.S[2,2] * Θ.SGP^2 *
    TG.wₛ[3]^2 + SGP[1]^2 * (TG.S[2,3]^2 + (-1) * TG.S[2,2] * TG.S[3,3] + TG.S[3,3] *
    Θ.SG * TG.wₛ[2]^2 + (-2) * TG.S[2,3] * Θ.SG * TG.wₛ[2] * TG.wₛ[3] + TG.S[2,2] * Θ.SG
    * TG.wₛ[3]^2) + (-2) * SGP[1] * (TG.S[1,3] * (( -1) * TG.S[2,2] * SGP[3] + SGP[3]
    * Θ.SG * TG.wₛ[2]^2 + TG.S[2,3] * (SGP[2] + (-1) * Θ.SGP * TG.wₛ[2]) + TG.S[2,2]
    * Θ.SGP * TG.wₛ[3] + (-1) * SGP[2] * Θ.SG * TG.wₛ[2] * TG.wₛ[3]) + TG.S[1,2] * ((-1)
    * SGP[2] * TG.S[3,3] + TG.S[2,3] * SGP[3] + TG.S[3,3] * Θ.SGP * TG.wₛ[2] + (-1) *
    TG.S[2,3] * Θ.SGP * TG.wₛ[3] + (-1) * SGP[3] * Θ.SG * TG.wₛ[2] * TG.wₛ[3] + SGP[2]
    * Θ.SG * TG.wₛ[3]^2) + TG.wₛ[1] * (TG.S[2,3]^2 * Θ.SGP + (-1) * TG.S[2,2] * TG.S[3,3]
    * Θ.SGP + SGP[2] * TG.S[3,3] * Θ.SG * TG.wₛ[2] + TG.S[2,2] * SGP[3] * Θ.SG * TG.wₛ[3]
    + (-1) * TG.S[2,3] * Θ.SG * (SGP[3] * TG.wₛ[2] + SGP[2] * TG.wₛ[3]))) + (-2) * TG.S[1,3]
    * ( TG.S[1,2] * (SGP[2] + (-1) * Θ.SGP * TG.wₛ[2]) * (SGP[3] + (-1) * Θ.SGP * TG.wₛ[3])
    + TG.wₛ[1] * (TG.S[2,3] * Θ.SGP * (( -1) * SGP[2] + Θ.SGP * TG.wₛ[2]) + SGP[2] *
    Θ.SG * ((-1) * SGP[3] * TG.wₛ[2] + SGP[2] * TG.wₛ[3]) + TG.S[2,2] * Θ.SGP * ( SGP[3]
    + (-1) * Θ.SGP * TG.wₛ[3])))))^(1/2))
  )

  return ω₁
end

function getω₁C1(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64})
  CSGP::Vector{ComplexF64} = (ComplexF64).(SGP) #this is necessary due to potential roundoff error
  CS::Matrix{ComplexF64} = (ComplexF64).(TG.S)
  Cwₛ::Vector{ComplexF64} = (ComplexF64).(TG.wₛ)

  @inbounds ω₁::Float64 = real(
    (CS[2,2] * (CSGP[3] * Cwₛ[1] + (-1) * CSGP[1] * Cwₛ[3])^2 + CSGP[2]^2 * (CS[3,3]
    * Cwₛ[1]^2 + Cwₛ[3] * ((-2) * CS[1,3] * Cwₛ[1] + CS[1,1] * Cwₛ[3])) + (-2) * CSGP[2]
    * (CSGP[1] * CS[3,3] * Cwₛ[1] * Cwₛ[2] + (-1) * CS[1,3] * CSGP[3] * Cwₛ[1] * Cwₛ[2]
    + (-1) * CS[1,2] * CSGP[3] * Cwₛ[1] * Cwₛ[3] + (-1) * CS[1,3] * CSGP[1] * Cwₛ[2]
    * Cwₛ[3] + CS[1,1] * CSGP[3] * Cwₛ[2] * Cwₛ[3] + CS[1,2] * CSGP[1] * Cwₛ[3]^2 + CS[2,3]
    * Cwₛ[1] * (CSGP[3] * Cwₛ[1] + (-1) * CSGP[1] * Cwₛ[3]) ) + Cwₛ[2] * (CSGP[3]^2 *
    ((-2) * CS[1,2] * Cwₛ[1] + CS[1,1] * Cwₛ[2]) + 2 * CSGP[1] * CSGP[3] * (CS[2,3] *
    Cwₛ[1] + ( -1) * CS[1,3] * Cwₛ[2] + CS[1,2] * Cwₛ[3]) + CSGP[1]^2 * (CS[3,3] * Cwₛ[2]
    + (-2) * CS[2,3] * Cwₛ[3])))^(-1) * (CSGP[2]^2 * CS[3,3] * Cwₛ[1] + (-2) * CS[2,3]
    * CSGP[2] * CSGP[3] * Cwₛ[1] + CS[2,2] * CSGP[3]^2 * Cwₛ[1] + (-1) * CSGP[1] * CSGP[2]
    * CS[3,3] * Cwₛ[2] + CSGP[1] * CS[2,3] * CSGP[3] * Cwₛ[2] + (-1) * CSGP[2] * CS[3,3]
    * Θ.SGP * Cwₛ[1] * Cwₛ[2] + CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[1] * Cwₛ[2] + CSGP[1]
    * CS[3,3] * Θ.SGP * Cwₛ[2]^2 + CSGP[1] * CS[2,3] * CSGP[2] * Cwₛ[3] + ( -1) * CSGP[1]
    * CS[2,2] * CSGP[3] * Cwₛ[3] + CS[2,3] * CSGP[2] * Θ.SGP * Cwₛ[1] * Cwₛ[3] + (-1)
    * CS[2,2] * CSGP[3] * Θ.SGP * Cwₛ[1] * Cwₛ[3] + (-2) * CSGP[1] * CS[2,3] * Θ.SGP
    * Cwₛ[2] * Cwₛ[3] + CSGP[1] * CS[2,2] * Θ.SGP * Cwₛ[3]^2 + CS[1,3] * CSGP[2] * (CSGP[3]
    * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3]) + (-1) * CS[1,2] * CSGP[3] * (CSGP[3] * Cwₛ[2]
    + (-1) * CSGP[2] * Cwₛ[3]) + (-1) * CS[1,3] * Θ.SGP * Cwₛ[2] * (CSGP[3] * Cwₛ[2]
    + (-1) * CSGP[2] * Cwₛ[3]) + CS[1,2] * Θ.SGP * Cwₛ[3] * (CSGP[3] * Cwₛ[2] + (-1)
    * CSGP[2] * Cwₛ[3]) + (-1) * (CSGP[3] * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3])^2 * ((CSGP[3]
    * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3])^(-2) * ((-1) * CS[1,1] * CSGP[2]^2 * CS[3,3]
    + 2 * CS[1,1] * CS[2,3] * CSGP[2] * CSGP[3] + CS[1,2]^2 * CSGP[3]^2 + (-1) * CS[1,1]
    * CS[2,2] * CSGP[3]^2 + (-2) * CS[1,2] * CSGP[2] * CS[3,3] * Θ.SGP * Cwₛ[1] + 2 *
    CS[1,2] * CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[1] + CSGP[2]^2 * CS[3,3] * Θ.SG * Cwₛ[1]^2
    + (-2) * CS[2,3] * CSGP[2] * CSGP[3] * Θ.SG * Cwₛ[1]^2 + CS[2,2] * CSGP[3]^2 * Θ.SG
    * Cwₛ[1]^2 + CS[2,3]^2 * Θ.SGP^2 * Cwₛ[1]^2 + (-1) * CS[2,2] * CS[3,3] * Θ.SGP^2
    * Cwₛ[1]^2 + 2 * CS[1,1] * CSGP[2] * CS[3,3] * Θ.SGP * Cwₛ[2] + (-2) * CS[1,1] *
    CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[2] + (-2) * CS[1,2] * CSGP[3]^2 * Θ.SG * Cwₛ[1] *
    Cwₛ[2] + 2 * CS[1,2] * CS[3,3] * Θ.SGP^2 * Cwₛ[1] * Cwₛ[2] + CS[1,1] * CSGP[3]^2
    * Θ.SG * Cwₛ[2]^2 + (-1) * CS[1,1] * CS[3,3] * Θ.SGP^2 * Cwₛ[2]^2 + CS[1,3]^2 * (CSGP[2]
    + (-1) * Θ.SGP * Cwₛ[2])^2 + (-2) * CS[1,1] * CS[2,3] * CSGP[2] * Θ.SGP * Cwₛ[3]
    + (-2) * CS[1,2]^2 * CSGP[3] * Θ.SGP * Cwₛ[3] + 2 * CS[1,1] * CS[2,2] * CSGP[3] *
    Θ.SGP * Cwₛ[3] + 2 * CS[1,2] * CSGP[2] * CSGP[3] * Θ.SG * Cwₛ[1] * Cwₛ[3] + (-2)
    * CS[1,2] * CS[2,3] * Θ.SGP^2 * Cwₛ[1] * Cwₛ[3] + (-2) * CS[1,1] * CSGP[2] * CSGP[3]
    * Θ.SG * Cwₛ[2] * Cwₛ[3] + 2 * CS[1,1] * CS[2,3] * Θ.SGP^2 * Cwₛ[2] * Cwₛ[3] + CS[1,1]
    * CSGP[2]^2 * Θ.SG * Cwₛ[3]^2 + CS[1,2]^2 * Θ.SGP^2 * Cwₛ[3]^2 + ( -1) * CS[1,1]
    * CS[2,2] * Θ.SGP^2 * Cwₛ[3]^2 + CSGP[1]^2 * (CS[2,3]^2 + (-1) * CS[2,2] * CS[3,3]
    + CS[3,3] * Θ.SG * Cwₛ[2]^2 + (-2) * CS[2,3] * Θ.SG * Cwₛ[2] * Cwₛ[3] + CS[2,2] *
    Θ.SG * Cwₛ[3]^2) + (-2) * CSGP[1] * (CS[1,3] * (( -1) * CS[2,2] * CSGP[3] + CSGP[3]
    * Θ.SG * Cwₛ[2]^2 + CS[2,3] * (CSGP[2] + (-1) * Θ.SGP * Cwₛ[2]) + CS[2,2] * Θ.SGP
    * Cwₛ[3] + (-1) * CSGP[2] * Θ.SG * Cwₛ[2] * Cwₛ[3]) + CS[1,2] * ((-1) * CSGP[2] *
    CS[3,3] + CS[2,3] * CSGP[3] + CS[3,3] * Θ.SGP * Cwₛ[2] + (-1) * CS[2,3] * Θ.SGP *
    Cwₛ[3] + (-1) * CSGP[3] * Θ.SG * Cwₛ[2] * Cwₛ[3] + CSGP[2] * Θ.SG * Cwₛ[3]^2) + Cwₛ[1]
    * (CS[2,3]^2 * Θ.SGP + (-1) * CS[2,2] * CS[3,3] * Θ.SGP + CSGP[2] * CS[3,3] * Θ.SG
    * Cwₛ[2] + CS[2,2] * CSGP[3] * Θ.SG * Cwₛ[3] + (-1) * CS[2,3] * Θ.SG * (CSGP[3] *
    Cwₛ[2] + CSGP[2] * Cwₛ[3]))) + (-2) * CS[1,3] * ( CS[1,2] * (CSGP[2] + (-1) * Θ.SGP
    * Cwₛ[2]) * (CSGP[3] + (-1) * Θ.SGP * Cwₛ[3]) + Cwₛ[1] * (CS[2,3] * Θ.SGP * (( -1)
    * CSGP[2] + Θ.SGP * Cwₛ[2]) + CSGP[2] * Θ.SG * ((-1) * CSGP[3] * Cwₛ[2] + CSGP[2]
    * Cwₛ[3]) + CS[2,2] * Θ.SGP * ( CSGP[3] + (-1) * Θ.SGP * Cwₛ[3])))))^(1/2))
  )

  return ω₁
end

function getω₁C2(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64})
  CSGP::Vector{ComplexF64} = (ComplexF64).(SGP) #this is necessary due to potential roundoff error
  CS::Matrix{ComplexF64} = (ComplexF64).(TG.S)
  Cwₛ::Vector{ComplexF64} = (ComplexF64).(TG.wₛ)

  @inbounds ω₁::Float64 = real(
    (CS[2,2] * (CSGP[3] * Cwₛ[1] + (-1) * CSGP[1] * Cwₛ[3])^2 + CSGP[2]^2 * (CS[3,3]
    * Cwₛ[1]^2 + Cwₛ[3] * ((-2) * CS[1,3] * Cwₛ[1] + CS[1,1] * Cwₛ[3])) + (-2) * CSGP[2]
    * (CSGP[1] * CS[3,3] * Cwₛ[1] * Cwₛ[2] + (-1) * CS[1,3] * CSGP[3] * Cwₛ[1] * Cwₛ[2]
    + (-1) * CS[1,2] * CSGP[3] * Cwₛ[1] * Cwₛ[3] + (-1) * CS[1,3] * CSGP[1] * Cwₛ[2]
    * Cwₛ[3] + CS[1,1] * CSGP[3] * Cwₛ[2] * Cwₛ[3] + CS[1,2] * CSGP[1] * Cwₛ[3]^2 + CS[2,3]
    * Cwₛ[1] * (CSGP[3] * Cwₛ[1] + (-1) * CSGP[1] * Cwₛ[3]) ) + Cwₛ[2] * (CSGP[3]^2 *
    ((-2) * CS[1,2] * Cwₛ[1] + CS[1,1] * Cwₛ[2]) + 2 * CSGP[1] * CSGP[3] * (CS[2,3] *
    Cwₛ[1] + ( -1) * CS[1,3] * Cwₛ[2] + CS[1,2] * Cwₛ[3]) + CSGP[1]^2 * (CS[3,3] * Cwₛ[2]
    + (-2) * CS[2,3] * Cwₛ[3])))^(-1) * (CSGP[2]^2 * CS[3,3] * Cwₛ[1] + (-2) * CS[2,3]
    * CSGP[2] * CSGP[3] * Cwₛ[1] + CS[2,2] * CSGP[3]^2 * Cwₛ[1] + (-1) * CSGP[1] * CSGP[2]
    * CS[3,3] * Cwₛ[2] + CSGP[1] * CS[2,3] * CSGP[3] * Cwₛ[2] + (-1) * CSGP[2] * CS[3,3]
    * Θ.SGP * Cwₛ[1] * Cwₛ[2] + CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[1] * Cwₛ[2] + CSGP[1]
    * CS[3,3] * Θ.SGP * Cwₛ[2]^2 + CSGP[1] * CS[2,3] * CSGP[2] * Cwₛ[3] + ( -1) * CSGP[1]
    * CS[2,2] * CSGP[3] * Cwₛ[3] + CS[2,3] * CSGP[2] * Θ.SGP * Cwₛ[1] * Cwₛ[3] + (-1)
    * CS[2,2] * CSGP[3] * Θ.SGP * Cwₛ[1] * Cwₛ[3] + (-2) * CSGP[1] * CS[2,3] * Θ.SGP
    * Cwₛ[2] * Cwₛ[3] + CSGP[1] * CS[2,2] * Θ.SGP * Cwₛ[3]^2 + CS[1,3] * CSGP[2] * (CSGP[3]
    * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3]) + (-1) * CS[1,2] * CSGP[3] * (CSGP[3] * Cwₛ[2]
    + (-1) * CSGP[2] * Cwₛ[3]) + (-1) * CS[1,3] * Θ.SGP * Cwₛ[2] * (CSGP[3] * Cwₛ[2]
    + (-1) * CSGP[2] * Cwₛ[3]) + CS[1,2] * Θ.SGP * Cwₛ[3] * (CSGP[3] * Cwₛ[2] + (-1)
    * CSGP[2] * Cwₛ[3]) + (CSGP[3] * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3])^2 * ((CSGP[3]
    * Cwₛ[2] + (-1) * CSGP[2] * Cwₛ[3])^(-2) * ((-1) * CS[1,1] * CSGP[2]^2 * CS[3,3]
    + 2 * CS[1,1] * CS[2,3] * CSGP[2] * CSGP[3] + CS[1,2]^2 * CSGP[3]^2 + (-1) * CS[1,1]
    * CS[2,2] * CSGP[3]^2 + (-2) * CS[1,2] * CSGP[2] * CS[3,3] * Θ.SGP * Cwₛ[1] + 2 *
    CS[1,2] * CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[1] + CSGP[2]^2 * CS[3,3] * Θ.SG * Cwₛ[1]^2
    + ( -2) * CS[2,3] * CSGP[2] * CSGP[3] * Θ.SG * Cwₛ[1]^2 + CS[2,2] * CSGP[3]^2 * Θ.SG
    * Cwₛ[1]^2 + CS[2,3]^2 * Θ.SGP^2 * Cwₛ[1]^2 + (-1) * CS[2,2] * CS[3,3] * Θ.SGP^2
    * Cwₛ[1]^2 + 2 * CS[1,1] * CSGP[2] * CS[3,3] * Θ.SGP * Cwₛ[2] + (-2) * CS[1,1] *
    CS[2,3] * CSGP[3] * Θ.SGP * Cwₛ[2] + (-2) * CS[1,2] * CSGP[3]^2 * Θ.SG * Cwₛ[1] *
    Cwₛ[2] + 2 * CS[1,2] * CS[3,3] * Θ.SGP^2 * Cwₛ[1] * Cwₛ[2] + CS[1,1] * CSGP[3]^2
    * Θ.SG * Cwₛ[2]^2 + (-1) * CS[1,1] * CS[3,3] * Θ.SGP^2 * Cwₛ[2]^2 + CS[1,3]^2 * (CSGP[2]
    + (-1) * Θ.SGP * Cwₛ[2])^2 + (-2) * CS[1,1] * CS[2,3] * CSGP[2] * Θ.SGP * Cwₛ[3]
    + (-2) * CS[1,2]^2 * CSGP[3] * Θ.SGP * Cwₛ[3] + 2 * CS[1,1] * CS[2,2] * CSGP[3] *
    Θ.SGP * Cwₛ[3] + 2 * CS[1,2] * CSGP[2] * CSGP[3] * Θ.SG * Cwₛ[1] * Cwₛ[3] + (-2)
    * CS[1,2] * CS[2,3] * Θ.SGP^2 * Cwₛ[1] * Cwₛ[3] + (-2) * CS[1,1] * CSGP[2] * CSGP[3]
    * Θ.SG * Cwₛ[2] * Cwₛ[3] + 2 * CS[1,1] * CS[2,3] * Θ.SGP^2 * Cwₛ[2] * Cwₛ[3] + CS[1,1]
    * CSGP[2]^2 * Θ.SG * Cwₛ[3]^2 + CS[1,2]^2 * Θ.SGP^2 * Cwₛ[3]^2 + (-1) * CS[1,1] *
    CS[2,2] * Θ.SGP^2 * Cwₛ[3]^2 + CSGP[1]^2 * (CS[2,3]^2 + (-1) * CS[2,2] * CS[3,3]
    + CS[3,3] * Θ.SG * Cwₛ[2]^2 + (-2) * CS[2,3] * Θ.SG * Cwₛ[2] * Cwₛ[3] + CS[2,2] *
    Θ.SG * Cwₛ[3]^2) + (-2) * CSGP[1] * (CS[1,3] * ((-1) * CS[2,2] * CSGP[3] + CSGP[3]
    * Θ.SG * Cwₛ[2]^2 + CS[2,3] * (CSGP[2] + (-1) * Θ.SGP * Cwₛ[2]) + CS[2,2] * Θ.SGP
    * Cwₛ[3] + (-1) * CSGP[2] * Θ.SG * Cwₛ[2] * Cwₛ[3]) + CS[1,2] * ((-1) * CSGP[2] *
    CS[3,3] + CS[2,3] * CSGP[3] + CS[3,3] * Θ.SGP * Cwₛ[2] + (-1) * CS[2,3] * Θ.SGP *
    Cwₛ[3] + (-1) * CSGP[3] * Θ.SG * Cwₛ[2] * Cwₛ[3] + CSGP[2] * Θ.SG * Cwₛ[3]^2) + Cwₛ[1]
    * ( CS[2,3]^2 * Θ.SGP + (-1) * CS[2,2] * CS[3,3] * Θ.SGP + CSGP[2] * CS[3,3] * Θ.SG
    * Cwₛ[2] + CS[2,2] * CSGP[3] * Θ.SG * Cwₛ[3] + ( -1) * CS[2,3] * Θ.SG * (CSGP[3]
    * Cwₛ[2] + CSGP[2] * Cwₛ[3]))) + (-2) * CS[1,3] * (CS[1,2] * (CSGP[2] + (-1) * Θ.SGP
    * Cwₛ[2]) * (CSGP[3] + (-1) * Θ.SGP * Cwₛ[3]) + Cwₛ[1] * (CS[2,3] * Θ.SGP * ((-1)
    * CSGP[2] + Θ.SGP * Cwₛ[2]) + CSGP[2] * Θ.SG * ((-1) * CSGP[3] * Cwₛ[2] + CSGP[2]
    * Cwₛ[3]) + CS[2,2] * Θ.SGP * (CSGP[3] + (-1) * Θ.SGP * Cwₛ[3])))))^(1/2))
  )

  return ω₁
end

function getω₂(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64}, ω₁::Float64)::Float64
  ω₂::Float64 = (
    -((-SGP[3] + SGP[3] * ω₁ * TG.wₛ[1] + Θ.SGP * TG.wₛ[3] - SGP[1] * ω₁ * TG.wₛ[3]) /
    (SGP[3] * TG.wₛ[2] - SGP[2] * TG.wₛ[3]))
  )

  return ω₂
end

function getω₃(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64}, ω₁::Float64)::Float64
  ω₃::Float64 = (
    -((-SGP[2] + SGP[2] * ω₁ * TG.wₛ[1] + Θ.SGP * TG.wₛ[2] - SGP[1] * ω₁ * TG.wₛ[2]) /
    (-SGP[3] * TG.wₛ[2] + SGP[2] * TG.wₛ[3]))
  )

  return ω₃
end

#builds the scaling vector ω
function getω(Θ::Parameters, TG::TrisectedPortfolio, SGP::Vector{Float64})::Vector{Float64}

  ωcandidate1::Vector{Float64} = Vector{Float64}(undef, 3)
  ωcandidate2::Vector{Float64} = Vector{Float64}(undef, 3)
  ωcandidate1[1]::Float64 = getω₁C1(Θ, TG, SGP)
  ωcandidate2[1]::Float64 = getω₁C2(Θ, TG, SGP)

  #fill the rest of the candidate vectors
  ωcandidate1[2] = getω₂(Θ, TG, SGP, ωcandidate1[1])
  ωcandidate1[3] = getω₃(Θ, TG, SGP, ωcandidate1[1])

  ωcandidate2[2] = getω₂(Θ, TG, SGP, ωcandidate2[1])
  ωcandidate2[3] = getω₃(Θ, TG, SGP, ωcandidate2[1])

  errcandidate1::Float64 = abs(ωcandidate1' * SGP - Θ.SGP)
  errcandidate2::Float64 = abs(ωcandidate2' * SGP - Θ.SGP)

  #pick the vector with the lowest error
  ω::Vector{Float64} = errcandidate1 < errcandidate2 ? ωcandidate1 : ωcandidate2

  return ω
end

getω(Θ::Parameters, TG::TrisectedPortfolio, P::Portfolio)::Vector{Float64} =
  getω(Θ, TG, cov(TG, P))

function rescaletrisected!(TG::TrisectedPortfolio, ω::Vector{Float64})::TrisectedPortfolio
  wₜ::Vector{Float64} = TG.P.wₜ

  for i ∈ 1:length(wₜ)
    wₜ[i] *= ω[TG.assignments[i]]
  end

  return refreshtrisected!(TG, wₜ)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using SetRounding
using Test

@testset "setrounding(Float64)" begin
    a = 0.1; b = 0.3

    c = setrounding(Float64, RoundDown) do
        a + b
    end

    d = setrounding(Float64, RoundUp) do
        a + b
    end

    @test c == 0.39999999999999997
    @test d == 0.4
end

@testset "setrounding(Float32)" begin
    a = Float32(0.1)
    b = Float32(0.3)

    c = setrounding(Float32, RoundDown) do
       a + b
    end

    d = setrounding(Float32, RoundUp) do
       a + b
    end

    @test c == 0.4f0
    @test d == 0.40000004f0
end


text:
"""
    isvalidtree(x, b)

Check if `tree` is a valid tree for the signal `x`.

# Arguments
- `x::AbstractVector` or `x::AbstractMatrix`: Signal (1D or 2D).
- `b::BitVector`: Tree to represent `x`.

# Returns
- `::Bool`: Whether tree is valid.
"""
function Wavelets.Util.isvalidtree(x::AbstractMatrix, b::BitVector)
    n, m = size(x)
    nb = length(b)
    gettreelength(n,m) == nb || return false    # Tree does not have correct expected length

    L₀ = getdepth(nb, :quad)
    ns = (1<<(2*L₀)-1)÷3                        # Index to iterate to
    for i in 1:ns
        @inbounds isnode = b[i]
        @inbounds haschild = b[getchildindex(i,:topleft)] || b[getchildindex(i,:topright)] ||
                             b[getchildindex(i,:bottomleft)] || b[getchildindex(i,:bottomright)]
        if !isnode && haschild
            return false
        end
    end
    return true
end

"""
    getchildindex(idx, child)

Get the child index of a parent index `idx`.

# Arguments
- `idx::T where T<:Integer`: Index of parent node.
- `child::Symbol`: Type of child. For binary trees, available children are `:left` and
  `:right`. For quadtrees, available children are `:topleft`, `:topright`, `:bottomleft`,
  `:bottomright`.

# Returns
- `::T`: Index of child node.

# Examples
```repl
using WaveletsExt

getchildindex(3,:left) == 6
getchildindex(3,:right) == 7
getchildindex(3,:topleft) == 22
getchildindex(3,:topright) == 23
getchildindex(3,:bottomleft) == 24
getchildindex(3,:bottomright) == 25
```
"""
function getchildindex(idx::Integer, child::Symbol)
    @assert child ∈ [:left, :right, :topleft, :topright, :bottomleft, :bottomright]

    if child == :left
        return idx<<1
    elseif child == :right
        return idx<<1 + 1
    elseif child == :topleft
        return 4*idx - 2
    elseif child == :topright
        return 4*idx - 1
    elseif child == :bottomleft
        return 4*idx
    elseif child == :bottomright
        return 4*idx + 1
    else
        throw(ArgumentError("Invalid child $child."))
    end
end

"""
    getparentindex(idx, tree_type)

Get the parent index of the child node `idx`.

# Arguments
- `idx::T where T<:Integer`: Index of child node.
- `tree_type::Symbol`: Tree type, ie. `:binary` tree or `:quad` tree.

# Returns
- `::T`: Index of parent node.
"""
function getparentindex(idx::T, tree_type::Symbol) where T<:Integer
    @assert tree_type ∈ [:binary, :quad]

    if tree_type == :binary
        return idx>>1
    elseif tree_type == :quad
        return floor(T, (idx+2)/4)
    else
        throw(ArgumentError("Invalid tree type $tree_type."))
    end
end

# Get leaf nodes in the form of a BitVector
"""
    getleaf(tree, tree_type)

Returns the leaf nodes of a tree.

# Arguments
- `tree::BitVector`: BitVector to represent binary tree.

# Returns
`::BitVector`: BitVector that can represent a binary tree, but only the leaves are labeled
1, the rest of the nodes are labeled 0.

# Examples
```@repl
using Wavelets, WaveletsExt

tree = maketree(4, 2, :dwt)     # [1,1,0]
getleaf(tree)                   # [0,0,1,1,1,0,0]
```
"""
function getleaf(tree::BitVector, tree_type::Symbol)
    # Setup and Sanity check
    @assert tree_type ∈ [:binary, :quad]
    nₜ = length(tree)                                               # Tree length
    L₀ = getdepth(nₜ, tree_type)                                    # Max depth in tree
    Enₜ = tree_type == :binary ? 1<<(L₀+1)-1 : (1<<(2*L₀+2)-1)÷3    # Expected tree length
    @assert Enₜ == nₜ
    n = tree_type == :binary ? 1<<(L₀+1) : 1<<(2*L₀+2)              # Nodes in depth L₀+1
    ns = 1<<(L₀+1)
    x = tree_type == :binary ? trues(ns) : trues(ns,ns)             # Test signal
    @assert isvalidtree(x, tree)                                    # Check tree validity

    result = falses(n+nₜ)
    @inbounds result[1] = true
    for i in eachindex(tree)
        if !tree[i]             # i-th node does not exist
            continue
        else
            # Current node is not leaf
            @inbounds result[i] = false
            # Set children node as leaves
            if tree_type == :binary
                @inbounds result[getchildindex(i,:left)] = true
                @inbounds result[getchildindex(i,:right)] = true
            elseif tree_type == :quad
                @inbounds result[getchildindex(i,:topleft)] = true
                @inbounds result[getchildindex(i,:topright)] = true
                @inbounds result[getchildindex(i,:bottomleft)] = true
                @inbounds result[getchildindex(i,:bottomright)] = true
            else
                throw(ArgumentError("Invalid tree type $tree_type."))
            end
        end
    end
    return result
end

# Build a quadtree
"""
    maketree(x[, s])
    maketree(n, m, L[, s])

Build quadtree for 2D wavelet transform. Indexing of the tree are as follows:

```
Level 0                 Level 1                 Level 2             ...
_________________       _________________       _________________
|               |       |   2   |   3   |       |_6_|_7_|10_|11_|
|       1       |       |_______|_______|       |_8_|_9_|12_|13_|   ...
|               |       |   4   |   5   |       |14_|15_|18_|19_|
|_______________|       |_______|_______|       |16_|17_|20_|21_|
```

# Arguments
- `x::AbstractMatrix{T} where T<:Number`: Input array.
- `n::Integer`: Number of rows in `x`.
- `m::Integer`: Number of columns in `x`.
- `L::Integer`: Number of decomposition levels.
- `s::Symbol`: (Default: `:full`) Type of quadtree. Available types are `:full` and `:dwt`.

# Returns
`::BitVector`: Quadtree representation.

# Examples
```@repl
using Wavelets, WaveletsExt

x = randn(16,16)
maketree(x)
```
"""
function Wavelets.Util.maketree(x::AbstractMatrix{T}, s::Symbol = :full) where T<:Number
    return maketree(size(x,1), size(x,2), maxtransformlevels(x), s)
end

function Wavelets.Util.maketree(n::Integer, m::Integer, L::Integer, s::Symbol = :full)
    @assert s ∈ [:full, :dwt]
    L₀ = maxtransformlevels(min(n,m))
    @assert 0 ≤ L ≤ L₀
    nq = gettreelength(n,m)

    # Construct quadtree
    tree = falses(nq)

    # Fill in true values depending on input s
    if s == :full
        rng = 4 .^ (0:(L-1)) |> sum |> x -> 1:x
        for i in rng
            @inbounds tree[i] = true
        end
    elseif s == :dwt
        tree[1] = true     # Root node
        for i in 0:(L-2)
            idx = (1<<(2*i+2)+2)÷3  # Subsequent LL subspace nodes (= 4^(0:i) |> sum |> x -> x+1)
            @inbounds tree[idx] = true
        end
    else
        throw(ArgumentError("Unknown symbol."))
    end
    return tree
end

# Get level of a particular index in a tree
"""
    getdepth(idx, tree_type)

Get depth of `idx` in a binary tree or quadtree.

# Arguments
- `idx::T where T<:Integer`: Index of a node.
- `tree_type::Symbol`: Tree type. Supported types are `:binary` and `:quad` trees.

# Returns
`::T`: Depth of node `idx`.

# Examples
```@repl
using WaveletsExt

getdepth(1,:binary)     # 0
getdepth(3,:binary)     # 1
getdepth(8,:binary)     # 3

getdepth(1,:quad)       # 0
getdepth(3,:quad)       # 1
getdepth(8,:quad)       # 2
```

**See also:** [`Wavelets.Util.maketree`](@ref)
"""
function getdepth(idx::T, tree_type::Symbol) where T<:Integer
    @assert idx > 0
    @assert tree_type ∈ [:binary, :quad]

    if tree_type == :binary
        return floor(T, log2(idx))
    elseif tree_type == :quad
        return floor(T, log(4, 3*idx-2))
    else
        throw(ArgumentError("Unsupported tree type $tree_type."))
    end
end

"""
    gettreelength(n[, m])

Compute the length of the binary/quad tree `BitVector`.

# Arguments
- `n::T where T<:Integer`: Array size at 1st dimension/Signal length for 1D-signals.
- `m::T where T<:Integer`: Array size at 2nd dimension.

# Returns
- `::T`: Length of binary/quad tree.

# Examples
```repl
using WaveletsExt

Utils.gettreelength(8)          # 7
Utils.gettreelength(8,8)        # 21
```
"""
function gettreelength(n::T) where T<:Integer
    L = maxtransformlevels(n)
    return 1<<L-1
end

function gettreelength(n::T, m::T) where T<:Integer
    L = maxtransformlevels(min(n,m))
    return (1<<(2*L)-1)÷3
end


text:
using BenchmarkTools, JLD2, Serialization

const BACKEND = JLD2.MmapIO
# const BACKEND = IOStream
const TEMPFILE = begin
    tmp, io = mktemp()
    close(io)
    tmp
end

function jld_benchmark_write(x)
    f = jldopen(TEMPFILE, true, true, true, BACKEND)
    write(f, "x", x)
    close(f)
end

function jld_benchmark_read()
    f = jldopen(TEMPFILE, false, false, false, BACKEND)
    read(f, "x")
    close(f)
end

function serialize_benchmark_write(x)
    f = open(TEMPFILE, "w")
    serialize(f, x)
    close(f)
end

function serialize_benchmark_read()
    f = open(TEMPFILE)
    deserialize(f)
    close(f)
end

function bench(title, data)
    GC.gc()
    println(title, " JLD2 write")
    show(stdout, MIME("text/plain"), @benchmark jld_benchmark_write($data))
    println("\n")
    GC.gc()
    println(title, " JLD2 read")
    show(stdout, MIME("text/plain"), @benchmark jld_benchmark_read())
    println("\n")
    GC.gc()
    println(title, " Serialization.serialize")
    show(stdout, MIME("text/plain"), @benchmark serialize_benchmark_write($data))
    println("\n")
    GC.gc()
    println(title, " Serialization.deserialize")
    show(stdout, MIME("text/plain"), @benchmark serialize_benchmark_read())
    println("\n\n")
    isa(data, Array) && empty!(data)
end

# Vector{Int}
bench("int_vector", rand(typemin(Int):typemax(Int), 100000000))

# Vector{Any} of Ints
bench("any_int_vector", convert(Vector{Any}, rand(typemin(Int):typemax(Int), 1000000)))

# Vector{Integer} of Ints
bench("integer_int_vector", convert(Vector{Integer}, rand(typemin(Int):typemax(Int), 1000000)))

# Vector of non-builtin immutable
bench("complex128_vector", [complex(rand(), rand()) for i = 1:10000000])

# Vector{Any} of non-builtin immutable
struct IntWrapper
    x::Int
end
bench("any_complex128_vector", Any[complex(rand(), rand()) for i = 1:1000000])

# Many empty arrays
bench("empty_arrays", Any[Int[] for i = 1:1000000])

# Many one-eleemnt arrays
bench("one_element_arrays", [[rand(typemin(Int):typemax(Int))] for i = 1:1000000])

# Equivalent benchmark from https://github.com/timholy/HDF5.jl/issues/170
# This is 1/10 of the amount of data, but we don't have a problem with
# scaling, so this should be fine
mutable struct Cell
    a::Array{Float64,1}
    b::Array{Float64,1}
    c::Array{Float64,2}
    d::Float64
    e::Array{Float64,1}
end
bench("cell", [Cell(rand(3), rand(3), rand(8, 3), rand(), rand(3)) for i = 1:118138])


text:
macro optstruct(s)
  if s.head !== :struct
      error("Not a struct def")
  end

  ismutable = s.args[1]
  if !ismutable
      error("Not mutable")
  end

  name = s.args[2]
  body = s.args[3]

  ctor = :(function $(name)(;kwargs...)
              K = new()
              for (key, value) in kwargs
                  field_type_key = typeof(getfield(K,key))
                  setfield!(K, key, convert(field_type_key, value))
              end
              return K
          end)

  newbody = [body.args; ctor]

  return Expr(s.head, ismutable, name, Expr(body.head, newbody...))
end

function rootmeansquare(A)
  s = 0.0
  for a in A
     s += a * a
  end
  return sqrt(s / length(A))
end

function httpgetgzip(uri)
  r = HTTP.get(uri)
  io = TranscodingStream(GzipDecompressor(), IOBuffer(r.body))
  str = read(io, String)

  return str
end

function savetofile(str, path)
  f = open(path, "w")
  write(f, str)
  close(f)
end


text:
using Distributed
if nprocs() == 1
    addprocs(3)
end

function foo(arg1)
    print("foo: $arg1")
end

@everywhere function bar(arg1)
    print("foo: $arg1")
end

foo("1")
bar("2")


text:
################################################################################
# Meta Data
# This file is currently experimental and under development.

################################################################################
# SimpleDict
# Simple dictionaries are dictionaries whose keys are Symbols. Simple dicts can
# be used in argument processing if their keys are valid tokens.
const AbstractSimpleDict{T} = AbstractDict{Symbol, T} where {T}
const SimpleDict{T} = Dict{Symbol, T} where {T}
const PSimpleDict{T} = PDict{Symbol, T} where {T}
SimpleDict() = SimpleDict{Any}()
PSimpleDict() = PSimpleDict{Any}()
# How we parse arguments in simpleflat.
_simplebuild!(d::SimpleDict{T}, k::Symbol, v::S, args...) where {T,S} = begin
    d[k] = v
    return args
end
_simplebuild!(d::SimpleDict{T}, p::Pair{Symbol,S}, args...) where {T,S} = begin
    d[p[1]] = p[2]
    return args
end
_simplebuild!(d::SimpleDict{T}, p::Tuple{Symbol,S}) where {T,S} = begin
    d[p[1]] = p[2]
    return args
end
_simplebuild!(d::SimpleDict{T}, p::AbstractSimpleDict{S}, args...) where {T,S} = begin
    for (k,v) in p
        d[k] = v
    end
    return args
end
_simplebuild!(d::SimpleDict{T}, a0, args...) where {T} =
    error("invalid simpleflat argument type: $(typeof(a0))")
"""
    simpleflat(args...)
    simpleflat(T, args...)

Yields a simple dictionary that results from flattening the arguments. This
flattening occurs as follows:
 1. First, the arguments are parsed. The arguments may include simple
    dictionaries or paired key-value arguments (e.g., the following call is
    valid: `simpleflat(sym1, val1, dict, sym2, val2)`. These arguments are
    collapsed right-to-left. Additionally, Pair{Symbol,_} or Tuple{Symbol,_}
    objects may be passed in.
 2. Next, the keyword values are merged into the result.

The return type is always SimpleDict{Any} unless the type argument T is
supplied as the first argument, in which case, it is of type SimpleDict{T}.
If the supplied type T is Union{}, then the type used is the result of
typejoining across all the values.
"""
simpleflat(::Type{T0}, args...; kws...) where {T0} = begin
    tjoin = false
    T = T0 === Union{} ? Any : T0
    # Create the buffer dictionary into which we gather these.
    if length(args) > 0 && isa(args[1], AbstractSimpleDict)
        d = SimpleDict{T}(args[1])
        args = args[2:end]
    else
        d = SimpleDict{T}()
    end
    # Parse each argument using the above _simplebuild!() methods.
    aa = args
    while !isempty(aa)
        aa = _simplebuild!(d, aa...)
    end
    # Merge the keywords in last.
    for (k,v) in kws
        d[k] = v
    end
    # Check if we still need to join types.
    if T0 === Union{}
        T = typejoin([typeof(v) for v in values(d)])
        (T === Any) || return SimpleDict{T}(d)
    end
    return d
end
simpleflat(args...; kw...) = simpleflat(Union{}, args...; kw...)
simpleflat() = SimpleDict{Any}()

################################################################################
# MetaData
# Now that we've defined PDict, we can also fill out the meta-data concept for
# Air. Meta-data is an optional feature of a type that can be tracked and used
# Details for handling meta-data and declaring types that carry meta-data.
const AbstractMetaData = AbstractSimpleDict
const empty_metadata = PDict{Symbol,Any}()
# A trait for tracking if an object has meta-data or not.
"""
    MetaDataTrait

The MetaDataTrait allows types to indicate at compile time that their objects
support a meta-data dictionary. Meta-data dictionaries are Dict or PDict objects
that map symbols to anything.

A type T that supports meta-data will overload the following methods:
 * metatrait(T) yields a trait-object of type WithMetaData <: MetaDataTrait
 * metadata(::WithMetaData, t::T) yields the meta-data dictionary of t
 * Optionally, one of:
   * setmetadata!(::WithMetaData, t::T, dict), which changes the meta-data
     dictionary of to the dictionary dict, or
   * withmetadata(::WithMetaData, t::T, dict), which yields a duplicate of t
     with the new meta-data dictionary dict.
"""
abstract type MetaDataTrait end
struct WithMetaData <: MetaDataTrait end
struct WoutMetaData <: MetaDataTrait end
# #metatrait
"""
    metatrait(T)

Yields the meta-data trait (WithMetaData() or WoutMetaData(), both of abstract
type MetaDataTrait) of the type T. Types that yield WithMetaData() track a
dictionary whose type T <: AbstractMetaData.

Any type T that declares itself a meta-data type by overloading metatrait(T) to
return WithMetaData() must also implement the metadata() function, which yields
the meta-data dictionary of the object passed to it.
"""
metatrait(::Type{T}) where {T} = WoutMetaData() # Default: No metadata tracked.
# #metadata
"""
    metadata(t)

Yields the meta-data dictionary of the given object t. If t does not implement
meta-data (see ?MetaDataTrait) then an empty persistent dictionary is returned.
"""
metadata(t::T) where {T} = metadata(metatrait(T), t)
metadata(::WoutMetaData, t::T) where {T} = empty_metadata
# #meta
"""
    meta(object, key, default)
    meta(object, key)

Yields the value associated with the given meta-data key in the given object.
If the key is not found in the object's meta-data (or if the object does not
track meta-data) yields default. If default is not provided, then an error is
raised when it would be returned.
"""
meta(t::T, k::Symbol, dflt) where {T} = get(metadata(metatrait(T), t), k, dflt)
meta(t::T, k::Symbol) where {T} = begin
    md = metadata(metatrait(T), t)
    v = get(md, k, nothing)
    (v === nothing) && !in(k, md) && throw(KeyError(k))
    return v
end
# #setmetadata!
"""
    setmetadata!(object, dict)

Sets the meta-data dictionary for the giveb object to be equal to the given
dictionary dict. If the object is not an object that tracks its meta-data, then
an exception is thrown.

On success, the new value of the meta-data dictionary is returned.
"""
setmetadata!(t::T, d::D) where {T, D <: AbstractMetaData} =
    setmetadata!(metatrait(T), t, d)
setmetadata!(::WoutMetaData, t::T, d::D) where {T, D <: AbstractMetaData} =
    error("type $(T) does not support meta-data")
setmetadata!(::WithMetaData, t::T, d::D) where {T, D <: AbstractMetaData} =
    error("type $(T) does not support setting meta-data of type $(D)")
# #withmetadata
"""
    withmetadata(object, dict)

Yields a duplicate of the given object except that its meta-data dictionary
will have been updated to be equal to dict.
"""
withmetadata(t::T, d::D) where {T, D <: AbstractMetaData} =
    withmetadata(metatrait(T), t, d)
withmetadata(::WoutMetaData, t::T, d::D) where {T, D <: AbstractMetaData} =
    error("type $(T) does not support meta-data")
withmetadata(::WithMetaData, t::T, d::D) where {T, D <: AbstractMetaData} =
    error("type $(T) does not support clone-updating meta-data of type $(D)")
# #setmeta!
"""
    withmeta(obj, k, v)
    withmeta(obj, k1, v1, k2, v2, ...)
    withmeta(obj, <key>=v)
    withmeta(obj, <key1>=v1, <key2>=v2, ...)
    withmeta(obj, k1a, v1a, k2a, v2a, ..., <key1b>=v1b, <key2b>=v2b, ...)
    withmeta(T, obj, args...)

Yields a duplicate of the given object in which the given key or keys have been
updated to be associated with new values. In order for this function to work,
the object obj must support the `withmetadata()` method and it must have a
dictionary type that wortks with Air's `setindex()` and `delete()` methods,
otherwise behavior is undefined (likely an exception will be raised).

Note that matched `k, v` argument pairs may be replaced with or interspersed
with arguments such as `k => v` and `(k,v)` so long as the associated keys
are symbols.
"""
withmeta(t::T, args...; kw...) where {T} =
    withmeta(metatrait(T), t, args...; kw...)
withmeta(::WoutMetaData, t::T, args...; kw...) where {T} =
    error("type $(T) does not support meta-data")
_buildmeta(md0::AbstractSimpleDict{V}, args...; kw...) where {V} =
    simpleflat(V, md0, args...; kw...)
withmeta(Q::WithMetaData, t::T, args...; kw...) where {T} = begin
    # Get the meta-data that we start with
    md0 = metadata(Q, t)
    # Process the arguments
    md = _buildmeta(md0, args...; kw...)
    # Make sure the type is correct
    MD0 = typeof(md0)
    (MD0 === typeof(md)) || (md = MD0(md))
    # And return the updated version.
    return withmetadata(Q, t, md)
end
"""
    woutmeta(obj, k)
    woutmeta(obj, k1, k2...)

Yields a duplicate of the given object in which the given meta-data keys have
been removed. If obj does not support meta-data or does not support
clone-updating the meta-data, an exception will be thrown.
"""
woutmeta(t::T, args...) where {T} = woutmeta(metatrait(T), t, args...)
woutmeta(::WoutMetaData, t::T, args...) where {T} =
    error("type $(T) does not support meta-data")
_trimmeta(md0::AbstractSimpleDict{V}, args...) where {V} = begin
    for arg in args
        md0 = delete(md0, arg)
    end
    return md0
end
woutmeta(Q::WithMetaData, t::T, args...) where {T} = begin
    # Get the meta-data that we start with
    md0 = metadata(Q, t)
    # Process the arguments
    md = _trimmeta(md0, args...)
    # Make sure the type is correct
    MD0 = typeof(md0)
    (MD0 === typeof(md)) || (md = MD0(md))
    # And return the updated version.
    return withmetadata(Q, t, md)
end


text:
function _timeline()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.1)

    typeof(time_instance)==RheoTimeData && time_instance.t==(0.0:0.1:15.0)
end
@test _timeline()

function _frequencyspec()
    freq_instance = frequencyspec(ω_start=0.0, ω_end=15.0, step=0.1)

    typeof(freq_instance)==RheoFreqData && freq_instance.ω==(0.0:0.1:15.0)
end
@test _frequencyspec()

function _strainfunction()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.1)

    strain_imposed = strainfunction(time_instance, t -> t)

    typeof(strain_imposed)==RheoTimeData && strain_imposed.t==(0.0:0.1:15.0) && strain_imposed.ϵ==(0.0:0.1:15.0)
end
@test _strainfunction()

function _stressfunction()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.1)

    stress_imposed = stressfunction(time_instance, t -> t)

    typeof(stress_imposed)==RheoTimeData && stress_imposed.t==(0.0:0.1:15.0) && stress_imposed.σ==(0.0:0.1:15.0)
end
@test _stressfunction()

function _hstep()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.1)
    imposed = stressfunction(time_instance, hstep(offset=1.0, amp=2.0))

    all(v -> v==0, imposed.σ[1:10]) && all(v -> v==2.0, imposed.σ[11:end])
end
@test _hstep()

function _ramp()
    dt = 0.5
    _offset = 1.0
    time_instance = timeline(t_start=0.0, t_end=15.0, step=dt)
    imposed = stressfunction(time_instance, ramp(offset=_offset, gradient=1.0))

    all(v -> v==0, imposed.σ[imposed.t.<_offset]) && imposed.σ[imposed.t.==(_offset+dt)]≈[dt] && imposed.σ[end]≈14.0
end
@test _ramp()

function _stairs()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    test1 = all(v -> v==0.0, imposed.σ[imposed.t.<3])
    test2 = all(v -> v==0.5, imposed.σ[(imposed.t.>=3) .& (imposed.t.<7)])
    test3 = all(v -> v==1.0, imposed.σ[(imposed.t.>=7) .& (imposed.t.<11)])
    test4 = all(v -> v==1.5, imposed.σ[(imposed.t.>=11) .& (imposed.t.<15)])
    test5 = all(v -> v==2.0, imposed.σ[imposed.t.>=15])

    test1 && test2 && test3 && test4 && test5
end
@test _stairs()

function _square()
    dt = 0.1
    time_instance = timeline(t_start=0.0, t_end=3.0, step=dt)
    imposed = stressfunction(time_instance, square(offset=1.0, amp=1.0, period=1.0, width=0.5))

    testresult = all(v -> v==0, imposed.σ[imposed.t.<1.0])
    ons = [1.0, 1.1, 1.2, 1.3, 1.4, 2.0, 2.1, 2.2, 2.3, 2.4, 3.0]
    offs = [1.5, 1.6, 1.7, 1.8, 1.9, 2.5, 2.6, 2.7, 2.8, 2.9]
    for (i,t) in enumerate(imposed.t)
        if t in ons && imposed.σ==0.0
            testresult = testresult && false
        elseif t in offs && imposed.σ==1.0
            testresult = testresult && false
        end
    end
    testresult
end
@test _square()

function _sawtooth()
    time_instance = timeline(t_start=0.0, t_end=5.0, step=0.1)
    _offset = 1.0
    imposed = stressfunction(time_instance, sawtooth(offset=_offset, amp=1.0, period=1.0))
    # this test could be improved
    all(v -> v==0, imposed.σ[imposed.t.<_offset]) && all(v -> v==0.0, imposed.σ[imposed.t.%1.0.==0.0])
end
@test _sawtooth()

function _triangle()
    time_instance = timeline(t_start=0.0, t_end=5.0, step=0.1)
    _amp = 1.0
    _offset = 1.0
    imposed = stressfunction(time_instance, triangle(offset=_offset, amp=1.0, period=1.0, width=0.5))

    all(v -> v==0, imposed.σ[imposed.t.<_offset]) && all(v -> v==0.0, imposed.σ[imposed.t.%1.0.==0.0]) && all(v -> v==_amp, imposed.σ[((imposed.t.%1.0).==0.5) .& (imposed.t.>_offset)])
end
@test _triangle()


text:
using Distributions
export SisMcmcSampler
export SisMcmcInsertDelete, SisMcmcInsertDeleteGibbs, SisMcmcSplitMerge
export SisMcmcInsertDeleteCenter

abstract type SisMcmcSampler end

# From Models

# Insert-Delete
# -------------

struct SisMcmcInsertDeleteGibbs<: SisMcmcSampler
    ν_gibbs::Int   # Maximum number of edit ops
    ν_td::Int  # Maximum increase or decrease in number of interactions
    path_dist::PathDistribution  # Distribution used to introduce new interactions
    β::Real  # Extra probability of Gibbs move
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions in Gibbs scan
    ind_add::Vector{Int} # Storage for indexing of additions in Gibbs scan
    vals::Vector{Int} # Storage for valuse to insert in Gibbs scan
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SisMcmcInsertDeleteGibbs(
        path_dist::T;
        K=100,
        ν_gibbs=4, ν_td=2,  β=0.6,
        desired_samples=1000, lag=1, burn_in=0,
        init=InitMode()
        ) where {T<:PathDistribution}
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_gibbs)
        ind_add = zeros(Int, ν_gibbs)
        vals = zeros(Int, ν_gibbs)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_gibbs] = "(maximum number of edit operations in iMCMC-within-Gibbs conditional updates)"
        par_info[:ν_td] = "(maximum number of interaction insertions or deletions)"
        par_info[:path_dist] = "(path distribution for insertions)"
        par_info[:β] = "(probability of Gibbs scan)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_gibbs, ν_td, path_dist, β, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals, ind_td
            )
    end
end

SisMcmcInsertDeleteGibbs(
    model::SIS;
    K=100,
    ν=4, β=0.0,
    desired_samples=1000, lag=1, burn_in=0, init=InitMode()
    ) = SisMcmcInsertDeleteGibbs(
        PathPseudoUniform(model.V, TrGeometric(0.8, 1, model.K_inner));
        K=K, ν=ν, β=β, desired_samples=desired_samples, lag=lag, burn_in=burn_in
        )

function Base.show(io::IO, sampler::SisMcmcInsertDeleteGibbs)
    title = "MCMC Sampler for SIS Models via iMCMC-within-Gibbs and Interaction Insertion/Deletion."
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

struct SisMcmcInsertDelete <: SisMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of update move
    len_dist::DiscreteUnivariateDistribution  # Dist. to sample path lengths
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    dist_curr::Vector{Float64}
    function SisMcmcInsertDelete(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.4, len_dist=TrGeometric(0.8,1,K),
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        dist_curr = [0.0]
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            len_dist, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers,
            ind_update, ind_td_add, ind_td_del,
            dist_curr
            )
    end
end

function Base.show(io::IO, sampler::SisMcmcInsertDelete)
    title = "MCMC Sampler for SIS Model"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: interaction insertion/deletion with edit allocation updates.")
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

# Centered Insert/Delete
# ----------------------

struct SisMcmcInsertDeleteCenter <: SisMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of update move
    p::Float64  # Par for geometric on size difference
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SisMcmcInsertDeleteCenter(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.4, p=0.7,
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:p] = "(parameter for Geometric dists. on size difference from mean)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            p, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers,
            ind_update, ind_td
            )
    end
end

function Base.show(io::IO, sampler::SisMcmcInsertDeleteCenter)
    title = "MCMC Sampler for SIS Models via with Multinomial Allocated Updates and Interaction Insertion/Deletion."
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

# Merge/split
# -----------

struct SisMcmcSplitMerge <: SisMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of update move
    η::Float64  # Noisy parameter for split/merge
    p::Float64  # Parameter for Geometric dist in split/merge
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values to insert in interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage for location of split/merges
    p_ins::Geometric
    function SisMcmcSplitMerge(
        ;ν_ed=2, ν_td=2, β=0.7,
        η=0.7, p=0.7,
        K=100,
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        p_ins = Geometric(p)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:η] = "(noise parameter for merge/splits)"
        par_info[:p] = "(parameter for Geometric dist. in merge/splits)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of paths, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            η, p,
            K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td,
            p_ins
            )
    end
end

function Base.show(io::IO, sampler::SisMcmcSplitMerge)
    title = "MCMC Sampler for SIS Model"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: edit allocation with merge and split\n")
    println(io, "Parameters:")
    num_of_pars = 6
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end


text:
function QZSSAlmanac(arg0::JString, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jint, arg14::AbsoluteDate, arg2::jint, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return QZSSAlmanac((JString, jint, jdouble, jdouble, jdouble, jint, AbsoluteDate, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function QZSSAlmanac(arg0::JString, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jint, arg2::jint, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return QZSSAlmanac((JString, jint, jdouble, jdouble, jdouble, jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_af0(obj::QZSSAlmanac)
    return jcall(obj, "getAf0", jdouble, ())
end

function get_af1(obj::QZSSAlmanac)
    return jcall(obj, "getAf1", jdouble, ())
end

function get_cic(obj::QZSSAlmanac)
    return jcall(obj, "getCic", jdouble, ())
end

function get_cis(obj::QZSSAlmanac)
    return jcall(obj, "getCis", jdouble, ())
end

function get_crc(obj::QZSSAlmanac)
    return jcall(obj, "getCrc", jdouble, ())
end

function get_crs(obj::QZSSAlmanac)
    return jcall(obj, "getCrs", jdouble, ())
end

function get_cuc(obj::QZSSAlmanac)
    return jcall(obj, "getCuc", jdouble, ())
end

function get_cus(obj::QZSSAlmanac)
    return jcall(obj, "getCus", jdouble, ())
end

function get_date(obj::QZSSAlmanac)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::QZSSAlmanac)
    return jcall(obj, "getE", jdouble, ())
end

function get_health(obj::QZSSAlmanac)
    return jcall(obj, "getHealth", jint, ())
end

function get_i0(obj::QZSSAlmanac)
    return jcall(obj, "getI0", jdouble, ())
end

function get_i_dot(obj::QZSSAlmanac)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_iodc(obj::QZSSOrbitalElements)
    return jcall(obj, "getIODC", jint, ())
end

function get_iode(obj::QZSSOrbitalElements)
    return jcall(obj, "getIODE", jint, ())
end

function get_m0(obj::QZSSAlmanac)
    return jcall(obj, "getM0", jdouble, ())
end

function get_mean_motion(obj::QZSSAlmanac)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_omega0(obj::QZSSAlmanac)
    return jcall(obj, "getOmega0", jdouble, ())
end

function get_omega_dot(obj::QZSSAlmanac)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::QZSSAlmanac)
    return jcall(obj, "getPa", jdouble, ())
end

function get_prn(obj::QZSSAlmanac)
    return jcall(obj, "getPRN", jint, ())
end

function get_sma(obj::QZSSAlmanac)
    return jcall(obj, "getSma", jdouble, ())
end

function get_source(obj::QZSSAlmanac)
    return jcall(obj, "getSource", JString, ())
end

function get_tgd(obj::QZSSOrbitalElements)
    return jcall(obj, "getTGD", jdouble, ())
end

function get_time(obj::QZSSAlmanac)
    return jcall(obj, "getTime", jdouble, ())
end

function get_week(obj::QZSSAlmanac)
    return jcall(obj, "getWeek", jint, ())
end


text:
module CcallMacros
export @ccall, @cdef, @disable_sigint, @nonzero_systemerror

"""
`calltoccall` is an implementation detail of @ccall and @cdef

takes and expression like :(printf("%d"::Cstring, value::Cuint)::Cvoid)
returns: a tuple of (function_name, return_type, arg_types, args

The above input outputs this:
(:printf, :Cvoid, :((Cstring, Cuint)), ["%d", :value])

Note that the args are in an array, not a quote block and have to be
appended to the ccall in a separate step.
"""
function calltoccall(expr)
    expr.head != :(::) &&
        error("@ccall needs a function signature with a return type")
    rettype = expr.args[2]

    call = expr.args[1]
    call.head != :call &&
        error("@ccall has to be a function call")

    if (f = call.args[1]) isa Expr
        lib = f.args[1]
        fname = f.args[2]
        func = :(($fname, $lib))
    else
        func = QuoteNode(f)
    end
    argtypes = :(())
    args = []
    for arg in call.args[2:end]
        varargs = false
        if arg.head == :...
            varargs = true
            arg = arg.args[1]
        end

        arg.head != :(::) &&
            error("args in @ccall must be annotated")
        value = arg.args[1]
        type_ = arg.args[2]
        # This currently doesn't work.
        if varargs
            value = :($value...)
            type_ = :($type_...)
        end
        push!(args, value)
        push!(argtypes.args, type_)
    end
    func, rettype, argtypes, args
end

"""
convert a julia-style function definition to a ccall:

`@ccall printf("%d"::Cstring, 10::Cint)::Cvoid`

same as:

`ccall(:printf, Cvoid, (Cstring, Cint), "%d", 10)`
"""
macro ccall(expr)
    func, rettype, argtypes, args = calltoccall(expr)
    output = :(ccall($func, $rettype, $argtypes))
    append!(output.args, args)
    esc(output)
end

"""
define a _very_ thin wrapper function on a ccall. Mostly for wrapping
libraries quickly as a foundation for a higher-level interface.

@cdef mkfifo(path::Cstring, mode::Cuint)::Cint

becomes:

mkfifo(path, mode) = ccall(:mkfifo, Cint, (Cstring, Cuint), path, mode)
"""
macro cdef(expr)
    func, rettype, argtypes, args = calltoccall(expr)
    call = :(ccall($func, $rettype, $argtypes))
    append!(call.args, args)
    name = func isa QuoteNode ? func.value : func.args[1].value
    definition = :($name())
    append!(definition.args, args)
    esc(:($definition = $call))
end

"""
disable SIGINT while expr is being executed. Mostly useful for calling
C functions that call back into Julia in a concurrent context because
memory corruption can occur and crash the whole program.
"""
macro disable_sigint(expr)
    out = quote
        disable_sigint() do
            $expr
        end
    end
    esc(out)
end

const comment = r"#=.*?=# "

"""
throw a system error if the expression returns a non-zero exit status.
"""
macro nonzero_systemerror(expr)
    str = replace(string(expr), comment => "")
    out = quote
        err = $expr
        systemerror($str, err != 0)
        err
    end
    esc(out)
end

end # module


text:
module vision

include("features.jl")
include("visualization.jl")

end # module


text:
# using Revise, Plots, Test
using BifurcationKit, LinearAlgebra, Setfield, SparseArrays, ForwardDiff, Parameters
const BK = BifurcationKit
norminf = x -> norm(x, Inf)

Fbp(x, p) = [x[1] * (3.23 .* p.μ - p.x2 * x[1] + p.x3 * 0.234 * x[1]^2) + x[2], -x[2]]
par = (μ = -0.2, ν = 0, x2 = 1.12, x3 = 1.0)
####################################################################################################
opt_newton = NewtonPar(tol = 1e-9, maxIter = 20)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01, pMax = 0.4, pMin = -0.5, detectBifurcation = 3, nev = 2, newtonOptions = opt_newton, maxSteps = 100, nInversion = 4, tolBisectionEigenvalue = 1e-8, dsminBisection = 1e-9)

br, = continuation(
	Fbp, [0.1, 0.1], par, (@lens _.μ),
	opts_br; plot = false, verbosity = 0, normC = norminf, printSolution = (x, p) -> x[1])

####################################################################################################
# normal form computation
jet = BK.get3Jet(Fbp, (x, p) -> ForwardDiff.jacobian(z -> Fbp(z, p), x))

bp = BK.computeNormalForm(jet..., br, 1; verbose=false)
@test BK.isTranscritical(bp) == true

jet = BK.get3Jet(Fbp, (x, p) -> BK.finiteDifferences(z -> Fbp(z, p), x))
bp = BK.computeNormalForm(jet..., br, 1; verbose=false)
@test BK.isTranscritical(bp) == true
show(bp)

# normal form
nf = bp.nf

@test norm(nf.a) < 1e-10
	@test norm(nf.b1 - 3.23) < 1e-8
	@test norm(nf.b2/2 - -1.12) < 1e-10
	@test norm(nf.b3/6 - 0.234) < 1e-10

####################################################################################################
# same but when the eigenvalues are not saved in the branch but computed on the fly
br_noev, = BK.continuation(
	Fbp, [0.1, 0.1], par, (@lens _.μ),
	printSolution = (x, p) -> norminf(x),
	(@set opts_br.saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)
bp = BK.computeNormalForm(jet..., br_noev, 1; verbose=false)
nf = bp.nf
@test norm(nf[1]) < 1e-10
	@test norm(nf[2] - 3.23) < 1e-8
	@test norm(nf[3]/2 - -1.12) < 1e-10
	@test norm(nf[4]/6 - 0.234) < 1e-10
####################################################################################################
# Automatic branch switching
br2, = continuation(jet..., br, 1, setproperties(opts_br; pMax = 0.2, ds = 0.01, maxSteps = 14); printSolution = (x, p) -> x[1], verbosity = 0)
@test br2 isa Branch
@test BK.haseigenvalues(br2) == true
@test BK.haseigenvector(br2) == true
BK.eigenvals(br2, 1, true)
BK.getfirstusertype(br2)
@test length(br2) == 12
# plot(br,br2)

br3, = continuation(jet..., br, 1, setproperties(opts_br; ds = -0.01); printSolution = (x, p) -> x[1], verbosity = 0, usedeflation = true)
# plot(br,br3)

# automatic bifurcation diagram (Transcritical)
bdiag = bifurcationdiagram(jet..., [0.1, 0.1], par,  (@lens _.μ), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .5, ds = 0.01, dsmax = 0.05, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 30, newtonOptions = (@set opt_newton.verbose=false), maxSteps = 15);
	printSolution = (x, p) -> x[1],
	tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)
####################################################################################################
# Case of the pitchfork
par_pf = @set par.x2 = 0.0
par_pf = @set par_pf.x3 = -1.0
brp, = BK.continuation(
	Fbp, [0.1, 0.1], par_pf, (@lens _.μ),
	printSolution = (x, p) -> x[1],
	opts_br; plot = false, verbosity = 0, normC = norminf)
bpp = BK.computeNormalForm(jet..., brp, 1; verbose=false)
show(bpp)

nf = bpp.nf
@test norm(nf[1]) < 1e-9
	@test norm(nf[2] - 3.23) < 1e-9
	@test norm(nf[3]/2 - 0) < 1e-9
	@test norm(nf[4]/6 + 0.234) < 1e-9

# test automatic branch switching
br2, = continuation(jet..., brp, 1, setproperties(opts_br; maxSteps = 19, dsmax = 0.01, ds = 0.001, detectBifurcation = 2, newtonOptions = (@set opt_newton.verbose=false)); printSolution = (x, p) -> x[1], verbosity = 0, ampfactor = 1)
	# plot(brp,br2, marker=:d)

# test methods for aBS
BK.from(br2) |> BK.type
BK.from(br2) |> BK.istranscritical
BK.type(nothing)
BK.show(stdout, br2)
BK.propertynames(br2)

# automatic bifurcation diagram (Pitchfork)
bdiag = bifurcationdiagram(jet..., [0.1, 0.1], par_pf,  (@lens _.μ), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .5, ds = 0.01, dsmax = 0.05, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 30, newtonOptions = (@set opt_newton.verbose=false), maxSteps = 15);
	printSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 2, normC = norminf)
####################################################################################################
function Fbp2d(x, p)
	return [ p.α * x[1] * (3.23 .* p.μ - 0.123 * x[1]^2 - 0.234 * x[2]^2),
			 p.α * x[2] * (3.23 .* p.μ - 0.456 * x[1]^2 - 0.123 * x[2]^2),
			 -x[3]]
end

jet = BK.get3Jet(Fbp2d, (x, p) -> ForwardDiff.jacobian(z -> Fbp2d(z, p), x))

par = (μ = -0.2, ν = 0, α = -1)

for α in (1,1)
	global par = @set par.α = α
	br, = BK.continuation(
		Fbp2d, [0.01, 0.01, 0.01], par, (@lens _.μ),
		printSolution = (x, p) -> norminf(x),
		setproperties(opts_br; nInversion = 2); plot = false, verbosity = 0, normC = norminf)
	# we have to be careful to have the same basis as for Fbp2d or the NF will not match Fbp2d
	bp2d = BK.computeNormalForm(jet..., br, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
	show(bp2d)

	BK.nf(bp2d)
	length(bp2d)
	bp2d(rand(2), 0.2)
	bp2d(Val(:reducedForm), rand(2), 0.2)

	@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -par.α * 0.123) < 1e-10
	@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -par.α * 0.234) < 1e-10
	@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -par.α * 0.0)   < 1e-10
	@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -par.α * 0.456) < 1e-10
	@test norm(bp2d.nf.b2, Inf) < 3e-6
	@test norm(bp2d.nf.b1 - par.α * 3.23 * I, Inf) < 1e-9
	@test norm(bp2d.nf.a, Inf) < 1e-6
end

##############################
# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br_noev, = BK.continuation(
	Fbp2d, [0.01, 0.01, 0.01], par, (@lens _.μ),
	printSolution = (x, p) -> norminf(x),
	setproperties(opts_br; nInversion = 2, saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)
bp2d = BK.computeNormalForm(jet..., br_noev, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -0.123) < 1e-15
@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -0.234) < 1e-15
@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -0.0)   < 1e-15
@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -0.456) < 1e-15
@test norm(bp2d.nf.b2, Inf) < 3e-15
@test norm(bp2d.nf.b1 - 3.23 * I, Inf) < 1e-9
@test norm(bp2d.nf.a, Inf) < 1e-15
####################################################################################################
# vector field to test nearby secondary bifurcations
FbpSecBif(u, p) = @. -u * (p + u * (2-5u)) * (p -.15 - u * (2+20u))
dFbpSecBif(x,p)=  ForwardDiff.jacobian( z-> FbpSecBif(z,p), x)
jet = BK.get3Jet(FbpSecBif, dFbpSecBif)

br_snd1, = BK.continuation(
	FbpSecBif, [0.0], -0.2, (@lens _),
	printSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	setproperties(opts_br; pMin = -1.0, pMax = .3, ds = 0.001, dsmax = 0.005, nInversion = 8, detectBifurcation=3); plot = false, verbosity = 0, normC = norminf)

# plot(br_snd1)

br_snd2, = BK.continuation(
	jet..., br_snd1, 1,
	printSolution = (x, p) -> x[1],
	setproperties(opts_br; pMin = -1.2, pMax = 0.2, ds = 0.001, detectBifurcation = 3, maxSteps=19, nInversion = 8, newtonOptions = NewtonPar(opts_br.newtonOptions),dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20); plot = false, verbosity = 0, normC = norminf,
	# tangentAlgo = BorderedPred(),
	# finaliseSolution = (z, tau, step, contResult) ->
	# 	(Base.display(contResult.eig[end].eigenvals) ;true)
	)

	# plot(plot(br_snd2.branch[1,:] |> diff, marker = :d),
	# plot(br_snd1,br_snd2, putbifptlegend=false))

bdiag = bifurcationdiagram(jet..., [0.0], -0.2, (@lens _), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .3, ds = 0.001, dsmax = 0.005, nInversion = 8, detectBifurcation = 3,dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));
	printSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)

# plot(bdiag; putbifptlegend=false, markersize=2, plotfold=false, title = "#branch = $(size(bdiag))")

# test calls for aBD
BK.level(bdiag)
BK.hasbranch(bdiag)
BK.from(bdiag.child[1].γ)
BK.getBranchesFromBP(bdiag, 2)
BK.getContResult(br)
BK.getContResult(getBranch(bdiag,(1,)).γ)
size(bdiag)
getBranch(bdiag, (1,))
show(stdout, bdiag)
####################################################################################################
# test of the pitchfork-D6 normal form
function FbpD6(x, p)
	return [ p.μ * x[1] + (p.a * x[2] * x[3] - p.b * x[1]^3 - p.c*(x[2]^2 + x[3]^2) * x[1]),
			 p.μ * x[2] + (p.a * x[1] * x[3] - p.b * x[2]^3 - p.c*(x[3]^2 + x[1]^2) * x[2]),
			 p.μ * x[3] + (p.a * x[1] * x[2] - p.b * x[3]^3 - p.c*(x[2]^2 + x[1]^2) * x[3])]
end
jet = BK.get3Jet(FbpD6, (x, p) -> ForwardDiff.jacobian(z -> FbpD6(z, p), x))

pard6 = (μ = -0.2, a = 0.3, b = 1.5, c = 2.9)

br, = BK.continuation(
	FbpD6, zeros(3), pard6, (@lens _.μ),
	printSolution = (x, p) -> norminf(x),
	setproperties(opts_br; nInversion = 6, ds = 0.001); plot = false, verbosity = 0, normC = norminf)

# plot(br;  plotfold = false)
# we have to be careful to have the same basis as for Fbp2d or the NF will not match Fbp2d
bp2d = BK.computeNormalForm(jet..., br, 1; ζs = [[1, 0, 0.], [0, 1, 0.], [0, 0, 1.]])
BK.nf(bp2d)

@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -pard6.b) < 1e-10
	@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -pard6.c) < 1e-10
	@test abs(bp2d.nf.b2[1,2,3] - pard6.a)   < 1e-10

# test the evaluation of the normal form
x0 = rand(3); @test norm(FbpD6(x0, set(pard6, br.lens, 0.001))  - bp2d(Val(:reducedForm), x0, 0.001), Inf) < 1e-12

br1, = BK.continuation(
	jet..., br, 1,
	setproperties(opts_br; nInversion = 4, dsmax = 0.005, ds = 0.001, maxSteps = 300, pMax = 1.); plot = false, verbosity = 0, normC = norminf, printSolution = (x, p) -> norminf(x))
	# plot(br1..., br, plotfold=false, putbifptlegend=false)

bdiag = bifurcationdiagram(jet..., zeros(3), pard6, (@lens _.μ), 3,
	(args...) -> setproperties(opts_br; pMin = -0.250, pMax = .4, ds = 0.001, dsmax = 0.005, nInversion = 4, detectBifurcation = 3, dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));
	printSolution = (x, p) -> norminf(x),
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)

# plot(bdiag; putbifptlegend=false, markersize=2,plotfold=false);title!("#branch = $(size(bdiag))")

####################################################################################################
# test of the Hopf normal form
function Fsl2!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]
	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

Fsl2(x, p) = Fsl2!(similar(x), x, p, 0.)
par_sl = (r = 0.5, μ = 0.132, ν = 1.0, c3 = 1.123, c5 = 0.2)
jet = BK.get3Jet(Fsl2, (x, p) -> ForwardDiff.jacobian(z -> Fsl2(z, p), x))

# detect hopf bifurcation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.01, pMax = 0.1, pMin = -0.3, detectBifurcation = 3, nev = 2, newtonOptions = (@set opt_newton.verbose = false), maxSteps = 100)

br, = BK.continuation(
	Fsl2, [0.0, 0.0], (@set par_sl.r = -0.1), (@lens _.r),
	printSolution = (x, p) -> norminf(x),
	opts_br; plot = false, verbosity = 0, normC = norminf)

hp = BK.computeNormalForm(jet..., br, 1)

nf = hp.nf
BK.type(hp)

@test abs(nf.a - 1) < 1e-9
@test abs(nf.b/2 - (-par_sl.c3 + im*par_sl.μ)) < 1e-14

##############################
# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br, _ = BK.continuation(
	Fsl2, [0.0, 0.0], (@set par_sl.r = -0.1), (@lens _.r),
	printSolution = (x, p) -> norminf(x),
	setproperties(opts_br, saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)

hp = BK.computeNormalForm(jet..., br, 1)

nf = hp.nf

@test abs(nf.a - 1) < 1e-9
@test abs(nf.b/2 - (-par_sl.c3 + im*par_sl.μ)) < 1e-14

show(hp)


text:
using PrettyTables

header = ["Return value" "DispatchedTuple" "DispatchedSet"
    "" "(non-unique keys allowed)" "(unique keys only)"]

col1 = ["Type", "Unregistered key (without default)", "Unregistered key (with default)", "Duplicative key"]

DT = ["Tuple", "()", "(default,)", "all registered values"]

DTS = ["Value", "error", "default", "one value"]

data = hcat(col1, DT, DTS)
pretty_table(
    data,
    header,
    header_crayon = crayon"yellow bold",
    crop = :none,
)


text:
struct NoDiag <: Diagnostics end

"""
Look for the Apparent Horizon during the evolution
"""
Base.@kwdef struct DiagAH <: Diagnostics
    # negative values turns off
    find_AH_every       :: Int  = -1
    find_AH_every_t     :: Float32  = -1.0
end

function (::NoDiag)(bulkevols, bulkconstrains, bulkderivs, boundary::Boundary,
                    gauge::Gauge, horizoncache::HorizonCache, systems::SystemPartition,
                    tinfo::Jecco.TimeInfo, evoleq::EvolutionEquations, io::InOut)
    # no-op
    () -> nothing
end

# initial guess for AH surface
function sigma_guess!(sigma, gaugecondition::ConstantAH)
    uAH = gaugecondition.u_AH
    fill!(sigma, 1/uAH)
end

function sigma_guess!(sigma, gaugecondition::GaugeCondition)
    # for other gauge conditions, just initialize sigma to 1, for lack of better
    # alternative
    fill!(sigma, 1)
end

function (diagnostics::DiagAH)(
    bulkevols, bulkconstrains, bulkderivs, boundary::Boundary, gauge::Gauge,
    horizoncache::HorizonCache, systems::SystemPartition, tinfo::Jecco.TimeInfo,
    evoleq::EvolutionEquations, io::InOut)

    # Apparent Horizon surface
    sigma = similar(gauge.xi)

    empty   = Cartesian{1}("u", systems[1].ucoord[1], systems[1].ucoord[1], 1)
    chart2D = Chart(empty, systems[end].xcoord, systems[end].ycoord)

    out_diag  = Jecco.Output(io.out_dir, "diagnostics_", tinfo)

    fields = Jecco.Field("sigma", sigma, chart2D)

    last_output_AH = -diagnostics.find_AH_every_t

    function ()
        it = tinfo.it
        tt = tinfo.t

        do_find_AH = false

        if (diagnostics.find_AH_every > 0 && it % diagnostics.find_AH_every == 0)
            do_find_AH = true
        end

        if (diagnostics.find_AH_every_t > 0 &&
            tt >= last_output_AH + diagnostics.find_AH_every_t - 1e-12)
            do_find_AH     = true
            last_output_AH = tt
        end

        if do_find_AH
            # it's useful to reset the initial guess, as the finder may sometimes fail
            # to find the AH, and the sigma surface would then be very displaced.
            sigma_guess!(sigma, evoleq.gaugecondition)
            find_AH!(sigma, bulkconstrains[end], bulkevols[end], bulkderivs[end], gauge,
                     horizoncache, systems[end], evoleq.ahf)

            # write sigma output
            out_diag(fields)
        end

        nothing
    end
end