text:
export reedsshepp, reedsshepp_length, reedsshepp_waypoints

reedsshepp_length(q0::StaticVector{3}, qf::StaticVector{3}; r=1) = reedsshepp(q0, qf, r=r).cost
function reedsshepp_waypoints(q0::StaticVector{3}, qf::StaticVector{3}, dt_or_N; v=1, r=1)
    waypoints(SimpleCarDynamics{0,0}(), SE2State(q0), reedsshepp(q0, qf, v=v, r=r).controls, dt_or_N)
end

const POST, POST_T, POST_R, POST_B, POST_R_T, POST_B_T, POST_B_R, POST_B_R_T = 0, 1, 2, 3, 4, 5, 6, 7

function reedsshepp((x0, y0, θ0)::StaticVector{3,T0}, (xf, yf, θf)::StaticVector{3,Tf}; v=1, r::R=1) where {T0,Tf,R}
    T = promote_type(T0, Tf, R)
    dx = (xf - x0)/r
    dy = (yf - y0)/r
    st, ct = sincos(θ0)
    target = SE2State{T}(dx*ct + dy*st, -dx*st + dy*ct, mod2piF(θf - θ0))

    tTarget   = timeflip(target)
    rTarget   = reflect(target)
    trTarget  = reflect(tTarget)
    bTarget   = backwards(target)
    btTarget  = timeflip(bTarget)
    brTarget  = reflect(bTarget)
    btrTarget = reflect(btTarget)

    c, ctrl, post = T(Inf), zeros(SVector{5,VelocityCurvatureStep{T}}), POST
    # (8.1) C S C
    b, c, ctrl = LpSpLp(target,       c, ctrl); b && (post = POST)
    b, c, ctrl = LpSpLp(tTarget,      c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpSpLp(rTarget,      c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpSpLp(trTarget,     c, ctrl); b && (post = POST_R_T)

    # (8.2) C S C
    b, c, ctrl = LpSpRp(target,       c, ctrl); b && (post = POST)
    b, c, ctrl = LpSpRp(tTarget,      c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpSpRp(rTarget,      c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpSpRp(trTarget,     c, ctrl); b && (post = POST_R_T)

    # (8.3) C|C|C
    b, c, ctrl = LpRmLp(target,       c, ctrl); b && (post = POST)
    # b, c, ctrl = LpRmLp(tTarget,      c, ctrl); b && (post = POST_T) # (redundant)
    b, c, ctrl = LpRmLp(rTarget,      c, ctrl); b && (post = POST_R)
    # b, c, ctrl = LpRmLp(trTarget,     c, ctrl); b && (post = POST_R_T) # (redundant)

    # (8.4) C|C C
    b, c, ctrl = LpRmLm(target,       c, ctrl); b && (post = POST)
    b, c, ctrl = LpRmLm(tTarget,      c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRmLm(rTarget,      c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRmLm(trTarget,     c, ctrl); b && (post = POST_R_T)
    b, c, ctrl = LpRmLm(bTarget,      c, ctrl); b && (post = POST_B)
    b, c, ctrl = LpRmLm(btTarget,     c, ctrl); b && (post = POST_B_T)
    b, c, ctrl = LpRmLm(brTarget,     c, ctrl); b && (post = POST_B_R)
    b, c, ctrl = LpRmLm(btrTarget,    c, ctrl); b && (post = POST_B_R_T)

    # (8.7) C Cu|Cu C
    b, c, ctrl = LpRpuLmuRm(target,   c, ctrl); b && (post = POST)
    b, c, ctrl = LpRpuLmuRm(tTarget,  c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRpuLmuRm(rTarget,  c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRpuLmuRm(trTarget, c, ctrl); b && (post = POST_R_T)

    # (8.8) C|Cu Cu|C
    b, c, ctrl = LpRmuLmuRp(target,   c, ctrl); b && (post = POST)
    b, c, ctrl = LpRmuLmuRp(tTarget,  c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRmuLmuRp(rTarget,  c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRmuLmuRp(trTarget, c, ctrl); b && (post = POST_R_T)

    # (8.9)
    b, c, ctrl = LpRmSmLm(target,     c, ctrl); b && (post = POST)
    b, c, ctrl = LpRmSmLm(tTarget,    c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRmSmLm(rTarget,    c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRmSmLm(trTarget,   c, ctrl); b && (post = POST_R_T)
    b, c, ctrl = LpRmSmLm(bTarget,    c, ctrl); b && (post = POST_B)
    b, c, ctrl = LpRmSmLm(btTarget,   c, ctrl); b && (post = POST_B_T)
    b, c, ctrl = LpRmSmLm(brTarget,   c, ctrl); b && (post = POST_B_R)
    b, c, ctrl = LpRmSmLm(btrTarget,  c, ctrl); b && (post = POST_B_R_T)

    # (8.10)
    b, c, ctrl = LpRmSmRm(target,     c, ctrl); b && (post = POST)
    b, c, ctrl = LpRmSmRm(tTarget,    c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRmSmRm(rTarget,    c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRmSmRm(trTarget,   c, ctrl); b && (post = POST_R_T)
    b, c, ctrl = LpRmSmRm(bTarget,    c, ctrl); b && (post = POST_B)
    b, c, ctrl = LpRmSmRm(btTarget,   c, ctrl); b && (post = POST_B_T)
    b, c, ctrl = LpRmSmRm(brTarget,   c, ctrl); b && (post = POST_B_R)
    b, c, ctrl = LpRmSmRm(btrTarget,  c, ctrl); b && (post = POST_B_R_T)

    # (8.11) C|Cpi/2 S Cpi/2|C
    b, c, ctrl = LpRmSmLmRp(target,   c, ctrl); b && (post = POST)
    b, c, ctrl = LpRmSmLmRp(tTarget,  c, ctrl); b && (post = POST_T)
    b, c, ctrl = LpRmSmLmRp(rTarget,  c, ctrl); b && (post = POST_R)
    b, c, ctrl = LpRmSmLmRp(trTarget, c, ctrl); b && (post = POST_R_T)

    ctrl = scalespeed.(scaleradius.(ctrl, r), v)
    if post == POST_T
        ctrl = timeflip.(ctrl)
    elseif post == POST_R
        ctrl = reflect.(ctrl)
    elseif post == POST_B
        ctrl = reverse(ctrl)
    elseif post == POST_R_T
        ctrl = reflect.(timeflip.(ctrl))
    elseif post == POST_B_T
        ctrl = timeflip.(ctrl)
        ctrl = reverse(ctrl)
    elseif post == POST_B_R
        ctrl = reflect.(ctrl)
        ctrl = reverse(ctrl)
    elseif post == POST_B_R_T
        ctrl = reflect.(timeflip.(ctrl))
        ctrl = reverse(ctrl)
    end
    (cost=c*r/v, controls=ctrl)
end

# Utilities (pedantic about typing to guard against problems)
@inline R(x::T, y::T) where {T} = hypot(x, y), atan(y, x)
@inline function M(t::T) where {T}
    m = mod2piF(t)
    ifelse(m > pi, m - 2*T(pi), m)
end
@inline function Tau(u::T, v::T, E::T, N::T) where {T}
    delta = u - v
    A = sin(u) - sin(delta)
    B = cos(u) - cos(delta) - 1
    r, θ = R(E*A + N*B, N*A - E*B)
    t = 2*cos(delta) - 2*cos(v) - 2*cos(u) + 3
    M(θ + ifelse(t < 0, pi, 0))
end
@inline Omega(u::T, v::T, E::T, N::T, t::T) where {T} = M(Tau(u, v, E, N) - u + v - t)
@inline timeflip(q::SE2State{T}) where {T} = SE2State(-q.x, q.y, -q.θ)
@inline reflect(q::SE2State{T}) where {T} = SE2State(q.x, -q.y, -q.θ)
@inline backwards(q::SE2State{T}) where {T} = SE2State(q.x*cos(q.θ) + q.y*sin(q.θ), q.x*sin(q.θ) - q.y*cos(q.θ), q.θ)

@inline function LpSpLp((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    r, θ = R(tx - sin(tθ), ty - 1 + cos(tθ))
    u = r
    t = mod2piF(θ)
    v = mod2piF(tθ - t)
    cnew = t + u + v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol(1, t),
        carsegment2stepcontrol(0, u),
        carsegment2stepcontrol(1, v),
        zero(VelocityCurvatureStep{T}),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpSpRp((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    r, θ = R(tx + sin(tθ), ty - 1 - cos(tθ))
    r*r < 4 && return false, c, ctrl
    u = sqrt(r*r - 4)
    r1, θ1 = R(u, T(2))
    t = mod2piF(θ + θ1)
    v = mod2piF(t - tθ)
    cnew = t + u + v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol( 0, u),
        carsegment2stepcontrol(-1, v),
        zero(VelocityCurvatureStep{T}),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmLp((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx - sin(tθ)
    N = ty + cos(tθ) - 1
    E*E + N*N > 16 && return false, c, ctrl
    r, θ = R(E, N)
    u = acos(1 - r*r/8)
    t = mod2piF(θ - u/2 + pi)
    v = mod2piF(pi - u/2 - θ + tθ)
    u = -u
    cnew = t - u + v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, u),
        carsegment2stepcontrol( 1, v),
        zero(VelocityCurvatureStep{T}),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmLm((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx - sin(tθ)
    N = ty + cos(tθ) - 1
    E*E + N*N > 16 && return false, c, ctrl
    r, θ = R(E, N)
    u = acos(1 - r*r/8)
    t = mod2piF(θ - u/2 + pi)
    v = mod2piF(pi - u/2 - θ + tθ) - 2*T(pi)
    u = -u
    cnew = t - u - v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, u),
        carsegment2stepcontrol( 1, v),
        zero(VelocityCurvatureStep{T}),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRpuLmuRm((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx + sin(tθ)
    N = ty - cos(tθ) - 1
    p = (2 + sqrt(E*E + N*N))/4
    (p < 0 || p > 1) && return false, c, ctrl
    u = acos(p)
    t = mod2piF(Tau(u, -u, E, N))
    v = mod2piF(Omega(u, -u, E, N, tθ)) - 2*T(pi)
    cnew = t + 2*u - v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, u),
        carsegment2stepcontrol( 1, -u),
        carsegment2stepcontrol(-1, v),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmuLmuRp((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx + sin(tθ)
    N = ty - cos(tθ) - 1
    p = (20 - E*E - N*N)/16
    (p < 0 || p > 1) && return false, c, ctrl
    u = -acos(p)
    t = mod2piF(Tau(u, u, E, N))
    v = mod2piF(Omega(u, u, E, N, tθ))
    cnew = t - 2*u + v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, u),
        carsegment2stepcontrol( 1, u),
        carsegment2stepcontrol(-1, v),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmSmLm((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx - sin(tθ)
    N = ty + cos(tθ) - 1
    D, β = R(E, N)
    D < 2 && return false, c, ctrl
    γ = acos(2/D)
    F = sqrt(D*D/4 - 1)
    t = mod2piF(pi + β - γ)
    u = 2 - 2*F
    u > 0 && return false, c, ctrl
    v = mod2piF(-3*T(pi)/2 + γ + tθ - β) - 2*T(pi)
    cnew = t + T(pi)/2 - u - v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, -T(pi)/2),
        carsegment2stepcontrol( 0, u),
        carsegment2stepcontrol( 1, v),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmSmRm((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx + sin(tθ)
    N = ty - cos(tθ) - 1
    D, β = R(E, N)
    D < 2 && return false, c, ctrl
    t = mod2piF(β + T(pi)/2)
    u = 2 - D
    u > 0 && return false, c, ctrl
    v = mod2piF(-T(pi) - tθ + β) - 2*T(pi)
    cnew = t + T(pi)/2 - u - v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, -T(pi)/2),
        carsegment2stepcontrol( 0, u),
        carsegment2stepcontrol(-1, v),
        zero(VelocityCurvatureStep{T})
    )
    true, cnew, ctrl
end

@inline function LpRmSmLmRp((tx, ty, tθ)::SE2State{T}, c::T, ctrl::SVector{5,VelocityCurvatureStep{T}}) where {T}
    E = tx + sin(tθ)
    N = ty - cos(tθ) - 1
    D, β = R(E, N)
    D < 2 && return false, c, ctrl
    γ = acos(2/D)
    F = sqrt(D*D/4 - 1)
    t = mod2piF(pi + β - γ)
    u = 4 - 2*F
    u > 0 && return false, c, ctrl
    v = mod2piF(pi + β - tθ - γ)
    cnew = t + pi - u + v
    c <= cnew && return false, c, ctrl
    ctrl = SVector(
        carsegment2stepcontrol( 1, t),
        carsegment2stepcontrol(-1, -T(pi)/2),
        carsegment2stepcontrol( 0, u),
        carsegment2stepcontrol( 1, -T(pi)/2),
        carsegment2stepcontrol(-1, v)
    )
    true, cnew, ctrl
end


text:
#initialize hyper-parameters of ground state and dynamic optimizations

#define structure type containing hyper-parameters of ground state optimization:
#number of Monte Carlo samples 'samples'
#number of iteration steps 'n_iter'
#size of the learning rate 'step_size'
struct GS_hp
    samples::Int64
    n_iter::Int64
    step_size::Float64
end

#define structure type containing hyper-parameters of dynamic optimization:
#number of Monte Carlo samples 'samples'
#total integration time 'tot_time'
#step size of the time-integration scheme 'step_size'
struct DYN_hp
    samples::Int64
    tot_time::Float64
    step_size::Float64
    minres_::Bool
end

ext: jl

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
export zoom_Sec

function zoom_Sec(h :: AbstractNLPModel;
				  kwargs...)

	(t, fₖ, normg, iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess) = trouve_intervalleA(h; direction = "Sec", kwargs...)
	return (t, fₖ, normg, iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)
end


text:
"""

The exported names are
$(EXPORTS)
"""
module ReferenceFEs

using Test
using DocStringExtensions
using LinearAlgebra
using Combinatorics

using Gridap.Helpers
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials
using Gridap.Integration

import Gridap.Arrays: kernel_cache
import Gridap.Arrays: apply_kernel!
import Gridap.Arrays: kernel_return_type
import Gridap.Fields: evaluate
import Gridap.Polynomials: MonomialBasis

import Gridap.Polynomials: get_order
import Gridap.Polynomials: get_orders

import Gridap.Integration: Quadrature
import Gridap.Integration: num_dims
import Gridap.Integration: num_point_dims

import Gridap.Io: to_dict
import Gridap.Io: from_dict

import Base: ==

export Polytope
export ExtrusionPolytope
export get_extrusion
export get_faces
export get_dimranges
export get_dimrange
export get_vertex_coordinates
export get_facet_normals
export get_facet_orientations
export get_edge_tangents
export get_vertex_permutations
export get_face_dimranges
export get_face_coordinates
export is_simplex
export is_n_cube
export simplexify
export num_dims
export num_cell_dims
export num_point_dims
export num_vertices
export num_faces
export num_facets
export num_edges
export get_facedims
export get_offsets
export get_offset
export get_face_vertices
export get_reffaces
export get_face_type
export get_bounding_box
export get_face_vertex_permutations
export test_polytope
export VERTEX
export SEGMENT
export TRI
export QUAD
export TET
export HEX
export WEDGE
export PYRAMID
export HEX_AXIS
export TET_AXIS
export INVALID_PERM

export Dof
export get_nodes
export get_face_moments
export get_face_nodes_dofs
export get_nodes
export evaluate_dof!
export evaluate_dof
export dof_cache
export dof_return_type
export test_dof
# export evaluate_dof_array

export ReferenceFE
export GenericRefFE
export get_polytope
export get_prebasis
export get_dof_basis
export get_default_conformity
export get_face_own_dofs
export get_face_own_dofs_permutations
export get_face_dofs
export get_own_dofs_permutations
export get_shapefuns
export compute_shapefuns
export test_reference_fe
export num_dofs

export LagrangianRefFE
export GenericLagrangianRefFE
export get_face_own_nodes_permutations
export get_own_nodes_permutations
export get_node_coordinates
export get_dof_to_node
export get_dof_to_comp
export get_node_and_comp_to_dof
export get_vertex_node
export num_nodes
export test_lagrangian_reference_fe

export LagrangianDofBasis
export compute_monomial_basis
export compute_own_nodes
export compute_face_orders
export compute_nodes
export compute_own_nodes_permutations
export compute_lagrangian_reffaces
export is_first_order
export is_affine
export is_Q
export is_P
export is_S

export MomentBasedDofBasis
export get_face_own_nodes
export get_face_nodes

export VERTEX1
export SEG2
export TRI3
export QUAD4
export TET4
export HEX8

export CONT
export DISC

export Conformity
export L2Conformity
export H1Conformity
export GradConformity
export CurlConformity
export DivConformity
export CDConformity

export SerendipityRefFE
export RaviartThomasRefFE
export NedelecRefFE

include("Polytopes.jl")

include("ExtrusionPolytopes.jl")

include("Dofs.jl")

include("MockDofs.jl")

include("LagrangianDofBases.jl")

include("ReferenceFEInterfaces.jl")

include("LagrangianRefFEs.jl")

include("CLagrangianRefFEs.jl")

include("SerendipityRefFEs.jl")

include("PDiscRefFEs.jl")

include("CDLagrangianRefFEs.jl")

include("RaviartThomasRefFEs.jl")

include("NedelecRefFEs.jl")

end # module


text:
using Documenter, DemoCards
using ImageFiltering

demos, demos_cb, demos_assets = makedemos("demos")

assets = []
isnothing(demos_assets) || (push!(assets, demos_assets))
format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = assets)

makedocs(
    modules  = [ImageFiltering, Kernel, KernelFactors, ImageFiltering.MapWindow],
    format   = format,
    sitename = "ImageFiltering",
    pages    = [
        "index.md",
        demos,
        "Function reference" => "function_reference.md"
    ]
)

demos_cb()

deploydocs(
    repo   = "github.com/JuliaImages/ImageFiltering.jl.git",
    push_preview = true
)


text:
using Parameters

@with_kw mutable struct Constants
    CellTyp::String = "LG M50" #"Doyle_94"
    T::Float64 = 298.15 # Cell Temperature (K)
    T_ref::Float64 = 298.15 # Reference Temperature (K)
    t_plus::Float64 = 0.2594 # Inital Transference Number
    tpf::Function = ce -> -0.1287*(ce/1000)^3+0.4106*(ce/1000)^2-0.4717*(ce/1000)+0.4492 # Transference Number Function - Requires ce in dm-3
    De::Float64 = 1.769e-10   # Inital Electrolyte Diffusivity (J mol⁻¹)
    Def::Function = ce -> 8.794e-11*(ce/1000)^2-3.972e-10*(ce/1000)+4.862e-10 # Electrolyte Diffusivity Function - Requires ce in dm-3
    SOC::Float64 = 0.8 # Initial State of Charge
    ce0::Float64 = 1000 # Initial Electrolyte Concentration (mol m⁻³)
    Ea_κ::Float64 = 0.
    Ea_De::Float64 = 0
    CC_A::Float64 =  0.1027  #Electrode Plate Area
    κ::Float64 = 0.9487
    κf::Function = ce -> 0.1297*(ce/1000)^3-2.51*(ce/1000)^1.5+3.329*(ce/1000) #Requires ce in dm-3
    Uocp::Function = (Electrode, θ) ->
        if Electrode == "Neg"
            Uocp = @. 1.97938*2.7182818284*exp(-39.3631*θ) + 0.2482 - 0.0909*tanh(29.8538*(θ - 0.1234)) - 0.04478*tanh(14.9159*(θ - 0.2769)) - 0.0205*tanh(30.4444*(θ - 0.6103))
        else
            Uocp = @. -0.8090*θ + 4.4875 - 0.0428*tanh(18.5138*(θ - 0.5542)) - 17.7326*tanh(15.7890*(θ - 0.3117)) + 17.5842*tanh(15.9308*(θ - 0.3120))
        end
    ∂Uocp::Function = (Electrode,θ) ->
        if Electrode == "Neg"
            #∂Uocp = @. 0.667934002(tanh(14.9159*θ - 4.13021271)^2) + 2.71371042(tanh(29.8538*θ - 3.68395892)^2) + 0.6241102(tanh(30.4444*θ - 18.58021732)^2) - 4.005754622 - (77.91453287630758(2.7182818284^(-39.3631*θ)))
            ∂Uocp = @. -0.62411*((sech(18.5802 - 30.4444*θ))^2 + 4.34813*(sech(3.68396 - 29.8538*θ))^2 + 1.07022*(sech(4.13021 - 14.9159*θ))^2) - 211.794*exp(-39.3631*θ)
        else
            ∂Uocp = @. 279.9800214*(tanh(15.789*θ - 4.9214313)^2) + 0.79239064*(tanh(18.5138*θ - 10.26034796)^2) - 1.4510386800000594 - (280.13037336*(tanh(15.9308*θ - 4.9704096)^2))
        end
    Ce_M::Int64 = 1     #(Rewritten)
    D1::Float64 = 1.0
    D2::Float64 = 1.0
    D3::Float64 = 1.0
    Ltot::Float64 = 0.
    Lnegsep::Float64 = 0.
    CeRootRange::Float64 = 10
end

@with_kw mutable struct Negative
    L::Float64 = 85.2e-6    # Negative Electrode Length
    Rs::Float64 = 5.86e-6   # Particle Radius [m]
    Ds::Float64 = 3.3e-14   # Solid Diffusivity [m^2/s]
    Ea_σ::Float64 = 0     # Activation Energy Solid Conductivity
    Ea_Ds::Float64 = 0    # Activation Energy Solid Diffusivity
    ϵ_s::Float64 = 0.75     # Active Material Volume Fraction
    ϵ_e::Float64 = 0.25     # Porosity of Negative Electrode
    De_brug::Float64 = 1.5  # Bruggeman Diffusion Exponent
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Exponent
    σ::Float64 = 215    # Solid Phase Conductivity
    σ_brug::Float64 = 1.5   # Bruggeman Solid Conductivity Exponent
    θ_100::Float64 = 0.910612 #0.9014 # Theta @ 100% Lithium Concentration
    θ_0::Float64 = 0.0263473 #0.0279 # Theta @ 0% Lithium Concentration
    cs_max::Float64 = 33133 # Max Electrode Concentration
    α::Float64 = 0.5    # Alpha Factor
    k_norm::Float64 = 6.716047e-12 #6.48e-7 #2.12e-10#1E-5 #7.226781e-7 # Reaction Rate
    Ea_κ::Float64 = 35000   # Activation Energy
    RFilm::Float64 = 0 # Film Resistance - Ωm²
    D1::Float64 = 1.   # Init Value
    D1f::Function = De -> De * ϵ_e^De_brug # Effective Diffusivity
    as::Float64 = 3.0*ϵ_s/Rs # Specific Interfacial Surf. Area
    β::Array{ComplexF64} = [0 - 0.0im] # Init Value
    β!::Function = (s) -> @. Rs*sqrt(s/Ds) # Negative β
end

@with_kw mutable struct Positive
    L::Float64 = 75.6e-6    # Positive Electrode Length
    Rs::Float64 = 5.22e-6    # Particle radius [m]
    Ds::Float64 = 4e-15   # Solid diffusivity [m^2/s]
    Ea_σ::Float64 = 0
    Ea_Ds::Float64 = 0
    ϵ_s::Float64 = 0.665    # Active Material Volume Fraction
    ϵ_e::Float64 = 0.335    # Porosity of positive electrode
    De_brug::Float64 = 1.5  # Bruggeman Diffusivity Exponent
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Exponent
    σ::Float64 = 0.18 #0.847  # Solid Phase Conductivity
    σ_brug::Float64 = 1.5   # Bruggeman Solid Conductivity Exponent
    θ_100::Float64 = 0.263849 #0.27 # Theta @ 100% Lithium Concentration
    θ_0::Float64 = 0.853974 #0.9084 # Theta @ 0% Lithium Concentration
    cs_max::Float64 = 63104 # Max Electrode Concentration
    α::Float64 = 0.5    # Alpha Factor
    k_norm::Float64 = 3.54458e-11 #3.42e-6 #1.12e-9#4E-05 #7.264272e-6 # Reaction Rate
    Ea_κ::Float64 = 17800   # Activation Energy
    RFilm::Float64 = 0 # Film Resistance - Ωm²
    D3::Float64 = 1   # Init Value
    D3f::Function = De -> De * ϵ_e^De_brug
    as::Float64 = 3.0*ϵ_s/Rs # Specific interfacial surf. area
    β::Array{ComplexF64} = [0 - 0.0im] # Init Value
    β!::Function = (s) -> @. Rs*sqrt(s/Ds) # Positive β
end

@with_kw mutable struct Seperator
    L::Float64 = 12e-6  # Seperator Length
    ϵ_e::Float64 = 0.47    # Porosity of separator
    De_brug::Float64 = 1.5  # Bruggeman Diffusivity Factor
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Factor
    D2::Float64 = 1   # Init Value
    D2f::Function = De -> De * ϵ_e^De_brug
end

@with_kw mutable struct Realisation
    Fs::Float64 = 4    # Sampling Frequency of Transfer Functions [Hz]
    SamplingT::Float64 = 0.25     # Final Model Sampling Time [s]
    M::Int64 = 4    # Model Order
    N::Int64 = 1    # Number of Inputs
    Tlen::Int64 = 16200 #Transfer Function Response Length [s] (Change to min)
    H1::Array{Int64,1} = 1:2500 #4000 #4612     # Hankel Dimensions 1
    H2::Array{Int64,1} = 1:2500 #4000 #4612     # Hankel Dimensions 2
    Outs::Int64 = 1    # Number of Outputs (Rewritten)
    Nfft::Int64 = 0
    f::UnitRange{Int64} = 0:1
    s::Array{ComplexF64} = [0 - 0.0im]
    Nfft!::Function = (Fs,Tlen) -> ceil(2^(log2(Fs*Tlen))) #2^(ceil(log2(Fs*Tlen))) Old way (constains to 2^ values)
    f!::Function = (Nfft) -> 0:Nfft-1
    s!::Function = (Fs,Nfft,f) -> transpose(((2im.*Fs)*tan.(pi.*f./Nfft)))
end

@with_kw mutable struct TransferFun
    #Electrolyte - Solid

    #6 - 4
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} = [Float64[0,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]

    #4 - 4
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} =  [Float64[0.0, 8.52e-5,9.72e-5,0.0001728], Float64[8.52e-5,9.72e-5,0.0001728], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]

    #8 - 4
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} = [Float64[0,2.84e-5,5.68e-5,8.52e-5,9.72e-5,1.224e-4,1.476e-4,1.728e-4], Float64[2.84e-5,5.68e-5,8.52e-5,9.72e-5,1.224e-4,1.476e-4,1.728e-4], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]

    # 9 - 6
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} = [Float64[0,2.84e-5,5.68e-5,8.52e-5,9.12e-5,9.72e-5,0.0001224,0.0001476,0.0001728], Float64[2.84e-5,5.68e-5,8.52e-5,9.12e-5,9.72e-5,0.0001224,0.0001476,0.0001728], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]

    # 4 - 2
    tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    Locs::Vector{Vector{Float64}} = [Float64[0, 8.52E-05, 9.72E-05, 1.728E-04], Float64[8.52E-05, 9.72E-05, 1.728E-04], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1]]

    # 4 - 6
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.2,0.4,0.6,0.8,1], Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1]]]

    # 6 - 2
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} = [Float64[0,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1]]

    # 6 - 6
    #tfs::Vector{Function} = [C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se]
    #Elec::Vector{String} = ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"]
    #Locs::Vector{Vector{Float64}} = [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.2,0.4,0.6,0.8,1], Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1]]

end

@with_kw mutable struct Params
    Const::Constants
    Neg::Negative
    Pos::Positive
    Sep::Seperator
    RA::Realisation
    Transfer::TransferFun
end

Cell = Params(Constants(),Negative(),Positive(),Seperator(),Realisation(),TransferFun())
Cell.Const.Lnegsep, Cell.Const.Ltot = Cell.Neg.L+Cell.Sep.L,Cell.Neg.L+Cell.Sep.L+Cell.Pos.L
Cell.Const.D1 = Cell.Const.De*Cell.Neg.ϵ_e^Cell.Neg.De_brug
Cell.Const.D2 = Cell.Const.De*Cell.Sep.ϵ_e^Cell.Sep.De_brug
Cell.Const.D3 = Cell.Const.De*Cell.Pos.ϵ_e^Cell.Pos.De_brug
Cell.Const.Ce_M = size(Cell.Transfer.Locs[1],1)
Cell.RA.Outs = sum([size(Cell.Transfer.Locs[i],1) for i in 1:length(Cell.Transfer.tfs)])

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
# Types for collections of arrays with restricted choice of keys and value types

module Multicollections

using ..Dorothy.Utils
using ..Dorothy.Graphs

export
		AbstractMulticollection, Multicollection, MulticollectionView,
		MulticollectionItem, collval, itemtocollprop, colltoitemprop,
		collitemname, treepaths, H2Hierarchy, H2Iterator, flattenh1, flattenh2,
		H3Hierarchy, H3IteratorH2, H3IteratorH1

abstract type AbstractMulticollection end

abstract type Multicollection <: AbstractMulticollection
	#=
	mutable struct
		header<:Header
		n::Int
		D<:Dict{Symbol,Any}
	end
	=#
end

struct MulticollectionView{T1<:Multicollection,
		T2<:AbstractVector{<:Integer}}  <: AbstractMulticollection
	C::T1
	I::T2

	@inline function MulticollectionView{T1,T2}(C::T1, I::T2) where
			{T1<:Multicollection,T2<:AbstractVector{<:Integer}}
		@boundscheck checkindexseries(C, I)
		new(C, I)
	end
end

struct MulticollectionItem{T<:Multicollection}
	C::T
	i::Int

	@inline function MulticollectionItem{T}(C::T, i::Integer) where
			{T<:Multicollection}
		@boundscheck checkbounds(C, i)
		new(C, i)
	end
end

@inline MulticollectionView(C::T1, I::T2) where
		{T1<:Multicollection,T2<:AbstractVector{<:Integer}} =
		MulticollectionView{T1,T2}(C, I)

@inline function MulticollectionView(C::MulticollectionView,
		I::AbstractVector{<:Integer})
	@boundscheck checkindexseries(C, I)
	@inbounds MulticollectionView(C.C, C.I[I])
end

MulticollectionItem(C::T, i::Integer) where {T<:Multicollection} =
		MulticollectionItem{T}(C, i)

MulticollectionItem(C::MulticollectionView{T}, i::Integer) where
		{T<:Multicollection} = MulticollectionItem{T}(C.C, C.I[i])

@inline Base.getproperty(C::AbstractMulticollection, name::Symbol) =
		getproperty(C, Val(name))

@inline Base.getproperty(C::AbstractMulticollection, ::Val{name}) where {name} =
		getindex(C, name)

@inline Base.getproperty(C::Multicollection, ::Val{:header}) =
		getfield(C, :header)

@inline Base.getproperty(C::Multicollection, ::Val{:n}) = getfield(C, :n)

@inline Base.getproperty(C::Multicollection, ::Val{:D}) = getfield(C, :D)

@inline Base.getproperty(C::MulticollectionView, ::Val{:C}) = getfield(C, :C)

@inline Base.getproperty(C::MulticollectionView, ::Val{:I}) = getfield(C, :I)

@inline Base.getproperty(C::MulticollectionView, ::Val{:header}) = C.C.header

@inline Base.getproperty(item::MulticollectionItem, name::Symbol) =
		Base.getproperty(item, Val(name))

@inline Base.getproperty(item::MulticollectionItem, ::Val{name}) where {name} =
		getindex(item, name)

@inline Base.getproperty(item::MulticollectionItem, ::Val{:C}) =
		getfield(item, :C)

@inline Base.getproperty(item::MulticollectionItem, ::Val{:i}) =
		getfield(item, :i)

@inline Base.getproperty(item::MulticollectionItem, ::Val{:header}) =
		item.C.header

@inline Base.setproperty!(C::AbstractMulticollection, name::Symbol, v) =
		setproperty!(C, Val(name), v)

@inline Base.setproperty!(C::AbstractMulticollection, ::Val{name}, v) where
		{name} = setindex!(C, v, name)

@inline Base.setproperty!(C::Multicollection, ::Val{:n}, v) =
		setfield!(C, :n, v)

@inline Base.setproperty!(item::MulticollectionItem, name::Symbol, v) =
		setindex!(item, v, name)

Base.propertynames(C::Multicollection, private::Bool = false) =
		private ? (:header, :n, :D, keys(C)...) : (:header, keys(C)...)

Base.propertynames(C::MulticollectionView, private::Bool = false) =
		private ? (:header, :C, :I, keys(C.C)...) : (:header, keys(C.C)...)

function Base.propertynames(item::MulticollectionItem, private::Bool = false)
	dynprops = [colltoitemprop(item, key) for key in keys(item.C)]
	private ? (:header, :C, :i, dynprops...) : (:header, dynprops...)
end

Base.show(io::IO, C::AbstractMulticollection) =
		print(io, "$(typeof(C))($(length(C)))")

function Base.show(io::IO, ::MIME"text/plain", C::Multicollection)
	n = length(C)
	nkeys = length(keys(C))
	itemname = collitemname(C)
	print(io, "$(n)-$(itemname) $(nkeys)-key $(typeof(C))")
	if nkeys > 0
		print(io, ":")
		for key in sort!(collect(keys(C)), by = lowercase∘string)
			print(io, "\n $(key)")
		end
	end
end

function Base.show(io::IO, ::MIME"text/plain", C::MulticollectionView)
	n = length(C)
	nkeys = length(keys(C))
	itemname = collitemname(C.C)
	print(io, "$(n)-$(itemname) $(nkeys)-key $(typeof(C.C)) view")
	if nkeys > 0
		print(io, ":")
		for key in sort!(collect(keys(C)), by = lowercase∘string)
			print(io, "\n $(key)")
		end
	end
end

Base.show(io::IO, item::MulticollectionItem) =
		print(io, "$(typeof(item.C))[$(item.i)]")

function Base.show(io::IO, ::MIME"text/plain", item::MulticollectionItem)
	n = length(item.C)
	itemkeys = keys(item)
	nkeys = length(itemkeys)
	itemname = collitemname(item.C)
	print(io, "$(nkeys)-key $(typeof(item.C)) $(itemname)")
	if nkeys > 0
		print(io, ":")
		for key in sort!(itemkeys, by = lowercase∘string)
			print(io, "\n $(key)")
		end
	end
end

Base.length(C::Multicollection) = C.n

Base.length(C::MulticollectionView) = length(C.I)

Base.firstindex(C::AbstractMulticollection) = 1

Base.lastindex(C::AbstractMulticollection) = length(C)

Base.eachindex(C::AbstractMulticollection) = Base.OneTo(length(C))

Base.checkbounds(::Type{Bool}, C::AbstractMulticollection, I) =
		checkbounds(Bool, 1:length(C), I)

Base.checkbounds(C::AbstractMulticollection, I) = checkbounds(1:length(C), I)

function Base.iterate(C::AbstractMulticollection, i::Int = 1)
	if i > length(C)
		nothing
	else
		MulticollectionItem(C, i), i + 1
	end
end

Base.eltype(::Type{T}) where {T<:AbstractMulticollection} =
		MulticollectionItem{T}

Base.similar(C::AbstractMulticollection, t::Type, n::Integer = length(C)) =
		similar(eachindex(C), t, n)

Base.parent(C::Multicollection) = C

Base.parent(C::MulticollectionView) = C.C

Base.parent(item::MulticollectionItem) = item.C

Base.parentindices(C::Multicollection) = eachindex(C)

Base.parentindices(C::MulticollectionView) = C.I

Base.parentindices(item::MulticollectionItem) = [item.i]

Base.view(C::Multicollection, I::AbstractVector{<:Integer}) =
		MulticollectionView(C, copy(I))

Base.view(C::MulticollectionView, I::AbstractVector{<:Integer}) =
		MulticollectionView(C, I)

Base.view(C::AbstractMulticollection, i::Integer) = MulticollectionView(C, [i])

@inline function Base.getindex(C::AbstractMulticollection)
	@boundscheck length(C) != 1 && throw(BoundsError(C, []))
	@inbounds C[1]
end

@inline Base.getindex(C::AbstractMulticollection, i::Integer) =
		MulticollectionItem(C, i)

@inline Base.getindex(C::Multicollection, I::AbstractVector{<:Integer}) =
		view(C, I)

@inline Base.getindex(C::Multicollection, key::Symbol) = C.D[key]

@inline Base.getindex(C::MulticollectionView, key::Symbol) = view(C.C[key], C.I)

@inline Base.getindex(item::MulticollectionItem, key::Symbol) =
		item.C[itemtocollprop(item, key)][item.i]

@inline function Base.setindex!(C::Multicollection, v, key::Symbol)
	if haskey(C, key)
		C[key] .= v
	else
		get!(C, key, v)
	end
end

@inline Base.setindex!(C::MulticollectionView, v, key::Symbol) = (C[key] .= v)

@inline Base.setindex!(item::MulticollectionItem, v, key::Symbol) =
		item.C[itemtocollprop(item, key)][item.i] = v

Base.haskey(C::Multicollection, key::Symbol) = haskey(C.D, key)

Base.haskey(C::MulticollectionView, key::Symbol) = haskey(C.C, key)

Base.haskey(item::MulticollectionItem, key::Symbol) =
		haskey(item.C, itemtocollprop(item, key))

Base.getkey(C::Multicollection, key::Symbol, default) =
		getkey(C.D, key, default)

Base.getkey(C::MulticollectionView, key::Symbol, default) =
		getkey(C.C, key, default)

Base.getkey(item::MulticollectionItem, key::Symbol, default) =
		haskey(item, key) ? key : default

Base.get(f, C::AbstractMulticollection, key::Symbol) =
		haskey(C, key) ? getindex(C, key) : f()

Base.get(C::AbstractMulticollection, key::Symbol, default) =
		(get(C, key) do; default; end)

Base.get(f, item::MulticollectionItem, key::Symbol) =
		haskey(item, key) ? getindex(item, key) : f()

Base.get(item::MulticollectionItem, key::Symbol, default) =
		(get(item, key) do; default; end)

function Base.get!(f, C::Multicollection, key::Symbol)
	get!(C.D, key) do
		collval(C, Val(key), length(C), f())
	end
end

Base.get!(C::Multicollection, key::Symbol, default) =
		(get!(C, key) do; default; end)

function Base.delete!(C::Multicollection, key::Symbol)
	delete!(C.D, key)
	C
end

Base.keys(C::Multicollection) = keys(C.D)

Base.keys(C::MulticollectionView) = keys(C.C)

Base.keys(item::MulticollectionItem) =
		[colltoitemprop(item, key) for key in keys(item.C)]

Base.values(C::AbstractMulticollection) = [C[key] for key in keys(C)]

Base.pairs(C::AbstractMulticollection) = [(key, C[key]) for key in keys(C)]

Base.:(==)(C1::Multicollection, C2::Multicollection) = (C1 === C2)

Base.:(==)(C1::MulticollectionView, C2::MulticollectionView) =
		(C1.C == C2.C && C1.I == C2.I)

Base.:(==)(item1::MulticollectionItem, item2::MulticollectionItem) =
		(item1.C == item2.C && item1.i == item2.i)

function Base.resize!(C::Multicollection, n::Integer)
	@boundscheck n >= 0 || error("expected positive length")
	for val in values(C.D)
		resize!(val, n)
	end
	C.n = n
	C
end

function Base.empty!(C::Multicollection)
	for val in values(C.D)
		empty!(val)
	end
	C.n = 0
	C
end

Base.deleteat!(C::Multicollection, i::Integer) = deleteat!(C, [i])

@inline function Base.deleteat!(C::Multicollection, I::AbstractArray{<:Integer})
	@boundscheck checkindexseries(C, I)
	for val in values(C.D)
		deleteat!(val, I)
	end
	C.n -= length(I)
	C
end

Base.splice!(C::Multicollection, i::Integer) = deleteat!(C, i)

Base.splice!(C::Multicollection, i::Integer, replacement::MulticollectionItem) =
		splice!(C, i:i, view(parent(replacement), parentindices(replacement)))

Base.splice!(C::Multicollection, i::Integer,
		replacement::AbstractMulticollection) = splice!(C, i:i, replacement)

Base.splice!(C::Multicollection, range::AbstractUnitRange{<:Integer}) =
		deleteat!(C, range)

Base.splice!(C::Multicollection, range::AbstractUnitRange{<:Integer},
		replacement::MulticollectionItem) =
		splice!(C, range, view(parent(replacement), parentindices(replacement)))

function Base.splice!(C::Multicollection, range::AbstractUnitRange{<:Integer},
		replacement::AbstractMulticollection)
	@boundscheck begin
		checkbounds(C, range)
		parent(replacement) !== C ||
				error("replacement must belong to a different $(typeof(C))")
		for key in keys(C)
			if ! haskey(replacement, key)
				error("missing $(key) key in replacement")
			end
		end
	end
	nreplaced = length(range)
	nreplacement = length(replacement)
	spliced = similar(C, nreplaced)
	if nreplacement == 0
		@inbounds deleteat!(C, range)
	else
		replaced = view(C, range)
		for (key, val) in pairs(C.D)
			get!(spliced, key, replaced[key])
			splice!(val, range, replacement[key])
		end
	end
	C.n += nreplacement - nreplaced
	spliced
end

Base.append!(C::Multicollection, src::MulticollectionItem) =
		append!(C, view(parent(src), parentindices(src)))

function Base.append!(C::Multicollection, src::AbstractMulticollection)
	n = length(C)
	splice!(C, n+1:n, src)
	C
end

Base.prepend!(C::Multicollection, src::MulticollectionItem) =
		prepend!(C, view(parent(src), parentindices(src)))

function Base.prepend!(C::Multicollection, src::AbstractMulticollection)
	splice!(C, 1:0, src)
	C
end

Base.insert!(C::Multicollection, i::Integer, src::MulticollectionItem) =
		insert!(C, i, view(parent(src), parentindices(src)))

function Base.insert!(C::Multicollection, i::Integer,
		src::AbstractMulticollection)
	splice!(C, i:i-1, src)
	C
end

function collval end

@inline itemtocollprop(item::MulticollectionItem, name::Symbol) =
		itemtocollprop(item, Val(name))

@inline itemtocollprop(::MulticollectionItem, ::Val{name}) where name = name

@inline colltoitemprop(item::MulticollectionItem, name::Symbol) =
		colltoitemprop(item, Val(name))

collitemname(::Multicollection) = "item"

@inline colltoitemprop(::MulticollectionItem, ::Val{name}) where name = name

function treepaths(II::Vector{<:Vector{Int}}, n::Integer)
	J = Vector{Int}(undef, n)
	for (j, I) in enumerate(II)
		for i in I
			J[i] = j
		end
	end
	J
end

struct H2Hierarchy
	tree::Vector{Vector{Int}}
	paths::Vector{Int}
end

H2Hierarchy(tree::Vector{Vector{Int}}, n::Integer) =
		H2Hierarchy(tree, treepaths(tree, n))

function Base.iterate(h2::H2Hierarchy, i::Integer = 1)
	if i > 2
		nothing
	else
		h2[i], i + 1
	end
end

Base.length(::H2Hierarchy) = 2

function Base.getindex(h2::H2Hierarchy, i::Integer)
	if i == 1
		h2.tree
	elseif i == 2
		h2.paths
	else
		throw(BoundsError(h2, i))
	end
end

Base.firstindex(::H2Hierarchy) = 1

Base.lastindex(h2::H2Hierarchy) = length(h2)

struct H2Iterator{T1<:AbstractMulticollection,T2<:Multicollection}
	C::T1
	h2::H2Hierarchy

	H2Iterator{T1,T2}(C::T1, h2::H2Hierarchy) where
			{T1<:AbstractMulticollection,T2<:Multicollection} = new(C, h2)
end

H2Iterator(C::T, h2::H2Hierarchy) where {T<:AbstractMulticollection} =
		H2Iterator{T,typeof(parent(C))}(C, h2)

@inline function Base.iterate(iter::H2Iterator, i::Int = 1)
	if i > length(iter)
		nothing
	else
		@inbounds iter[i], i + 1
	end
end

Base.eltype(::Type{<:H2Iterator{T1,T2}}) where
		{T1<:AbstractMulticollection,T2<:Multicollection} =
		MulticollectionView{T2,Vector{Int}}

Base.length(iter::H2Iterator) = length(iter.h2.tree)

@inline function Base.getindex(iter::H2Iterator{T1,T2}, i::Integer) where
		{T1<:AbstractMulticollection,T2<:Multicollection}
	@boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
	view(iter.C, iter.h2.tree[i])
end

Base.firstindex(::H2Iterator) = 1

Base.lastindex(iter::H2Iterator) = length(iter)

function flattenh1(III::Vector{<:Vector{Vector{Int}}})
	JJ = Vector{Int}[]
	for II in III
		J = Int[]
		for I in II
			for i in I
				push!(J, i)
			end
		end
		push!(JJ, J)
	end
	JJ
end

function flattenh2(III::Vector{<:Vector{Vector{Int}}})
	JJ = Vector{Int}[]
	for II in III
		for I in II
			push!(JJ, I)
		end
	end
	JJ
end

struct H3Hierarchy
	tree::Vector{Vector{Vector{Int}}}
	flath1::H2Hierarchy
	flath2::H2Hierarchy
end

function H3Hierarchy(tree::Vector{Vector{Vector{Int}}}, n::Integer)
	flath1tree = flattenh1(tree)
	flath2tree = flattenh2(tree)
	flath1 = H2Hierarchy(flath1tree, n)
	flath2 = H2Hierarchy(flath2tree, n)
	H3Hierarchy(tree, flath1, flath2)
end

function Base.iterate(h3::H3Hierarchy, i::Integer = 1)
	if i > 3
		nothing
	else
		h3[i], i + 1
	end
end

Base.IteratorEltype(::Type{<:H3Hierarchy}) = Base.EltypeUnknown()

Base.length(::H3Hierarchy) = 3

function Base.getindex(h3::H3Hierarchy, i::Integer)
	if i == 1
		h3.tree
	elseif i == 2
		h3.flath1
	elseif i == 3
		h3.flath2
	else
		throw(BoundsError(h3, i))
	end
end

Base.firstindex(::H3Hierarchy) = 1

Base.lastindex(h3::H3Hierarchy) = length(h3)

struct H3IteratorH2{T1<:AbstractMulticollection,T2<:Multicollection}
	C::T1
	h3::H3Hierarchy

	H3IteratorH2{T1,T2}(C::T1, h3::H3Hierarchy) where
			{T1<:AbstractMulticollection,T2<:Multicollection} = new(C, h3)
end

H3IteratorH2(C::T, h3::H3Hierarchy) where {T<:AbstractMulticollection} =
		H3IteratorH2{T,typeof(parent(C))}(C, h3)

@inline function Base.iterate(iter::H3IteratorH2, i::Int = 1)
	if i > length(iter)
		nothing
	else
		@inbounds iter[i], i + 1
	end
end

Base.eltype(::Type{<:H3IteratorH2{T1,T2}}) where
		{T1<:AbstractMulticollection,T2<:Multicollection} =
		MulticollectionView{T2,Vector{Int}}

Base.length(iter::H3IteratorH2) = length(iter.h3.tree)

@inline function Base.getindex(iter::H3IteratorH2{T1,T2}, i::Integer) where
		{T1<:AbstractMulticollection,T2<:Multicollection}
	@boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
	view(iter.C, iter.h3.flath1.tree[i])
end

Base.firstindex(::H3IteratorH2) = 1

Base.lastindex(iter::H3IteratorH2) = length(iter)

struct H3IteratorH1{T1<:AbstractMulticollection,T2<:Multicollection}
	C::T1
	i::Int
	h3::H3Hierarchy

	function H3IteratorH1{T1,T2}(C::T1, i::Integer, h3::H3Hierarchy) where
			{T1<:AbstractMulticollection,T2<:Multicollection}
		@boundscheck checkbounds(h3.tree, i)
		new(C, i, h3)
	end
end

H3IteratorH1(C::T, i::Integer, h3::H3Hierarchy) where
		{T<:AbstractMulticollection} =
		H3IteratorH1{T,typeof(parent(C))}(C, i, h3)

@inline function Base.iterate(iter::H3IteratorH1, i::Int = 1)
	if i > length(iter)
		nothing
	else
		@inbounds iter[i], i + 1
	end
end

Base.eltype(::Type{<:H3IteratorH1{T1,T2}}) where
		{T1<:AbstractMulticollection,T2<:Multicollection} =
		MulticollectionView{T2,Vector{Int}}

Base.length(iter::H3IteratorH1) = length(iter.h3.tree[iter.i])

@inline function Base.getindex(iter::H3IteratorH1{T1,T2}, i::Integer) where
		{T1<:AbstractMulticollection,T2<:Multicollection}
	@boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
	view(iter.C, iter.h3.tree[iter.i][i])
end

Base.firstindex(::H3IteratorH1) = 1

Base.lastindex(iter::H3IteratorH1) = length(iter)

end # module


text:
using Bifurcations
using LinearAlgebra
using OrdinaryDiffEq
using Setfield
using Parameters: @unpack

HAS_CUDA = try
    using CuArrays
    import CUDAnative
    true
catch err
    @error "CuArrays cannot be imported" exception = (err, catch_backtrace())
    global const CuArray = Nothing
    false
end

function phase_dynamics!(du, u, p, t)
    @unpack J, h, g = p
    tanh = u isa CuArray ? CUDAnative.tanh : Base.tanh
    du .= .- u .+ g .* (J * tanh.(u)) .+ h
    return du
end

function jacobian!(M, u, p, t)
    @unpack J, h, g = p
    cosh = u isa CuArray ? CUDAnative.cosh : Base.cosh
    M .= -I + g .* J .* Diagonal(inv.(cosh.(u)) .^ 2)
    return M
end

N = 10000
ode = ODEProblem(
    ODEFunction(phase_dynamics!; jac = jacobian!),
    randn(N),
    (0.0, 50.0),
    (
        J = randn(N, N) ./ .√(N),
        h = randn(N),
        g = 0.5,
    ),
)

ts = solve(ode, Tsit5())
ode0 = remake(ode, u0 = ts.u[end])

param_axis = @lens _.g
prob = BifurcationProblem(ode0, param_axis, (0.5, 1.5))
# prob = BifurcationProblem(ode, param_axis, (0.5, 1.5))

using Bifurcations.Continuations: ContinuationSolver, ContinuationOptions
opts = ContinuationOptions(
    atol = 0.05,
    rtol = 0.05,
    max_branches = 0,
    max_samples = 1,  # comment this out to do actual bifurcation analysis
    nominal_distance = 0.01N,
    # bidirectional_first_sweep = false,
    # start_from_nearest_root = true,
    verbose = true,
)
solver = ContinuationSolver(prob, opts)
# solve(prob; start_from_nearest_root = true)
# sol = @time solve!(solver)
# sol = @time solve!(solver)

if HAS_CUDA
    ode_gpu = remake(
        ode0;
        u0 = cu(ode0.u0),
        p = (
            J = cu(ode0.p.J),
            h = cu(ode0.p.h),
            g = ode0.p.g,
        ),
    )
    prob_gpu = BifurcationProblem(ode_gpu, param_axis, (0.5, 1.5))
    solver_gpu = ContinuationSolver(prob_gpu, opts)
    # solve(prob_gpu; start_from_nearest_root = true)
    # sol_gpu = @time solve!(solver_gpu)
    # sol_gpu = @time solve!(solver_gpu)
end


text:
module TwoVarianceComponentTest

using VarianceComponentModels, MathProgBase, Ipopt, Random, Test, LinearAlgebra

Random.seed!(123)

# generate data from a d-variate response variane component model
n = 100   # no. observations
d = 2     # no. categories
m = 2     # no. variance components
p = 2     # no. covariates
X = randn(n, p)
B = ones(p, d)
Σ = ntuple(x -> zeros(d, d), m)
for i in 1:m
  Σi = randn(d, d)
  copyto!(Σ[i], Σi' * Σi)
end
## make the first variance component 0 matrix
#fill!(Σ[1], 0.0)
V = ntuple(x -> zeros(n, n), m)
for i = 1:m-1
  Vi = randn(n, 50)
  copyto!(V[i], Vi * Vi')
end
copyto!(V[m], Matrix(I, n, n))
# form Ω
Ω = zeros(n*d, n*d)
for i = 1:m
  global Ω += kron(Σ[i], V[i])
end
Ωchol = cholesky(Ω)
Y = X * B + reshape(Ωchol.L * randn(n*d), n, d)

@info "Forming VarianceComponentModel from data"
@inferred VarianceComponentVariate(Y, X, V)
vcdata  = VarianceComponentVariate(Y, X, V)
vcmodel = VarianceComponentModel(vcdata)

@info "Pre-compute eigen-decomposition and rotate data"
vcdatarot  = TwoVarCompVariateRotate(vcdata)
vcmodelrot = TwoVarCompModelRotate(vcmodel)

@info "Evaluate log-pdf"
#@code_warntype logpdf(vcmodelrot, vcdatarot)
@inferred logpdf(vcmodelrot, vcdatarot)
@test logpdf(vcmodel, vcdata) == logpdf(vcmodelrot, vcdatarot)
@test (logpdf(vcmodelrot, [vcdatarot vcdatarot; vcdatarot vcdatarot]) -
    logpdf(vcmodel, [vcdata vcdata; vcdata vcdata])) ≈ 0.0

@info "Evaluate gradient"
∇ = zeros(2d^2)
#@code_warntype gradient!(∇, vcmodelrot, vcdatarot)
@inferred VarianceComponentModels.gradient!(∇, vcmodelrot, vcdatarot)
@inferred VarianceComponentModels.gradient!(∇, vcmodel, vcdatarot)
@inferred VarianceComponentModels.gradient!(∇, vcmodel, vcdata)
@test norm(VarianceComponentModels.gradient(vcmodel, vcdata) -
    VarianceComponentModels.gradient(vcmodelrot, vcdatarot)) ≈ 0.0
@test norm(VarianceComponentModels.gradient(vcmodel, vcdata) -
    VarianceComponentModels.gradient(vcmodel, vcdatarot)) ≈ 0.0
@test norm(VarianceComponentModels.gradient(vcmodel, [vcdata vcdata]) -
  2.0VarianceComponentModels.gradient(vcmodel, vcdata)) ≈ 0.0
@test norm(VarianceComponentModels.gradient(vcmodel, [vcdata vcdata]) -
  VarianceComponentModels.gradient(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

@info "Evaluate Fisher information matrix of Σ"
H = zeros(2d^2, 2d^2)
#@code_warntype fisher!(H, vcmodelrot, vcdatarot)
@inferred fisher_Σ!(H, vcmodelrot, vcdatarot)
@inferred fisher_Σ!(H, vcmodel, vcdata)
@test norm(fisher_Σ(vcmodel, vcdata) - fisher_Σ(vcmodelrot, vcdatarot)) ≈ 0.0
@test norm(fisher_Σ(vcmodel, vcdata) - fisher_Σ(vcmodel, vcdatarot)) ≈ 0.0
@test norm(fisher_Σ(vcmodel, [vcdata vcdata]) -
    2fisher_Σ(vcmodel, vcdata)) ≈ 0.0
@test norm(fisher_Σ(vcmodel, [vcdata vcdata]) -
    fisher_Σ(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

@info "Evaluate Fisher information matrix of B"
H = zeros(p * d, p * d)
#@code_warntype fisher_B!(H, vcmodelrot, vcdatarot)
@inferred fisher_B!(H, vcmodelrot, vcdatarot)
@inferred fisher_B!(H, vcmodel, vcdatarot)
@inferred fisher_B!(H, vcmodel, vcdata)
@test norm(fisher_B(vcmodel, vcdata) - fisher_B(vcmodelrot, vcdatarot)) ≈ 0.0
@test norm(fisher_B(vcmodel, vcdata) - fisher_B(vcmodel, vcdatarot)) ≈ 0.0
@test norm(fisher_B(vcmodel, [vcdata vcdata]) -
    2.0fisher_B(vcmodel, vcdata)) ≈ 0.0
@test norm(fisher_B(vcmodel, [vcdata vcdata]) -
    fisher_B(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

@info "Find MLE using Fisher scoring"
vcmfs = deepcopy(vcmodel)
#@code_warntype mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
#@inferred mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
logl_fs, _, _, Σcov_fs, Bse_fs, = mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
logl_fs_array, = mle_fs!(vcmfs, [vcdatarot vcdatarot]; solver = :Ipopt)
@show vcmfs.B
@show Bse_fs
@show B

@info "Find MLE using MM algorithm"
vcmmm = deepcopy(vcmodel)
#@code_warntype mle_mm!(vcmm, vcdatarot)
#@inferred mle_mm!(vcmmm, vcdatarot)
logl_mm, _, _, _, Bse_mm, = mle_mm!(vcmmm, vcdatarot)
@test abs(logl_fs - logl_mm) / (abs(logl_fs) + 1.0) < 1.0e-4
@show vcmmm.B
@show Bse_mm
@show B
vcmmm = deepcopy(vcmodel)
logl_mm_array, = mle_mm!(vcmmm, [vcdatarot vcdatarot])
@test abs(logl_fs_array - logl_mm_array) / (abs(logl_fs_array) + 1.0) < 1.0e-4

@info "Find MLE using Fisher scoring (linear equality + box constraints)"
vcmfs       = deepcopy(vcmodel)
vcmfs.A     = [1.0 -1.0 zeros(1, p*d-2)]
vcmfs.sense = '='
vcmfs.b     = 0.0
vcmfs.lb    = 0.0
vcmfs.ub    = 1.0
logl_fs, _, _, Σcov_fs, Bse_fs, = mle_fs!(vcmfs, vcdatarot;
    solver = :Ipopt, qpsolver = :Ipopt)
@show vcmfs.B
@test vcmfs.B[1] ≈ vcmfs.B[2]
@test all(vcmfs.B .≥ 0.0)
@test all(vcmfs.B .≤ 1.0)

@info "Find MLE using MM algorithm (linear equality + box constraints)"
vcmm = deepcopy(vcmodel)
vcmm.A = [1.0 -1.0 zeros(1, p*d-2)]
vcmm.sense = '='
vcmm.b = 0.0
vcmm.lb = 0.0
vcmm.ub = 1.0
logl_mm, _, _, Σcov_mm = mle_mm!(vcmm, vcdatarot; qpsolver = :Ipopt)
@show vcmm.B
@test vcmm.B[1] ≈ vcmm.B[2]
@test all(vcmm.B .≥ 0.0)
@test all(vcmm.B .≤ 1.0)
@test abs(logl_fs - logl_mm) / (abs(logl_fs) + 1.0) < 1.0e-4

@info "Heritability estimation"
h, h_se = heritability(vcmfs.Σ, Σcov_fs)
@show h, h_se
@test all(0.0 .≤ h .≤ 1.0)

@info "test fit_mle (FS)"
vcmmle = deepcopy(vcmodel)
logl_mle, _, _, Σcov_mle, Bse_mle, = fit_mle!(vcmmle, vcdata; algo = :FS)
@show vcmmle.B, Bse_mle, B
vcmle = deepcopy(vcmodel)
fit_mle!(vcmmle, [vcdata vcdata]; algo = :FS)

@info "test fit_mle (MM)"
vcmmle = deepcopy(vcmodel)
logl_mle, _, _, Σcov_mle, Bse_mle, = fit_mle!(vcmmle, vcdata; algo = :MM)
@show vcmmle.B, Bse_mle, B

@info "test fit_reml (FS)"
vcmreml = deepcopy(vcmodel)
logl_reml, _, _, Σcov_reml, Bse_reml, = fit_reml!(vcmreml, vcdata; algo = :FS)
@show vcmreml.B, Bse_reml, B

@info "test fit_reml (MM)"
vcmreml = deepcopy(vcmodel)
logl_reml, _, _, Σcov_reml, Bse_reml, = fit_reml!(vcmreml, vcdata; algo = :MM)
@show vcmreml.B, Bse_reml, B

end # module VarianceComponentTypeTest

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
"""
    Boxplot <: Drawing2D

Boxplot(s).
"""
@with_kw mutable struct Boxplot <: Drawing2D
    stats::Matrix{Float64}  # quantile etc
    nobj ::Int
    #
    boxstyles::Vector{BoxplotStyle}
    #
    horiz::Bool = false # vertical boxplots by default
end
Boxplot(d, n) = Boxplot(stats=d, nobj=n, boxstyles=nvec(n, BoxplotStyle))

"""
    Heatmap <: Drawing2D

Heatmap of a matrix.
"""
@with_kw mutable struct Heatmap <: Drawing2D
    data::Matrix{Int}
    zmin::Float64
    zmax::Float64
    #
    cmap::Vector{Color} = colormap("RdBu", 10)
    cmiss::Color = c"white" # box filling for missing values
    # transpose
    # whether to write the matrix as a transpose
    # this is useful because GLE can deal only with 1000-cols
    # files at most (at least with the way we do the heatmap now)
    transpose::Bool = false
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_export_constraints(m, p; _n="")

    ##Constraint (8e): Production export no greater than production
    @constraint(m, [t in p.elec_techs, ts in p.time_steps_with_grid],
        p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts]
        >= sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t])
        + m[Symbol("dvCurtail"*_n)][t, ts]
    )

    ##Constraint (8f): Total sales to grid no greater than annual allocation - storage tiers
    @constraint(m,
        p.hours_per_timestep * sum( m[Symbol("dvProductionToGrid"*_n)][t, u, ts]
            for t in p.elec_techs, u in p.export_bins_by_tech[t], ts in p.time_steps_with_grid
        ) <= p.max_grid_export_kwh
    )

    ### Constraint set (9): Net Meter Module
    ##Constraint (9c): Net metering only -- can't sell more than you purchase
    # note that hours_per_timestep is cancelled on both sides, but used for unit consistency (convert power to energy)
    NEM_techs = String[t for t in p.elec_techs if :NEM in p.export_bins_by_tech[t]]
    if !isempty(NEM_techs)
        @constraint(m,
            p.hours_per_timestep * sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts]
            for t in NEM_techs, ts in p.time_steps)
            <= p.hours_per_timestep * sum( m[Symbol("dvGridPurchase"*_n)][ts] for ts in p.time_steps)
        )
    end
end

function add_monthly_peak_constraint(m, p; _n="")
	## Constraint (11d): Monthly peak demand is >= demand at each hour in the month
	@constraint(m, [mth in p.months, ts in p.etariff.time_steps_monthly[mth]],
    m[Symbol("dvPeakDemandMonth"*_n)][mth] >= m[Symbol("dvGridPurchase"*_n)][ts]
    )
end

function add_tou_peak_constraint(m, p; _n="")
    ## Constraint (12d): Ratchet peak demand is >= demand at each hour in the ratchet`
    @constraint(m, [r in p.ratchets, ts in p.etariff.tou_demand_ratchet_timesteps[r]],
        m[Symbol("dvPeakDemandTOU"*_n)][r] >= m[Symbol("dvGridPurchase"*_n)][ts]
    )
end

function add_mincharge_constraint(m, p; _n="")
    @constraint(m,
        m[Symbol("MinChargeAdder"*_n)] >= m[Symbol("TotalMinCharge"*_n)] - ( m[Symbol("TotalEnergyChargesUtil"*_n)] +
        m[Symbol("TotalDemandCharges"*_n)] + m[Symbol("TotalExportBenefit"*_n)] + m[Symbol("TotalFixedCharges"*_n)] )
    )
end

function add_simultaneous_export_import_constraint(m, p; _n="")
    @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
          m[Symbol("dvGridPurchase"*_n)][ts]
        + sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.storage.types)
        - (1 - m[Symbol("binNoGridPurchases"*_n)][ts]) * 1.0E9 <= 0
    )
    @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
        sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.elec_techs, u in p.export_bins_by_tech[t])
        - m[Symbol("binNoGridPurchases"*_n)][ts] * 1.0E9 <= 0
    )
end

function add_elec_utility_expressions(m, p; _n="")

    if !isempty(p.etariff.export_bins) & !isempty(p.techs)
        # NOTE: levelization_factor is baked into dvProductionToGrid
        m[Symbol("TotalExportBenefit"*_n)] = @expression(m, p.pwf_e * p.hours_per_timestep *
            sum( sum(p.etariff.export_rates[u][ts] * m[Symbol("dvProductionToGrid"*_n)][t, u, ts]
                 for u in p.etariff.export_bins, t in p.techs_by_exportbin[u])
            for ts in p.time_steps)
        )
    else
        m[Symbol("TotalExportBenefit"*_n)] = 0
    end

    m[Symbol("TotalEnergyChargesUtil"*_n)] = @expression(m, p.pwf_e * p.hours_per_timestep *
        sum( p.etariff.energy_rates[ts] * m[Symbol("dvGridPurchase"*_n)][ts] for ts in p.time_steps)
    )

    if !isempty(p.etariff.tou_demand_rates)
        m[Symbol("DemandTOUCharges"*_n)] = @expression(m,
            p.pwf_e * sum( p.etariff.tou_demand_rates[r] * m[Symbol("dvPeakDemandTOU"*_n)][r] for r in p.ratchets)
        )
    else
        m[Symbol("DemandTOUCharges"*_n)] = 0
    end

    if !isempty(p.etariff.monthly_demand_rates)
        m[Symbol("DemandFlatCharges"*_n)] = @expression(m,
            p.pwf_e * sum( p.etariff.monthly_demand_rates[mth] * m[Symbol("dvPeakDemandMonth"*_n)][mth] for mth in p.months)
        )
    else
        m[Symbol("DemandFlatCharges"*_n)] = 0
    end

    m[Symbol("TotalDemandCharges"*_n)] = m[Symbol("DemandTOUCharges"*_n)] + m[Symbol("DemandFlatCharges"*_n)]

    m[Symbol("TotalFixedCharges"*_n)] = p.pwf_e * p.etariff.fixed_monthly_charge * 12

    if p.etariff.annual_min_charge > 12 * p.etariff.min_monthly_charge
        m[Symbol("TotalMinCharge"*_n)] = p.etariff.annual_min_charge
    else
        m[Symbol("TotalMinCharge"*_n)] = 12 * p.etariff.min_monthly_charge
    end

	if m[Symbol("TotalMinCharge"*_n)] >= 1e-2
		add_mincharge_constraint(m, p)
	else
		@constraint(m, m[Symbol("MinChargeAdder"*_n)] == 0)
	end

    m[Symbol("TotalElecBill"*_n)] = (
        m[Symbol("TotalEnergyChargesUtil"*_n)]
        + m[Symbol("TotalDemandCharges"*_n)]
        + m[Symbol("TotalExportBenefit"*_n)]
        + m[Symbol("TotalFixedCharges"*_n)]
        + 0.999 * m[Symbol("MinChargeAdder"*_n)]
    )
    #= Note: 0.999 * MinChargeAdder in Objective b/c when
        TotalMinCharge > (TotalEnergyCharges + TotalDemandCharges + TotalExportBenefit + TotalFixedCharges)
		it is arbitrary where the min charge ends up (eg. could be in TotalDemandCharges or MinChargeAdder).
		0.001 * MinChargeAdder is added back into LCC when writing to results.
    =#
    nothing
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
"""

"""
struct System <: AbstractSystem
	f::Function
	X::Hyperrectangle
	state_dim::Int

	function System(f::Function, X::Hyperrectangle{T,Vector{T},Vector{T}}, xdim::Int) where T
		new(f, X, xdim)
	end
end

function System(f::Function, xlb::Vector{T}, xub::Vector{T}) where T
	X = Hyperrectangle(low=xlb, high=xub)
	return System(f, X, dim(X))
end

struct UncertainSystem <: AbstractUncertainSystem
	f::Function
	X::Hyperrectangle
	W::Hyperrectangle
	state_dim::Int
	uncertainty_dim::Int

	function UncertainSystem(f::Function, X::Hyperrectangle{T,Vector{T},Vector{T}}, W::Hyperrectangle{T,Vector{T},Vector{T}}, xdim::Int, wdim::Int) where T
		new(f, X, W, xdim, wdim)
	end
end

function UncertainSystem(f::Function, xlb::Vector{T}, xub::Vector{T}, wlb::Vector{T}, wub::Vector{T}) where T
	X = Hyperrectangle(low=xlb, high=xub)
	W = Hyperrectangle(low=wlb, high=wub)
	return UncertainSystem(f, X, W, dim(X), dim(W))
end

struct ControlSystem <: AbstractControlSystem
	f::Function
	X::Hyperrectangle
	U::Hyperrectangle
	state_dim::Int
	input_dim::Int

	function ControlSystem(f::Function, X::Hyperrectangle{T,Vector{T},Vector{T}}, U::Hyperrectangle{T,Vector{T},Vector{T}}, xdim::Int, udim::Int) where T
		new(f, X, U, xdim, udim)
	end
end

function ControlSystem(f::Function, xlb::Vector{T}, xub::Vector{T}, ulb::Vector{T}, uub::Vector{T}) where T
	X = Hyperrectangle(low=xlb, high=xub)
	U = Hyperrectangle(low=ulb, high=uub)
	return ControlSystem(f, X, U, dim(X), dim(U))
end

struct UncertainControlSystem <: AbstractUncertainControlSystem
	f::Function
	X::Hyperrectangle
	U::Hyperrectangle
	W::Hyperrectangle
	state_dim::Int
	input_dim::Int
	uncertainty_dim::Int

	function UncertainControlSystem(f::Function, X::Hyperrectangle{T,Vector{T},Vector{T}}, U::Hyperrectangle{T,Vector{T},Vector{T}},
		W::Hyperrectangle{T,Vector{T},Vector{T}}, xdim::Int, udim::Int,wdim::Int) where T
		new(f, X, U, W, xdim, udim, wdim)
	end
end

function UncertainControlSystem(f::Function, xlb::Vector{T}, xub::Vector{T}, ulb::Vector{T}, uub::Vector{T}, wlb::Vector{T}, wub::Vector{T}) where T
	X = Hyperrectangle(low=xlb, high=xub)
	U = Hyperrectangle(low=ulb, high=uub)
	W = Hyperrectangle(low=wlb, high=wub)
	return UncertainControlSystem(f, X, U, W, dim(X), dim(U), dim(W))
end

function system(f::Function, xlb::Vector{T}, xub::Vector{T}; ulb::Union{Vector{T},Nothing}=nothing, uub::Union{Vector{T},Nothing}=nothing, wlb::Union{Vector{T},Nothing}=nothing, wub::Union{Vector{T},Nothing}=nothing) where T
	if isnothing(ulb) && isnothing(uub) && isnothing(wlb) && isnothing(wub)
		return System(f, xlb, xub)
	elseif isnothing(wlb) && isnothing(wub)
		return ControlSystem(f, xlb, xub, ulb, uub)
	elseif isnothing(ulb) && isnothing(uub)
		return UncertainSystem(f, xlb, xub, wlb, wub)
	else
		return UncertainControlSystem(f, xlb, xub, ulb, uub, wlb, wub)
	end
end


text:
# Conways game of life
#
# includes an RLE decoder
# UnicodePlots is used to display the animation in the terminal

using UnicodePlots

function count_neighbors(m::Matrix{Bool}, i::CartesianIndex)
    c = 0
    for j in CartesianIndex.((-1:1)', -1:1)
        j == CartesianIndex(0,0) && continue

        if checkbounds(Bool, m, i+j)
            c += m[i+j]
        end
    end
    return c
end

function gameoflife_step(m::Matrix{Bool}; live, appear)
    nextm = similar(m)
    for i in CartesianIndices(m)
        if m[i]
            nextm[i] = count_neighbors(m, i) ∈ live
        else
            nextm[i] = count_neighbors(m, i) ∈ appear
        end
    end
    return nextm
end

function animate(m::Matrix{Bool}, dt, steps; live=2:3, appear=3)
    print(repr(spy(m; labels=false, title="", canvas=DotCanvas)))
    height = div(size(m, 1), 2, RoundUp) + 1
    t = time()
    for _ in 1:steps
        m = gameoflife_step(m; live=live, appear=appear)
        while time() - t < dt end
        t = time()
        print("\r\u1b[K\u1b[A"^height) # remove previous graph
        print(repr(spy(m; labels=false, title="", canvas=DotCanvas)))
    end
    println()
end

function RLE_decoder(RLE; x, y)
    m = zeros(Bool, y, x)
    i = CartesianIndex(1,1)
    n = 0
    for j in collect(RLE)
        nn = n == 0 ? 1 : n
        if isspace(j)
            continue
        elseif j == '!' # end
            break
        elseif j == 'b' # dead
            i += CartesianIndex(0, nn)
            n = 0
        elseif j == '$' # new line
            i = CartesianIndex(i[1]+nn, 1)
            n = 0
        elseif j ∈ '0':'9'
            n = 10*n + j - '0'
        else # alive
            for _ in 1:nn
                m[i] = true
                i += CartesianIndex(0, 1)
            end
            n = 0
        end
    end
    return m
end

m = RLE_decoder(raw"16bo30b$16bobo16bo11b$16b2o17bobo9b$obo10bo21b2o10b$b2o11b2o31b$bo11b
2o32b3$10b2o20b2o13b$11b2o19bobo9b3o$10bo21bo11bo2b$27bo17bob$27b2o18b$26bobo!",
x = 47,
y = 14)

# add margins
T = 0 # top
B = 10 # bottom
R = 0 # right
L = 0 # left

m = Bool[
    zeros(T, L)         zeros(T, size(m,2)) zeros(T, R)
    zeros(size(m,1), L) m                   zeros(size(m,1), R)
    zeros(B, L)         zeros(B, size(m,2)) zeros(B, R)
]

animate(m, 0.01, 200)


text:
_with_vdw(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in keys(vanderwaalsradius)

_with_cov(a::PDBAtom, resname_a::String) = a.element in keys(covalentradius)

ishydrophobic(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in _hydrophobic

"""
Returns true if the atom, e.g. `("HIS","CG")`, is an aromatic atom in the residue.
"""
isaromatic(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in _aromatic

"""
Returns true if the atom, e.g. `("ARG","NE")`, is a cationic atom in the residue.
"""
iscationic(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in _cationic

"""
Returns true if the atom, e.g. `("GLU","CD")`, is an anionic atom in the residue.
"""
isanionic(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in _anionic

"""
Returns true if the atom, e.g. `("ARG","N")`, is a donor in H bonds.
"""
ishbonddonor(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in keys(_hbond_donor)

"""
Returns true if the atom, e.g. `("ARG","O")`, is an acceptor in H bonds.
"""
ishbondacceptor(a::PDBAtom, resname_a::String) = (resname_a, a.atom) in keys(_hbond_acceptor)

"""
`any(f::Function, a::PDBResidue, b::PDBResidue, criteria::Function)`

Test if the function `f` is true for any pair of atoms between the residues `a` and `b`.
This function only test atoms that returns `true` for the fuction `criteria`.
"""
function Base.any(f::Function, a::PDBResidue, b::PDBResidue, criteria::Function)
    resname_a, resname_b = a.id.name, b.id.name
    a_atoms = a.atoms
    b_atoms = b.atoms
    indices_a = _find(x -> criteria(x, resname_a), a_atoms)
    indices_b = _find(x -> criteria(x, resname_b), b_atoms)
    if length(indices_a) != 0 && length(indices_b) != 0
        @inbounds for i in indices_a
            for j in indices_b
                if f(a_atoms[i], b_atoms[j], resname_a, resname_b)
                    return(true)
                end
            end
        end
    end
    return(false)
end

# Interaction types
# =================

# van der Waals
# -------------

"""
Test if two atoms or residues are in van der Waals contact using:
`distance(a,b) <= 0.5 + vanderwaalsradius[a] + vanderwaalsradius[b]`.
It returns distance `<= 0.5` if the atoms aren't in `vanderwaalsradius`.
"""
function vanderwaals(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    return( distance(a,b) <= 0.5 +
               get(vanderwaalsradius, (resname_a, a.atom), 0.0) +
               get(vanderwaalsradius, (resname_b, b.atom), 0.0) )
end

vanderwaals(a::PDBResidue, b::PDBResidue) = any(vanderwaals, a, b, _with_vdw)

# van der Waals clash
# -------------------

"""
Returns `true` if the distance between the atoms is less than the sum of the
`vanderwaalsradius` of the atoms. If the atoms aren't on the list (i.e. `OXT`), the
`vanderwaalsradius` of the element is used. If there is not data in the dict,
distance `0.0` is used.
"""
function vanderwaalsclash(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    return( distance(a,b) <= get(vanderwaalsradius, (resname_a, a.atom),
                                 get(vanderwaalsradius, (resname_a, a.element), 0.0)) +
                             get(vanderwaalsradius, (resname_b, b.atom),
                                 get(vanderwaalsradius, (resname_b, b.element), 0.0)) )
end

vanderwaalsclash(a::PDBResidue, b::PDBResidue) = any(vanderwaalsclash, a, b, _with_vdw)

# Covalent
# --------

"""
Returns `true` if the distance between atoms is less than the sum of the `covalentradius`
of each atom.
"""
function covalent(a::PDBAtom, b::PDBAtom, resname_a, resname_b) # any(... calls it with the res names
    return( distance(a,b) <= get(covalentradius, a.element, 0.0) +
               get(covalentradius, b.element, 0.0) )
end

covalent(a::PDBAtom, b::PDBAtom) = covalent(a, b, "", "")

covalent(a::PDBResidue, b::PDBResidue) = any(covalent, a, b, _with_cov)

# Disulphide
# ----------

_issulphurcys(a::PDBAtom, resname_a) = resname_a == "CYS" && a.element == "S"

"Returns `true` if two `CYS`'s `S` are at 2.08 Å or less"
function disulphide(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if _issulphurcys(a, resname_a) && _issulphurcys(b, resname_b)
        return(squared_distance(a,b) <= (2.08 ^ 2))
    end
    return(false)
end

disulphide(a::PDBResidue, b::PDBResidue) = any(disulphide, a, b, _issulphurcys)

# Aromatic-Sulphur
# ----------------

_issulphur(a::PDBAtom) = a.element == "S"

"""
Returns `true` if an sulphur and an aromatic atoms are 5.3 Å or less"
"""
function aromaticsulphur(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if (_issulphur(a) && isaromatic(b, resname_b)) ||
            (_issulphur(b) && isaromatic(a, resname_a))
        return(squared_distance(a,b) <= 28.09) # 28.09 == 5.3 ^ 2
    end
    return(false)
end

_issulphuroraromatic(a::PDBAtom, resname_a) = _issulphur(a) || isaromatic(a, resname_a)

aromaticsulphur(a::PDBResidue, b::PDBResidue) = any(aromaticsulphur, a, b, _issulphuroraromatic)

# Π-Cation
# --------

"""
There's a Π-Cation interaction if a cationic and an aromatic atoms are at 6.0 Å or less
"""
function pication(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if (iscationic(a, resname_a) && isaromatic(b, resname_b)) ||
            (iscationic(b, resname_b) && isaromatic(a, resname_a))
        return(squared_distance(a,b) <= 36.0) # 36.0 == 6.0 ^ 2
    end
    return(false)
end

_iscationicoraromatic(a::PDBAtom, resname_a) = iscationic(a, resname_a) || isaromatic(a, resname_a)

pication(a::PDBResidue, b::PDBResidue) = any(pication, a, b, _iscationicoraromatic)

# Aromatic
# --------

"""
There's an aromatic interaction if centriods are at 6.0 Å or less.
"""
function aromatic(a::PDBResidue, b::PDBResidue)
    threshold = 36.0 # 6.0 ^ 2
    if (a.id.name in _aromatic_res) &&
            (b.id.name in _aromatic_res) &&
            (squared_distance(a, b) <= threshold)
        centres_a = _centre(_get_plane(a))
        centres_b = _centre(_get_plane(b))
        return(any(squared_distance(centroid_a,centroid_b) <= threshold for centroid_a in centres_a, centroid_b in centres_b))
    end
    return(false)
end

# Ionic
# -----

"""
There's an ionic interaction if a cationic and an anionic atoms are at 6.0 Å or less.
"""
function ionic(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if (iscationic(a, resname_a) && isanionic(b, resname_b)) ||
            (iscationic(b, resname_b) && isanionic(a, resname_a))
        return(squared_distance(a,b) <= 36.0) # 36.0 == 6.0 ^ 2
    end
    return(false)
end

_iscationicoranionic(a::PDBAtom, resname_a) = iscationic(a, resname_a) || isanionic(a, resname_a)

ionic(a::PDBResidue, b::PDBResidue) = any(ionic, a, b, _iscationicoranionic)

# Hydrophobic contact
# -------------------

"""
There's an hydrophobic interaction if two hydrophobic atoms are at 5.0 Å or less.
"""
function hydrophobic(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if ishydrophobic(a, resname_a) && ishydrophobic(b, resname_b)
        return(squared_distance(a,b) <= 25.0) # 5.0 ^ 2
    end
    return(false)
end

hydrophobic(a::PDBResidue, b::PDBResidue) = any(hydrophobic, a, b, ishydrophobic)

# Hydrogen bonds
# --------------

function _find_antecedent(res::PDBResidue, a::PDBAtom)
    ids = _hbond_acceptor[(res.id.name, a.atom)]
    _find(a -> a.atom in ids, res.atoms)
end

function _find_h(res::PDBResidue, a::PDBAtom)
    ids = _hbond_donor[(res.id.name, a.atom)]
    _find(a -> a.atom in ids, res.atoms)
end

function _hbond_kernel(donor, acceptor, indices_donor, indices_acceptor)
    @inbounds for i in indices_donor
        don = donor.atoms[i]
        indices_h = _find_h(donor, don)
        if length(indices_h) == 0
            continue
        end
        for j in indices_acceptor
            acc = acceptor.atoms[j]
            indices_ant = _find_antecedent(acceptor, acc)
            if squared_distance(don, acc) <= (3.9^2) && length(indices_ant) != 0
                for k in indices_h
                    hyd = donor.atoms[k]
                    if squared_distance(hyd, acc) <= 6.25 && angle(don, hyd, acc) >= 90.0 # 6.25 == 2.5²
                        for ant in indices_ant
                            if angle(don, acc, acceptor.atoms[ant]) >= 90.0 && angle(hyd, acc, acceptor.atoms[ant]) >= 90.0
                                return(true)
                            end
                        end
                    end
                end
            end
        end
    end
    return(false)
end

function _hydrogenbond_don_acc(donor::PDBResidue, acceptor::PDBResidue)
    if donor != acceptor
        indices_donor = findall(x -> ishbonddonor(x, donor.id.name), donor.atoms)
        indices_acceptor = findall(x -> ishbondacceptor(x, acceptor.id.name), acceptor.atoms)
        if length(indices_donor) != 0 && length(indices_acceptor) != 0
            return(_hbond_kernel(donor, acceptor, indices_donor, indices_acceptor))
        end
    end
    return(false)
end

"""
This function only works if there are hydrogens in the structure.
The criteria for a hydrogen bond are:

- d(Ai, Aj) < 3.9Å
- d(Ah, Aacc) < 2.5Å
- θ(Adon, Ah, Aacc) > 90°
- θ(Adon, Aacc, Aacc-antecedent) > 90°
- θ(Ah, Aacc, Aacc-antecedent) > 90°

Where Ah is the donated hydrogen atom, Adon is the hydrogen bond donor atom,
Aacc is the hydrogen bond acceptor atom and Aacc-antecednt is the atom antecedent to the
hydrogen bond acceptor atom.
"""
hydrogenbond(a::PDBResidue, b::PDBResidue) = _hydrogenbond_don_acc(a,b) || _hydrogenbond_don_acc(b,a)


text:
"""

    DumbAggregate{P<:Primitive} <: Aggregate

A dumb aggregate is the simplest way to store primitives: it's just a list and all the
primitives are checked for intersection one by one. It's mostly used for testing purposes,
as the `BVHAccelerator` will be faster for any reasonable number of primitives.

"""
struct DumbAggregate{P<:Primitive} <: Aggregate
    primitives::Array{P,1}
    bounds::BoundingBox
end

"""

    DumbAggregate(primitives::Array{P,1})

Construct a `DumbAggregate` from a list of primitives, computing the bounding box.

"""
function DumbAggregate(primitives::Array{P,1}) where P<:Primitive
    BB = BoundingBox([world_bounds(p) for p in primitives])
    return DumbAggregate(primitives, BB)
end

"""
    world_bounds(B::DumbAggregate)

Return the bounding box containing all the primitives in the `DumbAggregate`.

"""
world_bounds(A::DumbAggregate) = A.bounds

"""
    intersect(ray::Ray, A::DumbAggregate)

Find nearest intersection, if any, between `ray` and a `DumbAggregate`.

Returns `Intersection` or `nothing`.

"""
function intersect(ray::Ray, A::DumbAggregate)
    best_tmin = Inf
    best_isect = nothing
    if !intersectP(ray, A.bounds)
        return nothing
    end
    for primitive in A.primitives
        best_isect = update_isect(best_isect, intersect(ray, primitive))
    end
    return best_isect
end

"""
    intersectP(ray::Ray, A::DumbAggregate)

Returns `true` if the given ray intersects any primitive in the `DumbAggregate`, else
returns false. This is much faster to do than [`intersect`](@ref), since it's enough to
find any intersection, not specifically the nearest one.

Possibly the most often called function in the entire ray-tracer.

"""
function intersectP(R::Ray, A::DumbAggregate)
    if !intersectP(R, A.bounds)
        return false
    end
    for prim = A.primitives
        if intersectP(R, prim)
            return true
        end
    end
    return false
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using LinearAlgebra: norm, ⋅, ×
using StaticArrays: SVector

export
    sec2rad, rad2sec, normalize_angle, angle,
    dms2rad, rad2dms, sec2deg, deg2sec

"""
    sec2rad(sec)

Convert an angle in arcseconds to radians.

# Example

```jldoctest
julia> sec2rad(3600 * 30)
0.5235987755982988
```
"""
sec2rad(sec) = deg2rad(sec / 3600)

"""
    rad2sec(rad)
Convert an angle in radians to arcseconds.

# Example

```jldoctest
julia> rad2sec(0.5235987755982988)
107999.99999999999
```
"""
rad2sec(rad) = rad2deg(rad) * 3600

sec2deg(sec) = sec / 3600
deg2sec(deg) = deg * 3600

function normalize_angle(angle, center)
    angle - 2π * floor((angle + π - center) / 2π)
end

function dms2rad(deg, arcmin, arcsec)
    deg2rad(deg + arcmin/60 + arcsec/3600)
end

function rad2dms(rad)
    d = rad2deg(rad)
    deg = trunc(d)
    arcmin = trunc((d-deg)*60)
    arcsec = (d-deg-arcmin/60)*3600
    return deg, arcmin, arcsec
end

function angle(v1, v2)
    normprod = norm(v1) * norm(v2)
    if normprod == 0.0
        throw(DomainError())
    else
        v1v2 = v1 ⋅ v2
        threshold = normprod * 0.9999
        if v1v2 >= -threshold && v1v2 <= threshold
            return acos(v1v2 / normprod)
        else
            v3n = norm(v1 × v2)
            return v1v2 >= 0.0 ? (v3n / normprod) : π - asin(v3n / normprod)
        end
    end
end

function azimuth(v)
    atan(v[2], v[1])
end

function elevation(v)
    asin(v[3] / norm(v))
end

function vector_azel(az, el)
    saz, caz = sincos(az)
    sel, cel = sincos(el)
    SVector(caz * cel, saz * cel, sel)
end

function angular_velocity(ψ, δψ, θ, δθ, ϕ, δϕ)
    Ω₁ = δψ * sin(θ) * sin(ϕ) - δθ * cos(ϕ)
    Ω₂ = δψ * sin(θ) * cos(ϕ) + δθ * sin(ϕ)
    Ω₃ = δψ * cos(θ) + δϕ
    [Ω₁, Ω₂, Ω₃]
end


text:
"""
IterTolManager: convergence based on iterations and tolerance

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""

immutable IterTolManager <: IterationManager
    tol::Float64
    maxiter::Int
    verbose::Bool
    print_skip::Int
    print_prefix::ASCIIString
end
typealias DefaultManager IterTolManager

IterTolManager(tol::Float64, maxiter::Int) =
    IterTolManager(tol, maxiter, true, div(maxiter, 5), "")

function IterTolManager(;tol::Float64=1e-10, maxiter::Int=1000,
                         verbose::Bool=true, print_skip=div(maxiter, 5),
                         print_prefix="")
    IterTolManager(tol, maxiter, verbose, print_skip, print_prefix)
end

function Base.writemime(io::IO, ::MIME"text/plain", dm::DefaultManager)
    @printf(io,
            "DefaultManager. %2.4e tolerance level, %i max iterations (%s)",
            dm.tol, dm.maxiter, dm.verbose ? "verbose" : "not verbose")
end

finished(mgr::DefaultManager, istate::IterationState) =
    num_iter(istate) > mgr.maxiter || norm(istate) <= mgr.tol


text:
# Unit tests for implementations/molly/molly_result.jl

@testset "molly_result.jl" begin
    particles = [
        AtomsBase.Atom(:Ar, (@SVector [7.0, 7.0, 7.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"; meta=:data),
        AtomsBase.Atom(:Ar, (@SVector [7.0, 7.0, 21.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"; hello="world"),
        AtomsBase.Atom(:Ar, (@SVector [7.0, 21.0, 7.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [7.0, 21.0, 21.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21.0, 7.0, 7.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21.0, 7.0, 21.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21.0, 21.0, 7.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21.0, 21.0, 21.0])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au")
    ]
    box = [[28.0, 0.0, 0.0], [0.0, 28.0, 0.0], [0.0, 0.0, 28.0]]u"bohr"
    system = periodic_system(particles, box)

    simulator = MollySimulator(400, 10, t₀=1000, coupling=Molly.AndersenThermostat(94.4u"K", 0.1u"ps"))
    simulator2 = MollySimulator(400, 9, t₀=1000, stride=2)
    simulator3 = MollySimulator(400, 10, t₀=1000, stride=2)

    potential = InteratomicPotentials.LennardJones(1.657e-21u"J", 0.34u"nm", 0.765u"nm", [:Ar])

    result = simulate(system, simulator, potential)
    result2 = simulate(system, simulator2, potential)
    result3 = simulate(system, simulator3, potential)

    @test get_time_range(result) == (1000:400:5000)u"ħ_au / hartree"
    @test get_time_range(result2) == (1400:800:4600)u"ħ_au / hartree"
    @test get_time_range(result3) == (1400:800:5400)u"ħ_au / hartree"

    @test get_num_bodies(result) == 8
    @test get_bounding_box(result) == box
    @test get_boundary_conditions(result) == [Periodic(), Periodic(), Periodic()]

    @test reference_temperature(result) == 94.4u"K"
    @test ismissing(reference_temperature(result2))

    @test get_positions(result) isa AbstractVector{<:StaticVector{3,<:Unitful.Length}}
    @test get_velocities(result) isa AbstractVector{<:StaticVector{3,<:Unitful.Velocity}}
    @test get_particles(result) isa AbstractVector{<:AtomsBase.Atom}
    @test get_system(result) isa System{3}

    @test all(all(0.0u"bohr" ≤ c < 28.0u"bohr" for c ∈ p) for p ∈ get_positions(result))
    @test get_particles(result)[1].meta == :data && get_particles(result)[2].hello == "world"

    @test all(get_positions(result, t) isa AbstractVector{<:StaticVector{3,<:Unitful.Length}} for t ∈ 1:11)
    @test all(get_velocities(result, t) isa AbstractVector{<:StaticVector{3,<:Unitful.Velocity}} for t ∈ 1:11)
    @test all(get_particles(result, t) isa AbstractVector{<:AtomsBase.Atom} for t ∈ 1:11)

    @test all(all(all(0.0u"bohr" ≤ c < 28.0u"bohr" for c ∈ p) for p ∈ get_positions(result, t)) for t ∈ 1:11)
    @test all(get_particles(result, t)[1].meta == :data && get_particles(result, t)[2].hello == "world" for t ∈ 1:11)

    @test Atomistic.temperature(result) isa Unitful.Temperature
    @test Atomistic.kinetic_energy(result) isa Unitful.Energy
    @test Atomistic.potential_energy(result) isa Unitful.Energy

    @test all(isapprox(Atomistic.total_energy(result3, t), Atomistic.total_energy(result3), rtol=0.1) for t ∈ 2:6)
end


text:
#
# CR3BP systems in our solar system
#

"""
The Sun-Earth CR3BP system.
"""
const SunVenus = CR3BPParameters(get_μ(Sun) * massparamunit(Sun), get_μ(Venus) * massparamunit(Venus), 108.2e6u"km"; primary=:Sun, secondary=:Venus)

"""
The Sun-Earth CR3BP system.
"""
const SunEarth = CR3BPParameters(get_μ(Sun) * massparamunit(Sun), get_μ(Earth) * massparamunit(Earth), 1.0u"AU"; primary=:Sun, secondary=:Earth)

"""
The Earth-Moon CR3BP system.
"""
const EarthMoon = CR3BPParameters(get_μ(Earth) * massparamunit(Earth), get_μ(Moon) * massparamunit(Moon), 384400u"km"; primary=:Earth, secondary=:Moon)

"""
The Sun-Mars CR3BP system.
"""
const SunMars = CR3BPParameters(get_μ(Sun) * massparamunit(Sun), get_μ(Mars) * massparamunit(Mars), 227.9e6u"km"; primary=:Sun, secondary=:Mars)

"""
The Sun-Jupiter CR3BP system.
"""
const SunJupiter = CR3BPParameters(get_μ(Sun) * massparamunit(Sun), get_μ(Jupiter) * massparamunit(Jupiter), 778.6e6u"km"; primary=:Sun, secondary=:Jupiter)

"""
The Sun-Saturn CR3BP system.
"""
const SunSaturn = CR3BPParameters(get_μ(Sun) * massparamunit(Sun), get_μ(Saturn) * massparamunit(Saturn), 1433.5e6u"km"; primary=:Sun, secondary=:Saturn)


text:
export StaticModel

"""
    struct StaticModel{I,O} <: VariableTimeModel{I,O,O} end

A static model represents a variable that never changes its value.
The value is setup through an sfunc created by make_initial.
At any time point, it simply copies its previous value.
Because it can be used flexibly, we make it a subtype of VariableTimeModel.
"""
struct StaticModel{I,O} <: VariableTimeModel{I,Tuple{O},O}
    sf::SFunc{I,O}
end

function make_initial(m::StaticModel, time)
    return m.sf
end

function make_transition(::StaticModel{I,O}, parenttimes, time) where {I,O}
    return Det(Tuple{O}, O, x -> x[1])
end


text:
"""
    ocean_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{FreeSlip}, ::SWModel)

apply free slip boundary condition for velocity
sets reflective ghost point
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{FreeSlip},
    ::SWModel,
    ::TurbulenceClosure,
    q⁺,
    a⁺,
    n⁻,
    q⁻,
    a⁻,
    t,
    args...,
)
    q⁺.η = q⁻.η

    V⁻ = @SVector [q⁻.U[1], q⁻.U[2], -0]
    V⁺ = V⁻ - 2 * n⁻ ⋅ V⁻ .* SVector(n⁻)
    q⁺.U = @SVector [V⁺[1], V⁺[2]]

    return nothing
end

"""
    ocean_boundary_state!(::Union{NumericalFluxGradient, NumericalFluxSecondOrder}, ::Impenetrable{FreeSlip}, ::SWModel)

no second order flux computed for linear drag
"""
ocean_boundary_state!(
    ::Union{NumericalFluxGradient, NumericalFluxSecondOrder},
    ::VelocityBC,
    ::SWModel,
    ::LinearDrag,
    _...,
) = nothing

"""
    ocean_boundary_state!(::NumericalFluxGradient, ::Impenetrable{FreeSlip}, ::SWModel)

apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function ocean_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{FreeSlip},
    ::SWModel,
    ::ConstantViscosity,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    t,
    args...,
)
    V⁻ = @SVector [Q⁻.U[1], Q⁻.U[2], -0]
    V⁺ = V⁻ - n⁻ ⋅ V⁻ .* SVector(n⁻)
    Q⁺.U = @SVector [V⁺[1], V⁺[2]]

    return nothing
end

"""
    shallow_normal_boundary_flux_second_order!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::SWModel)

apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function ocean_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::SWModel,
    ::ConstantViscosity,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
    args...,
)
    Q⁺.U = Q⁻.U
    D⁺.ν∇U = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    ocean_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{NoSlip}, ::SWModel)

apply no slip boundary condition for velocity
sets reflective ghost point
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{NoSlip},
    ::SWModel,
    ::TurbulenceClosure,
    q⁺,
    α⁺,
    n⁻,
    q⁻,
    α⁻,
    t,
    args...,
)
    q⁺.η = q⁻.η
    q⁺.U = -q⁻.U

    return nothing
end

"""
    ocean_boundary_state!(::NumericalFluxGradient, ::Impenetrable{NoSlip}, ::SWModel)

apply no slip boundary condition for velocity
set numerical flux to zero for U
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{NoSlip},
    ::SWModel,
    ::ConstantViscosity,
    q⁺,
    α⁺,
    n⁻,
    q⁻,
    α⁻,
    t,
    args...,
)
    FT = eltype(q⁺)
    q⁺.U = @SVector zeros(FT, 2)

    return nothing
end

"""
    ocean_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::SWModel)

apply no slip boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for U
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::SWModel,
    ::ConstantViscosity,
    q⁺,
    σ⁺,
    α⁺,
    n⁻,
    q⁻,
    σ⁻,
    α⁻,
    t,
    args...,
)
    q⁺.U = -q⁻.U
    σ⁺.ν∇U = σ⁻.ν∇U

    return nothing
end

"""
    ocean_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{FreeSlip}, ::SWModel)

no mass boundary condition for penetrable
"""
ocean_boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{FreeSlip},
    ::SWModel,
    ::ConstantViscosity,
    _...,
) = nothing

"""
    ocean_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{FreeSlip}, ::SWModel)

apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function ocean_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Penetrable{FreeSlip},
    ::SWModel,
    ::ConstantViscosity,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
    args...,
)
    Q⁺.U = Q⁻.U
    D⁺.ν∇U = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    ocean_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Impenetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Impenetrable{KinematicStress},
    shallow::SWModel,
    turb::TurbulenceClosure,
    args...,
)
    return ocean_boundary_state!(
        nf,
        Impenetrable(FreeSlip()),
        shallow,
        turb,
        args...,
    )
end

"""
    ocean_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{KinematicStress},
    shallow::SWModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.U = Q⁻.U
    D⁺.ν∇U = n⁻ * kinematic_stress(shallow.problem, A⁻.y, 1000)'
    # applies windstress for now, will be fixed in a later PR

    return nothing
end

"""
    ocean_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
applies free slip conditions for first-order and gradient fluxes
"""
function ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{KinematicStress},
    shallow::SWModel,
    turb::TurbulenceClosure,
    args...,
)
    return ocean_boundary_state!(
        nf,
        Penetrable(FreeSlip()),
        shallow,
        turb,
        args...,
    )
end

"""
    ocean_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{KinematicStress}, ::HBModel)

apply kinematic stress boundary condition for velocity
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function ocean_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Penetrable{KinematicStress},
    shallow::SWModel,
    ::TurbulenceClosure,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.u = Q⁻.u
    D⁺.ν∇u = n⁻ * kinematic_stress(shallow.problem, A⁻.y, 1000)'
    # applies windstress for now, will be fixed in a later PR

    return nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
eq=Der(:a*:x,:x)≖:a
ex=Der(3*:x,:x)
md=matches(ex,eq.lhs)
@test !isempty(md)
m=matches(ex,eq)
@test m[1]==3

ex=Der(3*:x,:y)
md=matches(ex,eq.lhs)
@test isempty(md)

c1=Der(3,:x);pat1=Equation(Der(:a,:x),0)
@test c1&pat1==0
c2=Der(3*:x,:x);pat2=Der(:a*:x,:x)≖:a
@test c2&pat2==3
d=Der(:x^3,:x);pat=relations["Der"][3];
@test randeval(d&pat)==randeval(3*Pow(:x,2))
@test randeval((Der(:x^:n,:x)-Der(-0.1*:x^:m,:x)+1/:a*Der(:a*sqrt(:x),:x))&relations["Der"])==randeval(:n*Pow(:x,:n-1) + 0.1*:m* Pow(:x,:m -1) + 0.5*Pow(:x,-0.5))