text:
Trajectories.trajectory(tt, v::Type, D, ::Val{true}) = trajectory(tt, v, D)
Trajectories.trajectory(tt, v::Type, D, ::Val{false}) = trajectory(tt, v)

"""
    Trajectories.trajectory(tt, v::Type, D::Number)

Create a Trajectory with mutable states of dimension `D` along the time
collection `tt`.
"""
function Trajectories.trajectory(tt, v::Type, D::Number)
    ismutable(v) == Val(false) && return trajectory(tt, v)
    trajectory(collect(tt), [zeros(eltype(v),D) for _ in tt])
end

"""
    Trajectories.trajectory(tt, v::Type)

Create a Trajectory with immutable states along the time collection `tt`.
"""
function Trajectories.trajectory(tt, v::Type)
    ismutable(v)::Val{false}
    trajectory(collect(tt), zeros(v, length(tt)))
end

function Trajectories.trajectory(
        tt::AbstractArray{<:AbstractArray},
        P::DiffusionProcess,
        v::Type=default_type(P),
        w::Type=default_wiener_type(P),
    )
    (
        process = [_process_traj(t, P, v) for t in tt],
        wiener = [_wiener_traj(t, P, w) for t in tt],
    )
end

function Trajectories.trajectory(
        tt,
        P::DiffusionProcess,
        v::Type=default_type(P),
        w::Type=default_wiener_type(P)
    )
    (
        process = _process_traj(tt, P, v),
        wiener = _wiener_traj(tt, P, w)
    )
end

function _process_traj(tt, ::DiffusionProcess{T,DP}, v::Type) where {T,DP}
    trajectory(tt, v, eval(DP))
end

function _wiener_traj(tt, ::DiffusionProcess{T,DP,DW}, v::Type) where {T,DP,DW}
    trajectory(tt, v, eval(DW))
end

@recipe function f(tr::Trajectory, ::Val{:vs_time}; coords=1:length(tr.x[1]))
    for i in coords
        @series begin
            seriestype := :path
            tr.t, map(x->x[i], tr.x)
        end
    end
end

@recipe function f(tr::Trajectory, ::Val{:x_vs_y}; coords=1:length(tr.x[1]))
    @assert length(coords) == 2
    @series begin
        seriestype := :path
        map(x->x[coords[1]], tr.x), map(x->x[coords[2]], tr.x)
    end
end

@recipe function f(tr::Trajectory; vs_time=true, coords=1:length(tr.x[1]))
    plot(tr, vs_time ? Val(:vs_time) : ValVal(:x_vs_y); coords=coords)
end


text:
function FieldFunctionalDetector(arg0::Field)
    return FieldFunctionalDetector((Field,), arg0)
end

function g(obj::FieldFunctionalDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_function(obj::FieldFunctionalDetector)
    return jcall(obj, "getFunction", Function, ())
end

function with_function(obj::FieldFunctionalDetector, arg0::Function)
    return jcall(obj, "withFunction", FieldFunctionalDetector, (Function,), arg0)
end

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

using AbaqusReader: element_has_type, element_has_nodes, parse_abaqus,
                    parse_section, abaqus_read_mesh

datadir = first(splitext(basename(@__FILE__)))

@testset "parse beam.inp mesh" begin
    fn = joinpath(datadir, "beam.inp")
    mesh = open(parse_abaqus, fn)
    @test length(mesh["nodes"]) == 298
    @test length(mesh["elements"]) == 120
    @test length(mesh["element_sets"]["Body1"]) == 120
    @test length(mesh["node_sets"]["SUPPORT"]) == 9
    @test length(mesh["node_sets"]["LOAD"]) == 9
    @test length(mesh["node_sets"]["TOP"]) == 83
end

@testset "test read element section" begin
    data = """*ELEMENT, TYPE=C3D10, ELSET=BEAM
    1,       243,       240,       191,       117,       245,       242,       244,
    1,         2,       196
    2,       204,       199,       175,       130,       207,       208,       209,
    3,         4,       176
    """
    data = split(data, "\n")
    model = Dict{String, Any}()
    model["node_sets"] = Dict{String, Vector{Int}}()
    model["elements"] = Dict{Int, Vector{Int}}()
    model["element_sets"] = Dict{String, Vector{Int}}()
    model["element_types"] = Dict{Int, Symbol}()
    parse_section(model, data, :ELEMENT, 1, 5, Val{:ELEMENT})
    @test length(model["elements"]) == 2
    @test model["element_sets"]["BEAM"] == [1, 2]
    @test model["elements"][1] == [243, 240, 191, 117, 245, 242, 244, 1, 2, 196]
    @test model["elements"][2]== [204, 199, 175, 130, 207, 208, 209, 3, 4, 176]
end

@testset "parse cube_tet4.inp mesh" begin
    fn = joinpath(datadir, "cube_tet4.inp")
    mesh = open(parse_abaqus, fn)
    @test length(mesh["nodes"]) == 10
    @test length(mesh["elements"]) == 17
    @test haskey(mesh["elements"], 1)
    @test mesh["elements"][1] == [8, 10, 1, 2]
    @test mesh["element_types"][1] == :Tet4
    @test haskey(mesh["node_sets"], "SYM12")
    @test haskey(mesh["element_sets"], "CUBE")
    @test haskey(mesh["surface_sets"], "LOAD")
    @test haskey(mesh["surface_sets"], "ORDER")
    @test length(mesh["surface_sets"]["LOAD"]) == 2
    @test mesh["surface_sets"]["LOAD"][1] == (16, :S1)
    @test mesh["surface_types"]["LOAD"] == :ELEMENT
    @test length(Set(map(size, values(mesh["nodes"])))) == 1
end

@testset "parse nodes from abaqus .inp file to Mesh (NX export)" begin
    fn = joinpath(datadir, "nx_export_problem.inp")
    mesh = open(parse_abaqus, fn)
    @test length(mesh["nodes"]) == 3
end

@testset "parse abaqus .inp created using hypermesh" begin
    fn = joinpath(datadir, "hypermesh_model.inp")
    mesh = abaqus_read_mesh(fn)
    @test length(mesh["nodes"]) == 2
end

@testset "find element types and nodes" begin
    @test element_has_nodes(Val{:C3D4}) == 4
    @test element_has_type(Val{:C3D4}) == :Tet4
    @test element_has_nodes(Val{:C3D8}) == 8
    @test element_has_type(Val{:C3D8}) == :Hex8
    @test element_has_nodes(Val{:C3D10}) == 10
    @test element_has_type(Val{:C3D10}) == :Tet10
    @test element_has_nodes(Val{:C3D20}) == 20
    @test element_has_nodes(Val{:C3D20E}) == 20
    @test element_has_nodes(Val{:S3}) == 3
    @test element_has_type(Val{:S3}) == :Tri3
    @test element_has_nodes(Val{:STRI65}) == 6
    @test element_has_type(Val{:STRI65}) == :Tri6
    @test element_has_nodes(Val{:CPS4}) == 4
    @test element_has_type(Val{:CPS4}) == :Quad4
end

@testset "use GENERATE keyword" begin
    data = """
*NSET, NSET=testgen, GENERATE
7,13,2
"""
    fn = tempname() * ".inp"
    open(fn, "w") do fid write(fid, data) end
    mesh = open(parse_abaqus, fn)
    @test mesh["node_sets"]["testgen"] == [7, 9, 11, 13]
end

@testset "parse ELSET" begin
    lines = ["*ELSET, ELSET=TEST1", "1"]
    mesh = Dict("element_sets" => Dict{String, Vector{Int}}())
    parse_section(mesh, lines, :ELSET, 1, 2, Val{:ELSET})
    @test mesh["element_sets"]["TEST1"] == [1]
end

@testset "parse tet4.inp" begin
    fn = joinpath(datadir, "tet4.inp")
    mesh = abaqus_read_mesh(fn)
    @test length(mesh["nodes"]) == 116
end


text:
using Test
using LogicCircuits

@testset "node stats tests" begin

    lc = load_logical_circuit(zoo_lc_file("little_4var.circuit"));

    lstats = leaf_stats(lc);
    istats = inode_stats(lc);
    nstats = node_stats(lc);

    @test !(TrueNode in keys(lstats));
    @test !(FalseNode in keys(lstats));
    @test lstats[LiteralNode] == 16;

    @test istats[(⋀Node, 2)] == 9;
    @test istats[(⋁Node, 1)] == 10;
    @test istats[(⋁Node, 4)] == 2;

    for t in keys(lstats)
        @test lstats[t] == nstats[t]
    end
    for t in keys(istats)
        @test istats[t] == nstats[t]
    end

    @test num_nodes(lc) == 37;
    @test num_variables(lc) == 4;
    @test num_edges(lc) == 36;

end


text:
@testset "Zygote" begin
    grads = IdDict()
    grads[:x] = [-3.0 0.0 0.0; 4.0 0.0 0.0]
    ps = Zygote.Params([:x])
    gs = Zygote.Grads(grads, ps)
    clip_by_global_norm!(gs, ps, 4.0f0)
    @test isapprox(gs[:x], [-2.4 0.0 0.0; 3.2 0.0 0.0])

    gs.grads[:x] = [1.0 0.0 0.0; 1.0 0.0 0.0]
    clip_by_global_norm!(gs, ps, 4.0f0)
    @test isapprox(gs[:x], [1.0 0.0 0.0; 1.0 0.0 0.0])

    gs.grads[:x] = [0.0 0.0 0.0; 0.0 0.0 0.0]
    clip_by_global_norm!(gs, ps, 4.0f0)
    @test isapprox(gs[:x], [0.0 0.0 0.0; 0.0 0.0 0.0])
end

@testset "Distributions" begin
    @test isapprox(logpdf(Normal(), 2), normlogpdf(0, 1, 2))
    @test isapprox(
        logpdf.([Normal(), Normal()], [2, 10]),
        normlogpdf([0, 0], [1, 1], [2, 10]),
    )

    # Test numeric stability for 0 sigma
    @test isnan(normlogpdf(0, 0, 2, ϵ = 0))
    @test !isnan(normlogpdf(0, 0, 2))

    if CUDA.functional()
        cpu_grad = Zygote.gradient([0.2, 0.5]) do x
            sum(logpdf.([Normal(1, 0.1), Normal(2, 0.2)], x))
        end
        gpu_grad = Zygote.gradient(cu([0.2, 0.5])) do x
            sum(normlogpdf(cu([1, 2]), cu([0.1, 0.2]), x))
        end
        @test isapprox(cpu_grad[1], gpu_grad[1] |> Array)
    end
end


text:
import Setfield

struct Latch
    _data::UInt32
end

Base.UInt32(l::Latch) = getfield(l, :_data)

"""
True if a particle crossed scoring plane multiple times.
"""
const MASK_MULTICROSS       = 0b10000000000000000000000000000000

"""
Encodes the charge of a particle -1,0,+1.
"""
const MASK_CHARGE        = 0b01100000000000000000000000000000

"""
The bit region in which a particle was created. 0 for primary particles.
"""
const MASK_CREATION      = 0b00011111000000000000000000000000

"""
Marks the bit regions that were traversed by particle.
"""
const MASK_VISITED       = 0b00000000111111111111111111111110

"""
True if bremsstrahlung or annihilation even occurs in history.
"""
const MASK_BREMS         = 0b00000000000000000000000000000001

struct BitRegions <: AbstractVector{Bool}
    data::UInt32
end

Base.UInt32(o::BitRegions) = o.data

macro BitRegions(args...)
    pattern = zero(UInt32)
    for i in args
        pattern = setbit(pattern, true, i)
    end
    pattern
    esc(:(BitRegions($pattern)))
end

function Base.size(o::BitRegions)
    (23,)
end
function Base.getindex(o::BitRegions, index)
    @boundscheck checkbounds(o, index)
    data = UInt32(o)
    getbit(data, index)
end

function show_bitregions(io::IO, r::BitRegions)
    data = UInt32(r)
    print(io, "@BitRegions", "(")
    for i in eachindex(r)
        if r[i]
            print(io, i, ", ")
        end
    end
    print(io, ")")
end
Base.show(io::IO, ::MIME"text/plain", r::BitRegions) = show_bitregions(io, r)
Base.show(io::IO, r::BitRegions) = show_bitregions(io, r)

function Base.propertynames(latch::Latch)
    (:multicross, :charge, :creation, :visited, :brems)
end

function Base.getproperty(o::Latch, s::Symbol)
    if s == :multicross
        get_multicross(o)
    elseif s == :charge
        get_charge(o)
    elseif s == :creation
        get_creation(o)
    elseif s == :visited
        get_visited(o)
    elseif s == :brems
        get_brems(o)
    else
        @argcheck s in propertynames(o)
    end
end

function Setfield.setproperties(o::Latch, props::NamedTuple)
    Latch(;
        multicross=get(props, :multicross, o.multicross),
        charge=get(props, :charge, o.charge),
        creation=get(props, :creation, o.creation),
        visited=get(props, :visited, o.visited),
        brems=get(props, :brems, o.brems),
    )
end

function Latch(;
        multicross::Bool=false,
        charge::Int,
        creation::Int=0, # actually bitregion 0 does not exist
        visited::BitRegions=BitRegions(zero(UInt32)),
        brems::Bool=false,
    )
    l = zero(UInt32)
    l += brems_pattern(brems)
    l += visited_pattern(visited)
    l += creation_pattern(creation)
    l += charge_pattern(charge)
    l += multicross_pattern(multicross)
    Latch(l)
end

function brems_pattern(brems::Bool)
    l = zero(UInt32)
    setbit(l, brems, 0)
end
function visited_pattern(o::BitRegions)
    o.data
end
function creation_pattern(region)
    UInt32(region) << 24
end
function charge_pattern(charge)
    if charge == 0
        UInt32(0)
    elseif charge == -1
        UInt32(1 << 30)
    elseif charge == 1
        UInt32(1 << 29)
    else
        @argcheck charge in (-1, 0, 1)
    end
end
function multicross_pattern(mc::Bool)
    l = zero(UInt32)
    setbit(l, mc, 31)
end
function get_creation(latch::Latch)
    l = UInt32(latch) & MASK_CREATION
    Int(l >> 24)
end
function get_visited(latch::Latch)
    BitRegions(UInt32(latch) & MASK_VISITED)
end
function get_brems(latch::Latch)
    Bool(UInt32(latch) & MASK_BREMS)
end
function get_multicross(latch::Latch)
    getbit(UInt32(latch), 31)
end
function get_charge(latch::Latch)
    l = UInt32(latch)
    l &= MASK_CHARGE
    if l == charge_pattern(0)
        0
    elseif l == charge_pattern(-1)
        -1
    elseif l == charge_pattern(1)
        1
    else
        msg = "Unsupported charge pattern $l"
        throw(ArgumentError(msg))
    end
end

Base.show(io::IO, latch::Latch) = kwshow(io, latch)

#### EGSParticle
struct EGSParticle{ZLAST <: Union{Nothing,Float32}}
    _latch::Latch
    _E::Float32 # sign bit new histrory
    _x::Float32
    _y::Float32
    _u::Float32
    _v::Float32
    _weight::Float32 # sign bit sign w
    _zlast::ZLAST
end

Format(p::EGSParticle) = FormatEGS()

# TODO HACKY
Base.isapprox(p1::EGSParticle, p2::EGSParticle) = p1 === p2

function Base.propertynames(o::EGSParticle)
    (:latch, :new_history, :E, :x, :y, :u, :v, :w, :weight, :zlast)
end

function Setfield.setproperties(o::EGSParticle, props::NamedTuple)
    EGSParticle(
        get(props, :latch, o.latch),
        get(props, :new_history, o.new_history),
        get(props, :E, o.E),
        get(props, :x, o.x),
        get(props, :y, o.y),
        get(props, :u, o.u),
        get(props, :v, o.v),
        get(props, :w, o.w),
        get(props, :weight, o.weight),
        get(props, :zlast, o.zlast),
    )

end

function EGSParticle(latch::Latch, new_history::Bool, E, x, y, u, v, w, weight, zlast)
    @argcheck weight >= 0
    @argcheck E >= 0
    @argcheck Float32(u^2 + v^2 + w^2) ≈ 1
    charge = latch.charge
    E_rest = rest_energy_by_charge(charge)
    E_tot = kin2total(Float32(E), E_rest)
    _E = Float32((-1)^new_history * E_tot)
    _weight = Float32(sign(w) * weight)
    EGSParticle(latch, _E, Float32(x), Float32(y), Float32(u), Float32(v), _weight, zlast)
end

function EGSParticle(;latch , new_history=true, E,x,y,u,v,w,weight=1f0,zlast=nothing)
    EGSParticle(latch, new_history, E, x, y, u, v, w, weight, zlast)
end

Base.show(io::IO, o::EGSParticle) = kwshow(io, o)

@inline function Base.getproperty(o::EGSParticle, s::Symbol)
    if s == :latch
        get_latch(o)
    elseif s == :E
        get_E(o)
    elseif s == :x
        get_x(o)
    elseif s == :y
        get_y(o)
    elseif s == :u
        get_u(o)
    elseif s == :v
        get_v(o)
    elseif s == :w
        get_w(o)
    elseif s == :weight
        get_weight(o)
    elseif s == :zlast
        get_zlast(o)
    elseif s == :new_history
        get_new_history(o)
    else
        throw(ErrorException("$o does not have property $s"))
    end
end

function kin2total(Ekin::Float32, E_rest::Float64)
    Float32(Float64(Ekin) + E_rest)
end

function total2kin(Etotal::Float32, E_rest::Float64)
    Float32(Float64(Etotal) - E_rest)
end

function rest_energy_by_charge(charge::Int)
    ifelse(charge == 0, 0., 0.511)
end

get_latch(o::EGSParticle) = getfield(o, :_latch)
function get_E(o::EGSParticle)
    E_tot = abs(getfield(o, :_E))
    E_rest = rest_energy_by_charge(o.latch.charge)
    total2kin(E_tot, E_rest)
end

get_x(o::EGSParticle) = getfield(o, :_x)
get_y(o::EGSParticle) = getfield(o, :_y)
get_u(o::EGSParticle) = getfield(o, :_u)
get_v(o::EGSParticle) = getfield(o, :_v)

get_weight(o::EGSParticle) = abs(getfield(o, :_weight))
get_zlast(o::EGSParticle)  = getfield(o, :_zlast)
get_new_history(o::EGSParticle) = signbit(getfield(o, :_E))

function get_w(o::EGSParticle)
    u = Float64(get_u(o))
    v = Float64(get_v(o))
    sign_w = sign(getfield(o, :_weight))
    w64_square = 1 - u^2 - v^2
    w64 = if w64_square < 0
        @assert isapprox(w64_square,0,atol=100eps(Float32))
        0.
    else
        sign_w * sqrt(1 - u^2 - v^2)
    end
    Float32(w64)
end

function particle_type(p::EGSParticle)
    c = p.latch.charge
    if c == -1
        electron
    elseif c == 0
        photon
    else
        @assert c == 1
        positron
    end
end

isphoton(o::EGSParticle) = o.latch.charge == 0
iselectron(o::EGSParticle) = o.latch.charge == -1
ispositron(o::EGSParticle) = o.latch.charge == 1


text:
"""
simple!(FV,x)
This is the function for Figure 2.1
It also shows up in CI
"""
function simple!(FV, x)
    FV[1] = x[1] * x[1] + x[2] * x[2] - 2.0
    FV[2] = exp(x[1] - 1) + x[2] * x[2] - 2.0
end

function jsimple!(JacV, FV, x)
    JacV[1, 1] = 2.0 * x[1]
    JacV[1, 2] = 2.0 * x[2]
    JacV[2, 1] = exp(x[1] - 1)
    JacV[2, 2] = 2 * x[2]
end

"""
JVsimple(v, FV, x)

Jacobian-vector product for simple!. There is, of course, no reason
to use Newton-Krylov for this problem other than CI or demonstrating
how to call nsoli.jl.
"""
function JVsimple(v, FV, x)
    jvec = zeros(2)
    jvec[1] = 2.0 * x' * v
    jvec[2] = v[1] * exp(x[1] - 1.0) + 2.0 * v[2] * x[2]
    return jvec
end


text:
@deprecate k_nearest_neighbour{T <: AbstractFloat}(tree::KDTree, point::Vector{T}, k::Int) knn(tree::KDTree, point::Vector{T}, k::Int)

@deprecate query_ball_point{T <: AbstractFloat}(tree::KDTree{T}, point::Vector{T}, radius::T) inball(tree::KDTree{T}, point::Vector{T}, radius::T)

@deprecate  KDTree{T <: AbstractFloat}(data::Matrix{T}, ls::Int, reord::Bool = true)  KDTree(data, leafsize = ls, reorder = reord)


text:
# Web Scraping in Julia
# Getting nightly test information from web.
#
#
# Hongyang Zhou, hyzhou@umich.edu 11/22/2018 initial version
# 05/08/2019 add timings for different tests; add parameter input file

# I need a way to check if the test correctly passed!
# The logic imposed is really fragile...

using HTTP, Gumbo, Dates, DelimitedFiles, ProgressMeter, TOML
using Plots
using Statistics
using Printf
pyplot() # Set the backend to Plotly
#plotly()

"""
   getFuncTimings(dates,[testsuite,testcase])

Obtain the timings from html files for SWMF tests.

# Arguments
- `dates::Vector{Date}`: date interval.
- `testsuite::String="mesh"`: name of machine/platform.
"""
function getFuncTimings(dates::Array{Date,1}; testsuite::String = "mesh",
   verbose::Bool=false)
   url_front = "http://herot.engin.umich.edu/~gtoth/SWMF_TEST_RESULTS/"
   url_end   = "/test_swmf.html"
   ndates = size(dates,1)
   urls = Array{String}(undef,ndates)
   timings = zeros(ndates)
   tTotal = Vector{Float64}(undef,ndates)

   # Obtain data from html
   @showprogress 1 "Fetching data..." for i=1:ndates
      urls[i] = url_front*Dates.format(dates[i],"yyyy/mm/dd/")*testsuite*url_end

      try
         r = HTTP.request("GET", urls[i]; retries=0,verbose=0)
         rawText = String(r.body)[end-2850:end]
         rawText = rawText[ collect(findfirst("\n000",rawText))[1]+1:end ]

         timeTable = split(rawText,"\n")
         for row in timeTable[1:36]
            cols = split(row)
            try
               timings[i] = parse(Float64,cols[3])
            catch
               verbose && println("Test",cols[1]," failed on ",dates[i])
               continue
            end
         end

         tTotal[i] = parse(Float64,split(timeTable[38])[3])
      catch
         continue
      end

   end

   # Clear all large/small values
   for i=1:length(tTotal)
      if tTotal[i] < 1e-5 || tTotal[i] > 1e4
         tTotal[i] = 0.0
      elseif isnan(tTotal[i])
         tTotal[i] = 0.0
      end
   end

   return timings,tTotal
end

"""
   getTimings(dates,[testsuite,testcase])

Obtain the timings from html log files for SWMF tests.

# Arguments
- `dates::Vector{Date}`: date interval.
- `testsuite::String="mesh"`: name of machine/platform.
"""
function getTimings(dates::Array{Date,1};
   testsuite::String="gfortran",testcase::String="shocktube",
   verbose::Bool=false)
# Find test_shocktube_run
# PostProc.pl: TIMINGS from runlog (init, run) 0.01 1.28
url_front = "http://herot.engin.umich.edu/~gtoth/SWMF_TEST_RESULTS/"
url_end   = "/test_swmf.log"
ndates = size(dates,1)
urls = Array{String}(undef,ndates)
timings = zeros(ndates)
nFail,nError = 0,0

# Obtain data from log
@showprogress 1 "Fetching data..." for i=1:ndates
   urls[i] = url_front*Dates.format(dates[i],"yyyy/mm/dd/")*testsuite*url_end

   try
      r = HTTP.request("GET", urls[i]; retries=0,verbose=0)
      rawText = String(r.body)
      start_  = collect(findfirst(testcase*"_run\n",rawText))[end]
      end_    = collect(findfirst(testcase*".diff\n",rawText))[1] + 100

      if end_ - start_ > 5000 # something went wrong if indexes diff to much
         timings[i] = 0.0
         nFail += 1
         continue
      end

      rawText = rawText[start_:end_]

      start_  = findfirst("(init, run)",rawText)[end]+7
      end_    = start_ + 3

      timings[i] = parse(Float64,rawText[start_:end_])

      end_ = collect(findfirst(testcase*".diff\n",rawText))[end]
      end_ = collect(findnext(testcase*".diff\n",rawText,end_))[1]

      start_ = end_ - 21
      end_ -= 19

      if DoCheckResult && parse(Int,rawText[start_:end_]) != 0
         timings[i] = 0.0
         nFail += 1
      end

   catch
      nError += 1
      continue
   end

end

# Remove outliers
tMean,tStd = mean(timings),std(timings)
for i=1:ndates
   if timings[i] > tMean + 5*tStd
      timings[i] = 0.0
   end
end

@printf("test pass rate = %4.2f%%, no data rate = %4.2f%%",
   (1.0-nFail/(ndates-nError))*100, nError/ndates*100)

return timings

end

paramIn = TOML.parsefile("PARAM.toml")
verbose   = paramIn["Param"]["verbose"]
testsuite = paramIn["Param"]["testsuite"]
testcase  = paramIn["Param"]["testcase"]
DoCheckResult = paramIn["Param"]["DoCheckResult"]

dateStart = Date(paramIn["Param"]["dateStart"])
dateEnd   = Date(paramIn["Param"]["dateEnd"])
dateGap   = Day(paramIn["Param"]["dateGap"])
dr = collect(dateStart:dateGap:dateEnd)

#timings, tTotal = getFuncTimings(dr)
#plot(dr,tAll[:,4,3], linewidth=2, label=labels, legend=:best)
#p1 = plot(dr,tTotal, linewidth=2,xtickfontsize=16, ytickfontsize=16,
#   title="Total Test Func Timing, mesh nagfor mac",
#   xlabel="date",ylabel="execution time [s]")

#=
timings, tTotal = getFuncTimings(dr,testsuite="gfortran",verbose=verbose)

p2 = plot(dr,tTotal, linewidth=2, marker=(:dot), markersize=1,
   xtickfontsize=16, ytickfontsize=16,
   title="Total Test Func Timing, gfortran mac",
   xlabel="date",ylabel="execution time [s]")
=#
if paramIn["Param"]["timeType"] == 1

   timings, tTotal = getFuncTimings(dr,testsuite="pleiades",verbose=verbose)

   p3 = plot(dr,tTotal, linewidth=2, marker=(:dot), markersize=1,
      label="FuncTest total", xtickfontsize=16, ytickfontsize=16,
      title="Total Test Func Timing, pleiades ifort",
      xlabel="date",ylabel="execution time [s]")
end

if paramIn["Param"]["timeType"] == 2

   timings = getTimings(dr,testsuite=testsuite,testcase=testcase)

   p4 = plot(dr,timings, linewidth=2, marker=(:dot), markersize=1,
      label=testcase, xtickfontsize=16, ytickfontsize=16,
      title="Nightly test timing, "*testsuite,
      xlabel="date",ylabel="execution time [s]")
end


text:
abstract type LatticeFluidModel <: EoSModel end
abstract type SAFTModel <: LatticeFluidModel end
abstract type GCSAFTModel <: SAFTModel end
abstract type CubicModel <: EoSModel end
abstract type ABCubicModel <: CubicModel end #cubics that have an exact polynomial form to solve Z roots, this excludes CPA
abstract type ActivityModel <: EoSModel end
abstract type IdealModel <: EoSModel end
abstract type EmpiricHelmholtzModel <: EoSModel end
abstract type SatPureAproximation <: EoSModel end
abstract type AlphaModel <:EoSModel end
export SAFTModel,CubicModel,EmpiricHelmholtzModel
export IdealModel
export AlphaModel


text:
function RCI(name=:RCI)
    model = BondGraph(name)
    C = Component(:C)
    R = Component(:R)
    I = Component(:I)
    SS = Component(:SS)
    zero_law = EqualEffort()

    add_node!(model, [C, R, I, SS, zero_law])
    connect!(model, R, zero_law)
    connect!(model, C, zero_law)
    connect!(model, zero_law, I)
    connect!(model, zero_law, SS)

    model
end

# Based on https://bondgraphtools.readthedocs.io/en/latest/tutorials/RC.html
@testset "BondGraph Construction" begin
    model = BondGraph(:RC)
    C = Component(:C)
    R = Component(:R)
    zero_law = EqualEffort()

    add_node!(model, [R, C, zero_law])
    @test R in model.nodes
    @test C in model.nodes
    @test zero_law in model.nodes

    b1 = connect!(model, R, zero_law)
    b2 = connect!(model, zero_law, C)
    @test b1 in model.bonds
    @test b2 in model.bonds
end

@testset "BondGraph Modification" begin
    model = BondGraph(:RCI)
    C = Component(:C)
    R = Component(:R)
    I = Component(:I)
    SS = Component(:SS)
    zero_law = EqualEffort()
    one_law = EqualFlow()

    add_node!(model, [C, R, I, SS, zero_law, one_law])
    remove_node!(model, [SS, one_law])
    @test !(SS in model.nodes)
    @test !(one_law in model.nodes)

    connect!(model, R, zero_law)
    connect!(model, C, zero_law)

    @test I.freeports == [true]
    b1 = connect!(model, zero_law, I)
    @test ne(model) == 3
    @test b1 in model.bonds
    @test I.freeports == [false]

    disconnect!(model, zero_law, I)
    @test ne(model) == 2
    @test !(b1 in model.bonds)
    @test I.freeports == [true]

    connect!(model, zero_law, I)
    swap!(model, zero_law, one_law)
    @test I.freeports == [false]
    @test one_law in model.nodes
    @test inneighbors(model, one_law) == [R, C]
    @test outneighbors(model, one_law) == [I]
end

@testset "Construction Failure" begin
    model = BondGraph(:RC)
    C = new(:C)
    R = new(:R)
    zero_law = EqualEffort()

    add_node!(model, [R, C, zero_law])
    @test_logs (:warn, "R:R already in model") add_node!(model, R)
    @test_logs (:warn, "0 already in model") add_node!(model, zero_law)

    connect!(model, R, zero_law)
    @test_throws ErrorException connect!(model, R, zero_law)
    @test_throws ErrorException connect!(model, C, R)

    one_law = EqualFlow()
    @test_logs (:warn, "1 not in model") remove_node!(model, one_law)

    tf = new(:TF)
    add_node!(model, tf)
    @test_throws ErrorException swap!(model, tf, C)
end

@testset "Chemical reaction" begin
    model = BondGraph(:Chemical)
    A = Component(:C, :A)
    B = Component(:C, :B)
    C = Component(:C, :C)
    D = Component(:C, :D)
    Re = Component(:Re, :Reaction, numports=2)
    J_AB = EqualFlow()
    J_CD = EqualFlow()

    add_node!(model, [A, B, C, D, Re, J_AB, J_CD])
    connect!(model, A, J_AB)
    connect!(model, B, J_AB)
    connect!(model, C, J_CD)
    connect!(model, D, J_CD)

    @test freeports(Re) == [true, true]
    @test freeports(J_AB) == [false, false]

    # Connecting junctions to specific ports in Re
    connect!(model, Re, J_CD, srcportindex=2)
    @test freeports(Re) == [true, false]

    # connecting to a full port should fail
    @test_throws ErrorException connect!(model, J_AB, Re, dstportindex=2)

    connect!(model, J_AB, Re, dstportindex=1)
    @test freeports(Re) == [false, false]

    @test nv(model) == 7
    @test ne(model) == 6
end

@testset "Standard components" begin
    tf = new(:TF,:n)
    @test tf isa Component{2}
    @test tf.type == :TF
    @test numports(tf) == 2

    r = new(:R)
    @test r isa Component{1}
    @test r.type == :R
end

@parameters t
D = Differential(t)

@testset "Equations" begin
    c = new(:C)
    @parameters C
    @variables E[1](t) F[1](t) q(t)
    @test cr(c) == [
        0 ~ q/C - E[1],
        D(q) ~ F[1]
    ]
end

@testset "Parameters" begin
    tf = new(:TF)
    @parameters r
    @test iszero(BondGraphs.params(tf) - [r])

    Ce = new(:Ce)
    @parameters k R T
    @test iszero(BondGraphs.params(Ce) - [k, R, T])

    Re = new(:Re)
    @parameters r R T
    @test iszero(BondGraphs.params(Re) - [r, R, T])
end

@testset "State variables" begin
    r = new(:R)
    @test isempty(BondGraphs.state_vars(r))

    @variables q(t)
    c = new(:C)
    @test isequal(BondGraphs.state_vars(c), [q])

    ce = new(:ce)
    @test isequal(BondGraphs.state_vars(c), [q])
end

@testset "Inserting Nodes" begin
    bg = RCI()

    c, r, 𝟎 = bg.nodes[[1,2,5]]

    bondc0 = getbonds(bg, c, 𝟎)[1]
    bondr0 = getbonds(bg, r, 𝟎)[1]

    tf = Component(:TF, numports=2)
    insert_node!(bg, bondc0, tf)
    insert_node!(bg, bondr0, EqualFlow())

    @test tf in bg.nodes
    @test nv(bg) == 7
    @test ne(bg) == 6
end

@testset "Merging components" begin
    bg = RCI()

    newC = Component(:C, :newC)
    newR = Component(:R, :newR)
    add_node!(bg, [newC, newR])
    connect!(bg, newC, newR)

    C = getnodes(bg, "C")[1]
    merge_nodes!(bg, C, newC)

    R = getnodes(bg, "R")[1]
    merge_nodes!(bg, R, newR; junction=EqualFlow())

    @test isempty(getnodes(bg, "newC"))
    @test length(getnodes(bg, "1")) == 1
    @test length(getnodes(bg, "0")) == 2
    @test nv(bg) == 7
    @test ne(bg) == 7
end

@testset "Simplifying Junctions" begin
    bg = RCI()
    C, R, I, SS, 𝟎 = bg.nodes

    J0_new_1 = EqualEffort(;name=:new0_1)
    J0_new_2 = EqualEffort(;name=:new0_2)
    insert_node!(bg, (C, 𝟎), J0_new_1)
    insert_node!(bg, (R, 𝟎), J0_new_2)
    connect!(bg, J0_new_1, J0_new_2)

    J1_new_1 = EqualFlow(;name=:new1_1)
    J1_new_2 = EqualFlow(;name=:new1_2)
    add_node!(bg, J1_new_1)
    connect!(bg, 𝟎, J1_new_1)
    insert_node!(bg, (SS, 𝟎), J1_new_2)

    # Removing junction redundancies
    simplify_junctions!(bg, squash_identical=false)
    @test length(getnodes(bg, "1")) == 0
    @test nv(bg) == 7
    @test ne(bg) == 7

    # Squashing junction duplicates into a single junction
    #simplify_junctions!(bg)
    @test_broken length(getnodes(bg, :𝟎)) == 1
    @test_broken nv(bg) == 5
    @test_broken ne(bg) == 4
end

@testset "BondGraphNodes" begin
    C = Component(:C)
    bg1 = BondGraph(:BG1)
    bg2 = BondGraph(:BG2)
    bg3 = BondGraph(:BG3)
    main = BondGraph(:Main)

    bgn1 = BondGraphNode(bg1)
    bgn2 = BondGraphNode(bg2)
    bgn3 = BondGraphNode(bg3)

    @test bgn1.bondgraph === bg1
    @test bgn1.type === :BG
    @test bgn1.name === bg1.name
    @test freeports(bgn1) == Bool[]

    add_node!(bg1, C)
    add_node!(bg2, bgn1)
    add_node!(bg3, bgn2)
    add_node!(main, bgn3)

    @test main.BG3.BG2.BG1.C === C
end

# @testset "Expose component" begin

# end

# @testset "Nested BondGraphs" begin
#     bg_1 = RCI("RCI_1")
#     bg_2 = RCI("RCI_2")
#     bgn_1 = BondGraphNode(bg_1)
#     bgn_2 = BondGraphNode(bg_2)

#     bg_main = BondGraph("BG_Main")

#     #add_node!(bg_main, [bgn_1, bgn_2])
#     #connect!(bg_main, bgn_1, bgn_2)

# end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
export ConcreteDiscretePost

"""
    ConcreteDiscretePost <: DiscretePost

Textbook implementation of a discrete post operator, using concrete polyhedra
intersections.

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Notes

This operator requires that the `Polyhedra` library is loaded,
because some concrete operations between polytopes are used.

Currently, we assume that source invariants, target invariants and guards are
polytopes in constraint representation resp. half-spaces.

### Algorithm

The algorithm is based on [Flowpipe-Guard Intersection for Reachability
Computations with Support Functions](http://spaceex.imag.fr/sites/default/files/frehser_adhs2012.pdf).
"""
struct ConcreteDiscretePost <: DiscretePost
    options::Options

    function ConcreteDiscretePost(𝑂::Options)
        𝑂copy = copy(𝑂)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:check_invariant_intersection], false)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:overapproximation], Hyperrectangle)
        return new(𝑂copy)
    end
end

# convenience constructor from pairs of symbols
ConcreteDiscretePost(𝑂::Pair{Symbol,<:Any}...) = ConcreteDiscretePost(Options(Dict{Symbol,Any}(𝑂)))

# default options for the LazyDiscretePost discrete post operator
ConcreteDiscretePost() = ConcreteDiscretePost(Options())

init(𝒫::ConcreteDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝒮, copy(𝑂))

function init!(𝒫::ConcreteDiscretePost, 𝒮::AbstractSystem, 𝑂::Options)
    @assert isdefined(Main, :Polyhedra) "this algorithm needs the package " *
            "'Polyhedra' to be loaded"

    𝑂.dict[:n] = statedim(𝒮, 1)

    # solver-specific options (adds default values for unspecified options)
    𝑂out = validate_solver_options_and_add_default_values!(𝑂)

    return 𝑂out
end

function tube⋂inv!(𝒫::ConcreteDiscretePost,
                   reach_tube::Vector{<:ReachSet{<:LazySet{N}}},
                   invariant,
                   Rsets,
                   start_interval
                  ) where {N}
    # take intersection with source invariant

    # TODO First check for empty intersection, which can be more efficient.
    #      However, we need to make sure that the emptiness check does not just
    #      compute the concrete intersection; otherwise, we would do the work
    #      twice. This is currently the case for 'Polyhedra' polytopes.

    # counts the number of sets R⋂I added to Rsets
    count = 0
    for reach_set in reach_tube
        rs = reach_set.X
        if rs isa CartesianProductArray && length(array(rs)) == 1
            # TODO workaround for lazy X0
            rs_converted = Approximations.overapproximate(rs,
                Approximations.BoxDirections(dim(rs)))
        else
            rs_converted = HPolytope(constraints_list(rs))
        end
        R⋂I = intersection(invariant, rs_converted)
        if 𝒫.options[:check_invariant_intersection] && isempty(R⋂I)
            break
        end
        push!(Rsets, ReachSet{LazySet{N}}(R⋂I,
            reach_set.t_start + start_interval[1],
            reach_set.t_end + start_interval[2]))
        count = count + 1
    end

    return count
end

function post(𝒫::ConcreteDiscretePost,
              HS::HybridSystem,
              waiting_list::Vector{Tuple{Int, ReachSet{LazySet{N}}, Int}},
              passed_list,
              source_loc_id,
              tube⋂inv,
              count_Rsets,
              jumps,
              options
             ) where {N}
    jumps += 1
    for trans in out_transitions(HS, source_loc_id)
        info("Considering transition: $trans")
        target_loc_id = target(HS, trans)
        target_loc = HS.modes[target(HS, trans)]
        target_invariant = target_loc.X
        constrained_map = resetmap(HS, trans)
        guard = stateset(constrained_map)

        # perform jumps
        post_jump = Vector{ReachSet{LazySet{N}}}()
        sizehint!(post_jump, count_Rsets)
        for reach_set in tube⋂inv[length(tube⋂inv) - count_Rsets + 1 : end]
            # check intersection with guard
            R⋂G = intersection(guard, HPolytope(constraints_list(reach_set.X)))
            if isempty(R⋂G)
                continue
            end

            # apply assignment
            A⌜R⋂G⌟ = apply_assignment(𝒫, constrained_map, R⋂G)

            # intersect with target invariant
            A⌜R⋂G⌟⋂I = intersection(target_invariant, A⌜R⋂G⌟)
            if isempty(A⌜R⋂G⌟⋂I)
                continue
            end

            # store result
            push!(post_jump, ReachSet{LazySet{N}}(A⌜R⋂G⌟⋂I,
                                                     reach_set.t_start,
                                                     reach_set.t_end))
        end

        postprocess(𝒫, HS, post_jump, options, waiting_list, passed_list,
            target_loc_id, jumps)
    end
end

# --- handling assignments ---

function apply_assignment(𝒫::ConcreteDiscretePost,
                          constrained_map::ConstrainedLinearMap,
                          R⋂G::LazySet;
                          kwargs...)
    return linear_map(constrained_map.A, R⋂G)
end


text:
"""
Struct for VG Process

		vgProcess=VarianceGammaProcess(σ::num1,θ::num2,κ::num3) where {num1,num2,num3<: Number}

Where:\n
		σ =	volatility of the process.
		θ = variance of volatility.
		κ =	skewness of volatility.
"""
mutable struct VarianceGammaProcess{num <: Number, num1 <: Number, num2 <: Number, abstrUnderlying <: AbstractUnderlying, numtype <: Number} <: InfiniteActivityProcess{numtype}
    σ::num
    θ::num1
    κ::num2
    underlying::abstrUnderlying
    function VarianceGammaProcess(σ::num, θ::num1, κ::num2, underlying::abstrUnderlying) where {num <: Number, num1 <: Number, num2 <: Number, abstrUnderlying <: AbstractUnderlying}
        if σ <= 0
            error("volatility must be positive")
        elseif κ <= 0
            error("κappa must be positive")
        elseif 1 - σ * σ * κ / 2 - θ * κ < 0
            error("Parameters with unfeasible values")
        else
            zero_typed = zero(num) + zero(num1) + zero(num2)
            return new{num, num1, num2, abstrUnderlying, typeof(zero_typed)}(σ, θ, κ, underlying)
        end
    end
end

export VarianceGammaProcess;

function simulate!(X, mcProcess::VarianceGammaProcess, rfCurve::AbstractZeroRateCurve, mcBaseData::AbstractMonteCarloConfiguration, T::Number)
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    θ1 = mcProcess.θ
    κ1 = mcProcess.κ
    @assert T > 0

    dt = T / Nstep
    psi1 = -1 / κ1 * log(1 - σ^2 * κ1 / 2 - θ1 * κ1)
    drift = r - d - psi1

    gammaRandomVariable = Gamma(dt / κ1, κ1)

    simulate!(X, SubordinatedBrownianMotion(σ, drift, gammaRandomVariable), mcBaseData, T)

    f(x) = S0 * exp(x)
    broadcast!(f, X, X)

    nothing
end


text:
using Test
include("module-3-exercise-1.jl")

@testset "$(get_module_3_excersise_1_author_name()): matrix operations tests" begin
    @test get_module_3_excersise_1_author_name() != "YOUR NAME HERE"
    @test problem_1_1() == [672 2163 -1073; 134 -57 229; 50 -183 187]
    @test problem_1_2() == [-2 126 23; 176 -225 -93; 44 -60 30]
    @test problem_1_3() == [4 93 49; 100 125 -87; 70 26 36]
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using JuLIP, ACE, ZipFile, JSON

#---

species_list = [ [:Si], [:Al, :Ti] ] #, [:C, :O, :H] ]
degrees = [ 15, 14, 13, 12, 11, 12 ]
wLs = [ 1.5, 1.4, 1.3, 1.6, 1.55, 1.62 ]

zipname = dirname(pathof(ACE))[1:end-3] * "test/models/v08_compat.zip"
zipdir = ZipFile.Writer(zipname)

@info("Exporting several test basis sets and test potentials.")
for species in species_list, N = 1:length(degrees)
   local D, fptr
   @info("   ... species=$(species), N=$(N)");
   r0 = minimum(rnn.(species))
   _params = Dict("species" => species, "maxdeg" => degrees[N], "wL" => wLs[N],
                  "r0" => r0, "rcut" => 2.5 * r0, "N" => N,
                  "degreetype" => "SparsePSHDegree")
   basis = ACE.Testing.test_basis(_params)

   # file label
   flabel = "_$(N)_" * prod(string.(species)) * ".json"

   # write basis and basis tests
   D = write_dict(basis)
   D["_params"] = _params
   D["_tests"] = ACE.Testing.createtests(basis, 3; tests = ["E"])
   fptr = ZipFile.addfile(zipdir, "rpibasis" * flabel; method = ZipFile.Deflate)
   write(fptr, JSON.json(D))

   # write potential and potential tests
   V = ACE.Random.randcombine(basis; diff=1)
   D = write_dict(V)
   D["_tests"] = ACE.Testing.createtests(V, 3)
   fptr = ZipFile.addfile(zipdir, "rpipot" * flabel; method = ZipFile.Deflate)
   write(fptr, JSON.json(D))
end

close(zipdir)


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api;
#
# import typedb.client.api.database.Database;
# import typedb.protocol.SessionProto;
#
# import javax.annotation.CheckReturnValue;
#
# public interface TypeDBSession extends AutoCloseable {
#
#     @CheckReturnValue
#     boolean isOpen();
#
#     @CheckReturnValue
#     Type type();
#
#     @CheckReturnValue
#     Database database();
#
#     @CheckReturnValue
#     TypeDBOptions options();
#
#     @CheckReturnValue
#     TypeDBTransaction transaction(TypeDBTransaction.Type type);
#
#     @CheckReturnValue
#     TypeDBTransaction transaction(TypeDBTransaction.Type type, TypeDBOptions options);
#
#     void close();
#
#     enum Type {
#         DATA(0),
#         SCHEMA(1);
#
#         private final int id;
#         private final boolean isSchema;
#
#         Type(int id) {
#             this.id = id;
#             this.isSchema = id == 1;
#         }
#
#         public static Type of(int value) {
#             for (Type t : values()) {
#                 if (t.id == value) return t;
#             }
#             return null;
#         }
#
#         public int id() {
#             return id;
#         }
#
#         public boolean isData() { return !isSchema; }
#
#         public boolean isSchema() { return isSchema; }
#
#         public SessionProto.Session.Type proto() {
#             return SessionProto.Session.Type.forNumber(id);
#         }
#     }
# }


text:
module CUTENSOR

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cudaDataType
using ..CUDA: libcutensor,  @retry_reclaim

using CEnum: @cenum

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,Base.RefValue{cutensorHandle_t}}()

function handle()
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{handle::Base.RefValue{cutensorHandle_t}}
    states = get!(task_local_storage(), :CUTENSOR) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_handle = pop!(idle_handles, cuda.context) do
            handle = Ref{cutensorHandle_t}()
            cutensorInit(handle)
            handle
        end

        finalizer(current_task()) do task
            push!(idle_handles, cuda.context, new_handle) do
                # CUTENSOR doesn't need to actively destroy its handle
            end
        end

        (; handle=new_handle)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    return state.handle
end

end


text:
export solve_kepler

"""
    solve_kepler(mean_anomaly, eccentricity, desired_accuracy = 1e-6)

Solves Kepler's equation for the eccentric anomaly using Newton's method.

Arguments:
mean_anomaly    -- mean anomaly in radians
eccentricity    -- eccentricity of the ellipse defaulting to 1e-6

Returns: eccentric anomaly in radians.

M = E - e sin E

to get the time, t, or position, n, in orbit.

```
cos n = [cos E - e]/[1 - e cos E]
```

and

```
cos E = [cos n + e]/[1 + e cos n]
```

will be needed.

E = eccentric anomaly measured from perihelion about the center of the elliptical orbit
e = eccentricity of the orbit
M = mean anomaly = 2p t/P
t = time since perihelion passage
P = period of the orbit
n = true anomaly = angle of the object in obit relative to the perihelion
"""
function solve_kepler(mean_anomaly, eccentricity, desired_accuracy = 1e-6)
    eccentric_anomaly = mean_anomaly
    while true
        difference = eccentric_anomaly - eccentricity * sin(eccentric_anomaly) - mean_anomaly
        eccentric_anomaly -= difference / (1 - eccentricity * cos(eccentric_anomaly))
        if abs(difference) <= desired_accuracy
           break
        end
    end
    return eccentric_anomaly
end


text:
using ClusterComplex
using Plots
upscale = 1.0
default(size=(400*upscale,400*upscale))
pgfplots()
# gr()

X, H = ClusterComplex.dataset("TwoMoons")
s1 = scatter(X[1,:], X[2,:], 1, leg=false, markersize=1.7);
Plots.pdf(s1, "two-moons-ds")

X, H = ClusterComplex.dataset("Circles")
s2 = scatter(X[1,:], X[2,:], 1, leg=false, markersize=1.7);
Plots.pdf(s2, "circles-ds")

X, H = ClusterComplex.dataset("Sphere")
s3 = plot(X[1,:], X[2,:], X[3,:], seriestype=:scatter3d, legend=:none, markersize=1.7);
Plots.pdf(s3, "sphere-ds")

X, H = ClusterComplex.dataset("OIP300")
S = (rand(1:15000, 1200) |> unique)[1:1000]
s4 = plot(X[1,S], X[2,S], X[3,S], seriestype=:scatter3d, camera=(45, 15), legend=:none, markersize=1);
Plots.pdf(s4, "oip-300-ds")

X, H = ClusterComplex.dataset("OIP15")
S = (rand(1:15000, 1200) |> unique)[1:1000]
s5 = plot(X[1,S], X[2,S], X[3,S], seriestype=:scatter3d, camera=(45, 15), legend=:none, markersize=1);
Plots.pdf(s5, "oip-15-ds")


text:
struct Trans2CovBoolVariate{S<:AbstractFloat, T<:AbstractFloat} <: AbstractBoolVariate
    name::String
    basecovariate1::AbstractCovariate{S}
    basecovariate2::AbstractCovariate{T}
    transform::Function
end

Base.length(var::Trans2CovBoolVariate) = length(var.basecovariate1)

function slice(boolvar::Trans2CovBoolVariate{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractFloat} where {T<:AbstractFloat}
    basecov1 = boolvar.basecovariate1
    basecov2 = boolvar.basecovariate2
    f = boolvar.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip2(slice(basecov1, fromobs, toobs, slicelength), slice(basecov2, fromobs, toobs, slicelength))
    mapslice2(f, slices, slicelength, Bool)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using MPI
using P4est_wrapper
using Test

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

#############################################################################
# Definition of data structures and function callbacks
#############################################################################

# Point struct
struct point_t x::Cdouble; y::Cdouble; end

# Refine callback
function my_refine( ::Ptr{p4est_t}, which_tree::p4est_topidx_t, quadrant::Ptr{p4est_quadrant_t})
    @assert which_tree == 0
    q = quadrant[]
    ((q.level == 0) || (q.level < 5 && p4est_quadrant_child_id(quadrant) == 1)) && return Cint(1)
    return Cint(0)
end

# C-callable refine callback
const my_refine_c = @cfunction(my_refine, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}))

# Search callback function
function search_point_fn(::Ptr{p4est_t}, ::p4est_topidx_t, ::Ptr{p4est_quadrant_t}, local_num::p4est_locidx_t, point::Ptr{Cvoid})
    @info local_num, unsafe_wrap(Array, Ptr{point_t}(point), 1)[1]
    return Cint(1)
end

# C-callable search callback function
const search_point_fn_c = @cfunction(search_point_fn, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}, p4est_locidx_t, Ptr{Cvoid}))

#############################################################################
# Main program
#############################################################################

mpicomm = P4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)

# Create a connectivity structure for the unit square.
unitsquare_connectivity = p4est_connectivity_new_unitsquare()
@test unitsquare_connectivity != C_NULL

# Create a new forest
unitsquare_forest = p4est_new(mpicomm, unitsquare_connectivity, 0, C_NULL, C_NULL)
@test unitsquare_forest != C_NULL

# Register callback function to decide for refinement.
p4est_refine(unitsquare_forest, 1, my_refine_c, C_NULL)

# Create a new sc_array
ptr_to_points_sc_array = sc_array_new_count(sizeof(point_t), 2)
@test ptr_to_points_sc_array != C_NULL

# Check correct number of elements and its size in sc_array
points_sc_array = ptr_to_points_sc_array[]
@test points_sc_array.elem_count == 2 && points_sc_array.elem_size == sizeof(point_t)

# Obtain reference to first point_t
ptr_to_points_array = Ptr{point_t}(sc_array_index(ptr_to_points_sc_array,0))
@test ptr_to_points_array != C_NULL

# Overwrite point_t values in sc_array
unsafe_store!(ptr_to_points_array, point_t(0.24, 0.24), 1)
unsafe_store!(ptr_to_points_array, point_t(0.51, 0.256), 2)

# Search through quadrants, the local part of a forest, or the partition.
p4est_search(unitsquare_forest, C_NULL, search_point_fn_c, ptr_to_points_sc_array)
# Write forest to VTK file
p4est_vtk_write_file(unitsquare_forest, C_NULL, "my_step")

# Destroy sc_array structure
sc_array_destroy(ptr_to_points_sc_array)
# Destroy the forest
p4est_destroy(unitsquare_forest)
# Destroy the connectivity
p4est_connectivity_destroy(unitsquare_connectivity)

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
"""
sdm(
    holo::AbstractMatrix{T},
    windows_size::Tuple{Integer,Integer},
    interval::Tuple{Integer,Integer},
    P::FFTW.FFTWPlan,
    scale::Integer;
    shift = false
    ) -> AbstractMatrix where T <: AbstractFloat

`空域掩膜法`的算法实现，具体可见`doi: 10.1109/JDT.2015.2479646`

# Arguments
- holo: 全息图（非RGB类型而是AbstractFloat类型）
- windows_size: 掩膜窗口大小
- interval: 窗口偏移量
- holo: 原始全息图
- P: 其类型为FFTW.FFTWPlan，表示重建采用`fft`傅立叶变换
- scale: 由于零级像强度过大，故需要用`scale`拉升±1级像的强度
- shift: 是否旋转图像，默认为否（除非是想展示完整的重建像，否则不建议将`shift`设为`true`，因为其内存和时间开销比较大）
"""
function sdm(holo::AbstractMatrix{T}, windows_size::Tuple{Integer,Integer},
    interval::Tuple{Integer,Integer}, P::FFTW.FFTWPlan, scale::Integer; shift = false) where T <: AbstractFloat

    # 计算子全息图的个数并进行预分配内存
    nx, ny = (size(holo) .- windows_size) .÷ interval
    S = (nx + 1) * (ny + 1)
    holo_tensor = zeros(T, size(holo, 1), size(holo, 2), S)

    # 获取子全息图
    make_sub_holo!(holo_tensor, holo, windows_size, interval, (nx, ny))

    if typeof(P) <: FFTW.rFFTWPlan
        re_tensor = similar(holo, size(holo, 1) ÷ 2 + 1, size(holo, 2), S)
        sdm_img = similar(holo, size(holo, 1) ÷ 2 + 1, size(holo, 2))
    elseif typeof(P) <: FFTW.cFFTWPlan
        re_tensor = similar(holo, size(holo, 1), size(holo, 2), S)
        sdm_img = similar(holo, size(holo, 1), size(holo, 2))
    else
        @error "P的类型不正确..."
    end

    # 对子全息图进行重建
    reconst_tensor!(re_tensor, holo_tensor, P, scale; shift = shift)

    # 求其均值
    sdm_core!(sdm_img, re_tensor)

    return sdm_img

end

"""
    sdm_core!(sdm_img::AbstractMatrix{T}, re_tensor::AbstractArray{T,3}) -> AbstractMatrix where T <: AbstractFloat

`空域掩膜法`的核心算法，对子再现像序列进行叠加平均，并存进`sdm_img`中

# Arguments
- sdm_img: 存放结果的矩阵
- re_tensor: 子再现像序列
"""
function sdm_core!(sdm_img::AbstractMatrix{T}, re_tensor::AbstractArray{T,3}) where T <: AbstractFloat
    dims1, dims2, _ = size(re_tensor)
    @inbounds Threads.@threads for ind in CartesianIndices((1:dims1, 1:dims2))
        i, j = ind.I
        sdm_img[i, j] = mean(@view re_tensor[i, j, :])
    end

    return sdm_img
end


text:
using InfiniteOptSymbolics
using Test

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


text:
using Test
# using Revise
using CutCellDG
include("../useful_routines.jl")

l,m = 1.,2.
stiffness = CutCellDG.plane_strain_voigt_hooke_matrix(l,m)
teststiffness = [5. 1. 0.
                 1. 5. 0.
                 0. 0. 2.]
@test allapprox(teststiffness,stiffness)

stiffness = CutCellDG.HookeStiffness(l,m,l,m)
@test allapprox(teststiffness,stiffness[-1])
@test allapprox(teststiffness,stiffness[+1])

l1,m1 = 1.,2.
l2,m2 = 3.,4.
stiffness = CutCellDG.HookeStiffness(l1,m1,l2,m2)
teststiffness2 = [11.  3.  0.
                   3.  11. 0.
                   0.   0. 4.]
teststiffness1 = [5. 1. 0.
                  1. 5. 0.
                  0. 0. 2.]
@test allapprox(teststiffness1,stiffness[1])
@test allapprox(teststiffness2,stiffness[-1])

@test allapprox((l1,m1),CutCellDG.lame_coefficients(stiffness,1))
@test allapprox((l2,m2),CutCellDG.lame_coefficients(stiffness,-1))


text:
include("Chipmunk.jl")

space = Chipmunk.cpSpaceNew()
Chipmunk.cpSpaceSetGravity(space, Chipmunk.cpVect(0, 10))
body = Chipmunk.cpBodyNew(1, 1)
Chipmunk.cpSpaceAddBody(space, body)

for i = 1:10
    Chipmunk.cpSpaceStep(space, 1 / 60)

    p = Chipmunk.cpBodyGetPosition(body)
    v = Chipmunk.cpBodyGetVelocity(body)

    println("position = ($(p.x), $(p.y)), velocity = ($(v.x), $(v.y))")
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SLIC(k, m; tol=1e-4, maxiter=10)

A method for clustering geospatial data into approximately `k`
clusters using Simple Linear Iterative Clustering (SLIC).
The method produces clusters of samples that are spatially
connected based on a distance `dₛ` and that, at the same
time, are similar in terms of `vars` with distance `dᵥ`.
The tradeoff is controlled with a hyperparameter parameter
`m` in an additive model `dₜ = √(dᵥ² + m²(dₛ/s)²)`.

## Parameters

* `k`       - Approximate number of clusters
* `m`       - Hyperparameter of SLIC model
* `tol`     - Tolerance of k-means algorithm (default to `1e-4`)
* `maxiter` - Maximum number of iterations (default to `10`)
* `weights` - Dictionary with weights for each attribute (default to `nothing`)

## References

* Achanta et al. 2011. [SLIC superpixels compared to state-of-the-art
  superpixel methods](https://ieeexplore.ieee.org/document/6205760)
"""
struct SLIC{W} <: ClusteringMethod
  k::Int
  m::Float64
  tol::Float64
  maxiter::Int
  weights::W
end

function SLIC(k::Int, m::Real; tol=1e-4, maxiter=10, weights=nothing)
  @assert tol > 0 "invalid tolerance"
  @assert maxiter > 0 "invalid number of iterations"
  SLIC{typeof(weights)}(k, m, tol, maxiter, weights)
end

function partition(data, method::SLIC)
  # normalize atributes
  𝒯 = TableDistances.normalize(values(data))
  Ω = georef(first(𝒯), domain(data))

  # weights for each attribute
  weights = method.weights

  # SLIC hyperparameter
  m = method.m

  # initial spacing of clusters
  s = slic_spacing(Ω, method)

  # initialize cluster centers
  c = slic_initialization(Ω, s)

  # ball neighborhood search
  searcher = BallSearch(Ω, NormBall(s))

  # pre-allocate memory for label and distance
  l = fill(0, nelements(Ω))
  d = fill(Inf, nelements(Ω))

  # performance parameters
  tol     = method.tol
  maxiter = method.maxiter

  # k-means algorithm
  err, iter = Inf, 0
  while err > tol && iter < maxiter
    o = copy(c)

    slic_assignment!(Ω, searcher, weights, m, s, c, l, d)
    slic_update!(Ω, c, l)

    err = norm(c - o) / norm(o)
    iter += 1
  end

  orphans = findall(iszero, l)
  if length(orphans) > 0
    assigned = findall(!iszero, l)
    Ω₀ = view(domain(Ω), assigned)
    csearcher = KNearestSearch(Ω₀, 1)

    for orphan in orphans
      p = centroid(Ω, orphan)
      i = search(p, csearcher)[1]
      l[orphan] = l[assigned[i]]
    end
  end

  subsets = [findall(isequal(k), l) for k in 1:length(c)]

  Partition(data, subsets)
end

function slic_spacing(data, method)
  V = measure(boundingbox(data))
  d = embeddim(data)
  k = method.k
  (V/k) ^ (1/d)
end

function slic_initialization(data, s)
  # efficient neighbor search
  searcher = KNearestSearch(data, 1)

  # bounding box properties
  bbox = boundingbox(data)
  lo, up = coordinates.(extrema(bbox))

  # cluster centers
  clusters = Vector{Int}()
  neighbor = Vector{Int}(undef, 1)
  ranges = [(l+s/2):s:u for (l, u) in zip(lo, up)]
  for x in Iterators.product(ranges...)
    search!(neighbor, Point(x), searcher)
    push!(clusters, neighbor[1])
  end

  unique(clusters)
end

function slic_assignment!(data, searcher, weights, m, s, c, l, d)
  for (k, cₖ) in enumerate(c)
    pₖ = centroid(data, cₖ)
    inds = search(pₖ, searcher)

    # distance between points
    X  = (coordinates(centroid(data, ind)) for ind in inds)
    xₖ = [coordinates(pₖ)]
    dₛ = pairwise(Euclidean(), X, xₖ)

    # distance between variables
    𝒮ᵢ = view(data, inds)
    𝒮ₖ = view(data, [cₖ])
    V  = values(𝒮ᵢ)
    vₖ = values(𝒮ₖ)
    td = TableDistance(normalize=false, weights=weights)
    dᵥ = pairwise(td, V, vₖ)

    # total distance
    dₜ = @. √(dᵥ^2 + m^2 * (dₛ/s)^2)

    @inbounds for (i, ind) in enumerate(inds)
      if dₜ[i] < d[ind]
        d[ind] = dₜ[i]
        l[ind] = k
      end
    end
  end
end

function slic_update!(data, c, l)
  for k in 1:length(c)
    inds = findall(isequal(k), l)
    X  = (coordinates(centroid(data, ind)) for ind in inds)
    μ  = [mean(X)]
    dₛ = pairwise(Euclidean(), X, μ)
    @inbounds c[k] = inds[argmin(vec(dₛ))]
  end
end


text:
export sSmsPolicy

struct sSPolicy end
function (::sSPolicy)(item::AbstractItem, s::T, S::T) where T<:Number
    ip = inventory_position(item)
    max(0,(ip <= s)*(S - ip))
end
(p::sSPolicy)(item::AbstractItem, v::AbstractVector) = p(item, v...)

action_size(::sSPolicy) = 2
print_action(::sSPolicy) = ["s", "S"]

struct RQPolicy end

function (::RQPolicy)(item::AbstractItem, R, Q)
    max(0,(inventory_position(item) <= R)*Q)
end
(p::RQPolicy)(item::AbstractItem, v::AbstractVector) = p(item, v...)

action_size(::RQPolicy) = 2
print_action(::RQPolicy) = ["R", "Q"]

struct QPolicy end

function (::QPolicy)(::AbstractItem, Q)
    max(0,Q)
end
(p::QPolicy)(item::AbstractItem, v::AbstractVector) = p(item, v...)

action_size(::QPolicy) = 1
print_action(::QPolicy) = ["Q"]

struct YQPolicy end

function (::YQPolicy)(::AbstractItem, Y, Q)
    max(0,(Y > 0)*Q)
end
(p::YQPolicy)(item::AbstractItem, v::AbstractVector) = p(item, v...)
action_size(::YQPolicy) = 2
print_action(::YQPolicy) = ["Y", "Q"]

struct sSmsPolicy end
function (::sSmsPolicy)(item::AbstractItem, s::T, Sms::T) where T<:Number
    ip = inventory_position(item)
    max(0,(ip <= s)*(Sms+s - ip))
end
(p::sSmsPolicy)(item::AbstractItem, v::AbstractVector) = p(item, v...)

action_size(::sSmsPolicy) = 2
print_action(::sSmsPolicy) = ["s", "S-s"]

#(R, nQ)
#StaticStaticPolicy
#StaticDynamicPolicy

max_stars_repo_path: src/elementary blocks/policies.jl

max_issues_repo_path: src/elementary blocks/policies.jl

max_forks_repo_path: src/elementary blocks/policies.jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEMDocs.jl/blob/master/LICENSE

module JuliaFEMDocs

"""
    copydocs(package)

Copy documentation of some package to the src dir of JuliaFEMDocs/docs.
"""
function copydocs(package)
    # source
    pkg_path = abspath(pathof(package), "..", "..")
    src_docs = joinpath(pkg_path, "docs", "src")
    src_readme = joinpath(pkg_path, "README.md")
    # destination
    dst_root = abspath(pathof(JuliaFEMDocs), "..", "..", "docs", "src", "packages")
    dst_path = joinpath(dst_root, string(package))
    dst_readme = joinpath(dst_path, "README.md")
    # copy
    isdir(dst_root) || mkdir(dst_root)
    cp(src_docs, dst_path)
    cp(src_readme, dst_readme)
    return nothing
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsReaction(;
        count=nothing,
        name=nothing,
        users=nothing,
    )

    - count::Int32
    - name::String
    - users::Vector{DefsUserId}
"""
mutable struct ObjsReaction <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    users::Any # spec type: Union{ Nothing, Vector{DefsUserId} } # spec name: users

    function ObjsReaction(;count=nothing, name=nothing, users=nothing)
        o = new()
        validate_property(ObjsReaction, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(ObjsReaction, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ObjsReaction, Symbol("users"), users)
        setfield!(o, Symbol("users"), users)
        o
    end
end # type ObjsReaction

const _property_map_ObjsReaction = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"), Symbol("name")=>Symbol("name"), Symbol("users")=>Symbol("users"))
const _property_types_ObjsReaction = Dict{Symbol,String}(Symbol("count")=>"Int32", Symbol("name")=>"String", Symbol("users")=>"Vector{DefsUserId}")
Base.propertynames(::Type{ ObjsReaction }) = collect(keys(_property_map_ObjsReaction))
Swagger.property_type(::Type{ ObjsReaction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsReaction[name]))}
Swagger.field_name(::Type{ ObjsReaction }, property_name::Symbol) =  _property_map_ObjsReaction[property_name]

function check_required(o::ObjsReaction)
    (getproperty(o, Symbol("count")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("users")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsReaction }, name::Symbol, val)
end


text:
ops = (Rotate180(),Crop(5:200,200:500),Rotate90(1),Crop(1:250, 1:150))
@testset "$(str_showcompact(ops))" begin
    @test_throws ArgumentError augmentbatch!(similar(camera, 250, 150, 3), cameras, ops)
    @test_throws DimensionMismatch augmentbatch!(similar(camera, 250, 140, 2), cameras, ops)
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(out, cameras, ops)) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(out, cameras, AugmentorBase.ImmutablePipeline(ops))) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(out, collect.(obsview(cameras)), ops)) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    outs = [similar(camera, 250, 150), similar(camera, 250, 150)]
    @test @inferred(augmentbatch!(outs, cameras, ops)) === outs
    @test typeof(outs) <: Vector
    @test eltype(outs) <: Array{eltype(camera)}
    @test_reference "reference/rot_crop_either_crop.txt" outs[1]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2]
    out = similar(camera, 2, 250, 150)
    cameras_t = permutedims(cameras, (3,1,2))
    @test @inferred(augmentbatch!(out, cameras_t, ops, ObsDim.First())) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[1,:,:]
    @test_reference "reference/rot_crop_either_crop.txt" out[2,:,:]
end
@testset "$(str_showcompact(ops)) multiple images" begin
    outs = (similar(N0f8.(camera), 250, 150, 2), similar(camera, 250, 150, 2))
    @test @inferred(augmentbatch!(outs, (N0f8.(cameras), cameras), ops)) === outs
    @test typeof(outs) <: Tuple{Array{N0f8,3},Array{Gray{N0f8},3}}
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,1])
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,2])
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,2]
    @test @inferred(augmentbatch!(outs, obsview((N0f8.(cameras), cameras)), ops)) === outs
    @test typeof(outs) <: Tuple{Array{N0f8,3},Array{Gray{N0f8},3}}
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,1])
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,2])
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,2]
    outs = (similar(N0f8.(camera), 2, 250, 150), similar(camera, 2, 250, 150))
    cameras_t = permutedims(cameras, (3,1,2))
    @test @inferred(augmentbatch!(outs, (N0f8.(cameras_t), cameras_t), ops, ObsDim.First())) === outs
    @test typeof(outs) <: Tuple{Array{N0f8,3},Array{Gray{N0f8},3}}
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][1,:,:])
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][2,:,:])
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][1,:,:]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][2,:,:]
end

@testset "Multithreaded: $(str_showcompact(ops))" begin
    @test_throws ArgumentError augmentbatch!(CPUThreads(), similar(camera, 250, 150, 3), cameras, ops)
    # doesn't work because exception is thrown in thread
    # @test_throws DimensionMismatch augmentbatch!(CPUThreads(), similar(camera, 250, 140, 2), cameras, ops)
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(CPUThreads(), out, cameras, ops)) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(CPUThreads(), out, cameras, AugmentorBase.ImmutablePipeline(ops))) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    out = similar(camera, 250, 150, 2)
    @test @inferred(augmentbatch!(CPUThreads(), out, collect.(obsview(cameras)), ops)) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" out[:,:,2]
    outs = [similar(camera, 250, 150), similar(camera, 250, 150)]
    @test @inferred(augmentbatch!(CPUThreads(), outs, cameras, ops)) === outs
    @test typeof(outs) <: Vector
    @test eltype(outs) <: Array{eltype(camera)}
    @test_reference "reference/rot_crop_either_crop.txt" outs[1]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2]
    out = similar(camera, 2, 250, 150)
    cameras_t = permutedims(cameras, (3,1,2))
    @test @inferred(augmentbatch!(CPUThreads(), out, cameras_t, ops, ObsDim.First())) === out
    @test typeof(out) <: Array
    @test eltype(out) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" out[1,:,:]
    @test_reference "reference/rot_crop_either_crop.txt" out[2,:,:]
end
@testset "Multithreaded: $(str_showcompact(ops)) multiple images" begin
    outs = (similar(N0f8.(camera), 250, 150, 2), similar(camera, 250, 150, 2))
    @test @inferred(augmentbatch!(CPUThreads(), outs, (N0f8.(cameras), cameras), ops)) === outs
    @test typeof(outs) <: Tuple{Array{N0f8,3},Array{Gray{N0f8},3}}
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,1])
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][:,:,2])
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,1]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][:,:,2]
    outs = (similar(N0f8.(camera), 2, 250, 150), similar(camera, 2, 250, 150))
    cameras_t = permutedims(cameras, (3,1,2))
    @test @inferred(augmentbatch!(CPUThreads(), outs, (N0f8.(cameras_t), cameras_t), ops, ObsDim.First())) === outs
    @test typeof(outs) <: Tuple{Array{N0f8,3},Array{Gray{N0f8},3}}
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][1,:,:])
    @test_reference "reference/rot_crop_either_crop.txt" Gray.(outs[1][2,:,:])
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][1,:,:]
    @test_reference "reference/rot_crop_either_crop.txt" outs[2][2,:,:]
end

ops = Rotate90()
@testset "$(str_showcompact(ops))" begin
    out = similar(camera, 512, 512, 2)
    @test @inferred(augmentbatch!(out, cameras, ops)) === out
    outs = (similar(camera, 512, 512, 2), similar(camera, 512, 512, 2))
    @test @inferred(augmentbatch!(outs, (N0f8.(cameras), cameras), ops)) === outs
end


text:
using Documenter
using SubstitutionTilings

makedocs(
    sitename = "SubstitutionTilings",
    format = Documenter.HTML(prettyurls = false),
    modules = [SubstitutionTilings],
)

deploydocs(
    repo = "github.com/danrocag/SubstitutionTilings.jl.git"
)


text:
"""
Bridges `CP.AllEqual` to a series of `CP.EqualTo`.
"""
struct AllEqual2EqualToBridge{T} <: MOIBC.AbstractBridge
    # For an AllEqual set of dimension `d`, vector of size `d-1`: first
    # variable equal to the `i+1`th.
    cons::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{AllEqual2EqualToBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.AllEqual,
) where {T}
    return MOIBC.bridge_constraint(
        AllEqual2EqualToBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{AllEqual2EqualToBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.AllEqual,
) where {T}
    f_scalars = MOIU.scalarize(f)
    dim = MOI.output_dimension(f)

    cons = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_scalars[1] - f_scalars[i],
            MOI.EqualTo(zero(T))
        )
        for i in 2:dim
    ]

    return AllEqual2EqualToBridge(cons)
end

function MOI.supports_constraint(
    ::Type{AllEqual2EqualToBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.AllEqual},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{AllEqual2EqualToBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{AllEqual2EqualToBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(
    b::AllEqual2EqualToBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons)
end

function MOI.get(
    b::AllEqual2EqualToBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return collect(values(b.cons))
end


text:
function white_box_jsr(A,d=2)
    s = discreteswitchedsystem(A)
    optimizer_constructor = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
    soslyapb(s, d, optimizer_constructor=optimizer_constructor, tol=1e-4, verbose=0)
    seq = sosbuildsequence(s, d, p_0=:Primal)
    psw = findsmp(seq)
    return psw.growthrate
end
#=
function white_box_stabilization_quad(A,B,tol=1e-4)
    normAmax = 0
    dim,dimIn = size(B)
    for Ai in A
        normAi = opnorm(Ai)
        if normAi>normAmax
            normAmax = normAi
        end
    end
    ub = normAmax
    lb = 0
    K = zeros(dimIn,dim)
    P = Matrix(I,dim,dim)
    while ub-lb > tol
        gamma = (ub + lb)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, Q[1:dim, 1:dim] in PSDCone())
        @variable(model, Y[1:dimIn,1:dim])
        @variable(model, s>=0)
        @SDconstraint(model, Q >= Matrix(I,dim,dim))
        @objective(model, Min, s)
        for Ai in A
            @SDconstraint(model, [gamma^2*Q+s*Matrix(I,dim,dim) Q*Ai'+(B*Y)';Ai*Q+B*Y Q+s*Matrix(I,dim,dim)] >= 0)
        end
        JuMP.optimize!(model)
        if value.(s) < 1e-10
            ub = gamma
            P = inv(value.(Q))
            K = value.(Y)*P
        else
            lb = gamma
        end
    end
    return ub, K, P
end
=#

function white_box_stabilization_quad(A,B,tol=1e-4)
    normAmax = 0
    dim,dimIn = size(B)
    for Ai in A
        normAi = opnorm(Ai)
        if normAi>normAmax
            normAmax = normAi
        end
    end
    ub = normAmax
    lb = 0
    K = zeros(dimIn,dim)
    P = Matrix(I,dim,dim)
    while ub-lb > tol
        gamma = (ub + lb)/2
        solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
        model = Model(solver)
        @variable(model, Q[1:dim, 1:dim] in PSDCone())
        @variable(model, Y[1:dimIn,1:dim])
        @SDconstraint(model, Q >= Matrix(I,dim,dim))
        @objective(model, Max, 0)
        for Ai in A
            @SDconstraint(model, [gamma^2*Q Q*Ai'+(B*Y)';Ai*Q+B*Y Q] >= 0)
        end
        JuMP.optimize!(model)
        if termination_status(model) == MOI.OPTIMAL
            ub = gamma
            P = inv(value.(Q))
            K = value.(Y)*P
        else
            lb = gamma
        end
    end
    return ub, K, P
end

function white_box_LQR(A,B,Q,R)
    (dim,dimIn) = size(B)
    lower_triangular(P) = [P[i, j] for i = 1:size(P, 1) for j = 1:i]
    solver = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)
    model = Model(solver)
    @variable(model, S[1:dim, 1:dim] in PSDCone())
    @variable(model, Y[1:dimIn,1:dim])
    @variable(model, t)
    @SDconstraint(model, S >= 0)
    @constraint(model, [t; 1; lower_triangular(S)] in MOI.LogDetConeTriangle(dim))
    @objective(model, Max, t)
    for Ai in A
        @SDconstraint(model, [S S*Ai'+(B*Y)' S Y';Ai*S+B*Y S zeros(dim,dim) zeros(dim,dimIn);S zeros(dim,dim) inv(Q) zeros(dim,dimIn);Y zeros(dimIn,dim) zeros(dimIn,dim) inv(R)] >= 0)
    end
    JuMP.optimize!(model)
    if termination_status(model) == MOI.OPTIMAL
        P = inv(value.(S))
        K = value.(Y)*P
        return K, P
    else
        println("The LQR problem is infeasible!")
        return zeros(dimIn,dim), zeros(dim,dim)
    end
end

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
module MultiFieldCellFieldsTests

using BlockArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using Gridap.CellData
using Gridap.MultiField
using Gridap.TensorValues
using Test

domain = (0,1,0,1)
cells = (2,2)
model = CartesianDiscreteModel(domain,cells)

trian = Triangulation(model)

u1(x) = sin(x[1])
cf1 = CellField(u1,trian)

u2(x) = cos(x[2])
cf2 = CellField(u2,trian)

cf = MultiFieldCellField([cf1,cf2])

@test cf1 === cf[1]
@test cf2 === cf[2]

_cf1, _cf2 = cf

@test cf1 === _cf1
@test cf2 === _cf2

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

degree = order

trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)
q = get_coordinates(quad)

trian_Γ = SkeletonTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,degree)
q_Γ = get_coordinates(quad_Γ)

V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=VectorValue{2,Float64})
Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

dv, dq = get_cell_basis(Y)
du, dp = get_cell_basis(X)

n = VectorValue(1,2)

cellmat = integrate( (n⋅dv)*dp + dq*dp, trian, quad)
cellvec = integrate( n⋅dv, trian, quad)
@test isa(cellmat,VectorOfBlockArrayCoo)
@test isa(cellvec,VectorOfBlockArrayCoo)

cellmat1 = integrate( ((n⋅dv) - dq)*((n⋅du) + dp), trian, quad)
cellmat2 = integrate( (n⋅dv)*(n⋅du) + (n⋅dv)*dp - dq*(n⋅du) - dq*dp, trian, quad)
test_array(cellmat1,cellmat2,≈)

cellmat1 = integrate( (n⋅dv)*2, trian, quad)
cellmat2 = integrate( (n⋅dv)*fill(2,num_cells(trian)), trian, quad)
test_array(cellmat1,cellmat2,≈)

α = CellField(2,trian)
op(u,∇u,v,∇v,α) = α*(u⋅v) + ∇u⊙∇v
cellmat1 = integrate( operate(op,du,∇(du),dv,∇(dv),α) , trian, quad)
cellmat2 = integrate( α*(du⋅dv) + ∇(du)⊙∇(dv) , trian, quad)
test_array(cellmat1,cellmat2,≈)
@test isa(cellmat2,VectorOfBlockArrayCoo)
@test isa(cellmat1,VectorOfBlockArrayCoo)

α = CellField(2,trian)
op2(u,∇u,α) = α*(∇u⋅u)
cellmat1 = integrate( dv⋅operate(op2,du,∇(du),α) , trian, quad)
cellmat2 = integrate( dv⋅(α*(∇(du)⋅du)), trian, quad)
test_array(cellmat1,cellmat2,≈)
@test isa(cellmat2,VectorOfBlockArrayCoo)
@test isa(cellmat1,VectorOfBlockArrayCoo)

conv(u,∇u,α) = α*(u⋅∇u)
dconv(du,∇du,u,∇u,α) = conv(u,∇du,α)+conv(du,∇u,α)

u = zero(U)
cellvec2 = integrate(dv⊙(α*(u⋅∇(u))),trian,quad)
cellvec1 = integrate(dv⊙operate(conv,u,∇(u),α),trian,quad)
test_array(cellvec1,cellvec2,≈)
@test isa(cellvec2,VectorOfBlockArrayCoo)
@test isa(cellvec1,VectorOfBlockArrayCoo)

cellmat1 = integrate( dv⋅operate(dconv,du,∇(du),u,∇(u),α) , trian, quad)
cellmat2 = integrate( dv⋅( α*(du⋅∇(u)) + α*(u⋅∇(du))), trian, quad)
test_array(cellmat1,cellmat2,≈)
@test isa(cellmat2,VectorOfBlockArrayCoo)
@test isa(cellmat1,VectorOfBlockArrayCoo)

dv_Γ, dq_Γ = restrict(get_cell_basis(Y), trian_Γ)
du_Γ, dp_Γ = restrict(get_cell_basis(X), trian_Γ)

cellmat_Γ = integrate(  jump(n⋅dv_Γ)*dp_Γ.⁺ + mean(dq_Γ)*jump(dp_Γ), trian_Γ,quad_Γ)
cellvec_Γ = integrate(  jump(n⋅dv_Γ) + mean(dq_Γ), trian_Γ,quad_Γ)
L = 1
R = 2
@test isa(cellmat_Γ,VectorOfBlockArrayCoo)
@test isa(cellmat_Γ[Block(L,R)],VectorOfBlockArrayCoo)
@test isa(cellvec_Γ,VectorOfBlockArrayCoo)
@test isa(cellvec_Γ[Block(L)],VectorOfBlockArrayCoo)

cell = 1
@test isa(cellmat_Γ[cell][Block(L,R)],BlockArrayCoo)
@test isa(cellvec_Γ[cell][Block(L)],BlockArrayCoo)

cellmat1_Γ = integrate(((n⋅dv_Γ.⁺)-dq_Γ.⁻)*((n⋅du_Γ.⁺)+dp_Γ.⁻),trian_Γ,quad_Γ)
cellmat2_Γ = integrate((n⋅dv_Γ.⁺)*(n⋅du_Γ.⁺)+(n⋅dv_Γ.⁺)*dp_Γ.⁻-dq_Γ.⁻*(n⋅du_Γ.⁺)-dq_Γ.⁻*dp_Γ.⁻,trian_Γ,quad_Γ)
test_array(cellmat1_Γ,cellmat2_Γ,≈)

#a = cellmat_Γ
#using BenchmarkTools
#cache = array_cache(a)
#@btime getindex!($cache,$a,2)

end # module


text:
"""
    AnchoredInterval{P,T,L,R}

`AnchoredInterval` is a subtype of `AbstractInterval` that represents a non-iterable range
or span of values defined not by two endpoints but instead by a single `anchor` point and
the value type `P` which represents the size of the range. When `P` is positive, the
`anchor` represents the lesser endpoint (the beginning of the range); when `P` is negative,
the `anchor` represents the greater endpoint (the end of the range).

The interval represented by an `AnchoredInterval` value may be closed (both endpoints are
included in the interval), open (neither endpoint is included), or half-open. This openness
is defined by the bounds types `L` and `R`, which defaults to half-open (with the lesser
endpoint included for positive values of `P` and the greater endpoint included for negative
values).

### Why?

`AnchoredIntervals` are most useful in cases where a single value is used to stand in for a
range of values. This happens most often with dates and times, where "HE15" is often used as
shorthand for (14:00..15:00].

To this end, `HourEnding` is a type alias for `AnchoredInterval{Hour(-1)}`. Similarly,
`HourBeginning` is a type alias for `AnchoredInterval{Hour(1)}`.

### Rounding

While the user may expect an `HourEnding` or `HourBeginning` value to be anchored to a
specific hour, the constructor makes no guarantees that the anchor provided is rounded:

```jldoctest; setup = :(using Intervals, Dates)
julia> HourEnding(DateTime(2016, 8, 11, 2, 30))
AnchoredInterval{-1 hour,DateTime,Open,Closed}(2016-08-11T02:30:00)
```

The `HE` and `HB` pseudoconstructors round the input up or down to the nearest hour, as
appropriate:

```jldoctest; setup = :(using Intervals, Dates)
julia> HE(DateTime(2016, 8, 11, 2, 30))
AnchoredInterval{-1 hour,DateTime,Open,Closed}(2016-08-11T03:00:00)

julia> HB(DateTime(2016, 8, 11, 2, 30))
AnchoredInterval{1 hour,DateTime,Closed,Open}(2016-08-11T02:00:00)
```

### Example

```jldoctest; setup = :(using Intervals, Dates)
julia> AnchoredInterval{Hour(-1)}(DateTime(2016, 8, 11, 12))
AnchoredInterval{-1 hour,DateTime,Open,Closed}(2016-08-11T12:00:00)

julia> AnchoredInterval{Day(1)}(DateTime(2016, 8, 11))
AnchoredInterval{1 day,DateTime,Closed,Open}(2016-08-11T00:00:00)

julia> AnchoredInterval{Minute(5),Closed,Closed}(DateTime(2016, 8, 11, 12, 30))
AnchoredInterval{5 minutes,DateTime,Closed,Closed}(2016-08-11T12:30:00)
```

See also: [`Interval`](@ref), [`HE`](@ref), [`HB`](@ref)
"""
struct AnchoredInterval{P, T, L <: Bounded, R <: Bounded} <: AbstractInterval{T,L,R}
    anchor::T

    function AnchoredInterval{P,T,L,R}(anchor::T) where {P, T, L <: Bounded, R <: Bounded}
        # A valid interval requires that neither endpoints or the span are nan. Typically,
        # we use `left <= right` to ensure a valid interval but for `AnchoredInterval`s
        # computing the other endpoint requires `anchor + P` which may fail with certain
        # types (e.g. ambiguous or non-existent ZonedDateTimes).
        #
        # We can skip computing the other endpoint if both the anchor and span are finite as
        # this ensures the computed endpoint is also finite.
        if !_isfinite(anchor) || !_isfinite(P)
            left, right = sign(P) < 0 ? (anchor + P, anchor) : (anchor, anchor + P)

            if !(left <= right)
                msg = if sign(P) < 0
                    "Unable to represent a right-anchored interval where the " *
                    "left ($anchor + $P) > right ($anchor)"
                else
                    "Unable to represent a left-anchored interval where the " *
                    "left ($anchor) > right ($anchor + $P)"
                end
                throw(ArgumentError(msg))
            end
        end

        return new{P,T,L,R}(anchor)
    end
end

_isfinite(x) = iszero(x - x)
_isfinite(x::Real) = Base.isfinite(x)

function AnchoredInterval{P,T,L,R}(anchor) where {P,T,L,R}
    AnchoredInterval{P,T,L,R}(convert(T, anchor))
end

AnchoredInterval{P,L,R}(anchor::T) where {P,T,L,R} = AnchoredInterval{P,T,L,R}(anchor)

# When an interval is anchored to the lesser endpoint, default to Inclusivity(false, true)
# When an interval is anchored to the greater endpoint, default to Inclusivity(true, false)
function AnchoredInterval{P,T}(anchor) where {P,T}
    L = bound_type(P ≥ zero(P))
    R = bound_type(P ≤ zero(P))
    return AnchoredInterval{P,T,L,R}(anchor)
end

AnchoredInterval{P}(anchor::T) where {P,T} = AnchoredInterval{P,T}(anchor)

# Note: Ideally we would define the restriction `T <: TimeType` but doing so interferes with
# the `HourEnding{L,R}` constructor.
"""
    HourEnding{T<:TimeType, L, R} <: AbstractInterval{T}

A type alias for `AnchoredInterval{Hour(-1), T}` which is used to denote a 1-hour period of
time which ends at a time instant (of type `T`).
"""
const HourEnding{T,L,R} = AnchoredInterval{Hour(-1), T, L, R} where {T, L <: Bounded, R <: Bounded}
HourEnding(anchor::T) where T = HourEnding{T}(anchor)

# Note: Ideally we would define the restriction `T <: TimeType` but doing so interferes with
# the `HourBeginning{L,R}` constructor.
"""
    HourBeginning{T<:TimeType, L, R} <: AbstractInterval{T}

A type alias for `AnchoredInterval{Hour(1), T}` which is used to denote a 1-hour period of
time which begins at a time instant (of type `T`).
"""
const HourBeginning{T,L,R} = AnchoredInterval{Hour(1), T, L, R} where {T, L <: Bounded, R <: Bounded}
HourBeginning(anchor::T) where T = HourBeginning{T}(anchor)

"""
    HE(anchor) -> HourEnding

`HE` is a pseudoconstructor for [`HourEnding`](@ref) that rounds the anchor provided up to the
nearest hour.
"""
HE(anchor) = HourEnding(ceil(anchor, Hour))

"""
    HB(anchor) -> HourBeginning

`HB` is a pseudoconstructor for [`HourBeginning`](@ref) that rounds the anchor provided down to the
nearest hour.
"""
HB(anchor) = HourBeginning(floor(anchor, Hour))

function Base.copy(x::AnchoredInterval{P,T,L,R}) where {P,T,L,R}
    return AnchoredInterval{P,T,L,R}(anchor(x))
end

##### ACCESSORS #####

# We would typically compute `first` and `last` using `min` and `max` respectively, but we
# can get unexpected behaviour if adding the span to the anchor endpoint produces a value
# that is no longer comparable (e.g., `NaN`).

function Base.first(interval::AnchoredInterval{P}) where P
    P < zero(P) ? (interval.anchor + P) : (interval.anchor)
end

function Base.last(interval::AnchoredInterval{P}) where P
    P < zero(P) ? (interval.anchor) : (interval.anchor + P)
end

anchor(interval::AnchoredInterval) = interval.anchor
span(interval::AnchoredInterval{P}) where P = abs(P)

##### CONVERSION #####

# Allows an interval to be converted to a scalar when the set contained by the interval only
# contains a single element.
function Base.convert(::Type{T}, interval::AnchoredInterval{P,T}) where {P,T}
    if isclosed(interval) && (iszero(P) || first(interval) == last(interval))
        return first(interval)
    else
        # Remove deprecation in version 2.0.0
        depwarn(
            "`convert(T, interval::AnchoredInterval{P,T})` is deprecated for " *
            "intervals which are not closed with coinciding endpoints. " *
            "Use `anchor(interval)` instead.",
            :convert,
        )
        return anchor(interval)

        # TODO: For when deprecation is removed
        # throw(DomainError(interval, "The interval is not closed with coinciding endpoints"))
    end
end

function Base.convert(::Type{Interval}, interval::AnchoredInterval{P,T,L,R}) where {P,T,L,R}
    return Interval{T,L,R}(first(interval), last(interval))
end

function Base.convert(::Type{Interval{T}}, interval::AnchoredInterval{P,T,L,R}) where {P,T,L,R}
    return Interval{T,L,R}(first(interval), last(interval))
end

# Conversion methods which currently aren't needed but could prove useful. Commented out
# since these are untested.

#=
function Base.convert(::Type{AnchoredInterval{P,T}}, interval::Interval{T}) where {P,T}
    @assert abs(P) == span(interval)
    anchor = P < zero(P) ? last(interval) : first(interval)
    AnchoredInterval{P,T}(last(interval), inclusivity(interval))
end

function Base.convert(::Type{AnchoredInterval{P}}, interval::Interval{T}) where {P,T}
    @assert abs(P) == span(interval)
    anchor = P < zero(P) ? last(interval) : first(interval)
    AnchoredInterval{P,T}(anchor, inclusivity(interval))
end
=#

function Base.convert(::Type{AnchoredInterval{Ending}}, interval::Interval{T,L,R}) where {T,L,R}
    if !isbounded(interval)
        throw(ArgumentError("Unable to represent a non-bounded interval using a `AnchoredInterval`"))
    end
    AnchoredInterval{-span(interval), T, L, R}(last(interval))
end

function Base.convert(::Type{AnchoredInterval{Beginning}}, interval::Interval{T,L,R}) where {T,L,R}
    if !isbounded(interval)
        throw(ArgumentError("Unable to represent a non-bounded interval using a `AnchoredInterval`"))
    end
    AnchoredInterval{span(interval), T, L, R}(first(interval))
end

##### DISPLAY #####

function Base.show(io::IO, interval::T) where T <: AnchoredInterval
    if get(io, :compact, false)
        print(io, interval)
    else
        print(io, "$T(")
        show(io, anchor(interval))
        print(io, ")")
    end
end

function Base.print(io::IO, interval::AnchoredInterval{P,T}) where {P, T <: Union{Date, AbstractDateTime}}
    # Print to io in order to keep properties like :limit and :compact
    if get(io, :compact, false)
        io = IOContext(io, :limit=>true)
    end

    print(io, description(interval))
end

##### ARITHMETIC #####

Base.:+(a::T, b) where {T <: AnchoredInterval} = T(anchor(a) + b)

Base.:+(a, b::AnchoredInterval) = b + a
Base.:-(a::AnchoredInterval, b) = a + -b

# Required for StepRange{<:AnchoredInterval}
Base.:-(a::AnchoredInterval, b::AnchoredInterval) = anchor(a) - anchor(b)

Base.:-(a::T, b::AnchoredInterval{P,T}) where {P, T <: Number} = a + -b

function Base.:-(a::AnchoredInterval{P,T,L,R}) where {P, T <: Number, L, R}
    return AnchoredInterval{-P, T, R, L}(-anchor(a))
end

##### EQUALITY #####

function Base.:(==)(a::AnchoredInterval{P,T}, b::AnchoredInterval{P,T}) where {P,T}
    return anchor(a) == anchor(b) && bounds_types(a) == bounds_types(b)
end

# Required for min/max of AnchoredInterval{LaxZonedDateTime} when the anchor is AMB or DNE
function Base.isless(a::AnchoredInterval{P,T,L1}, b::AnchoredInterval{P,T,L2}) where {P,T,L1,L2}
    return (
        anchor(a) < anchor(b) ||
        (anchor(a) == anchor(b) && L1 === Closed && L2 === Open)
    )
end

##### RANGE #####

function Base.:(:)(start::AnchoredInterval{P,T}, step::S, stop::AnchoredInterval{P,T}) where {P,T,S}
    return StepRange{AnchoredInterval{P,T}, S}(start, step, stop)
end

# Infer step for two-argument StepRange{<:AnchoredInterval}
function Base.:(:)(start::AnchoredInterval{P,T}, stop::AnchoredInterval{P,T}) where {P,T}
    return (:)(start, abs(P), stop)
end

function Base.steprange_last(start::T, step, stop::AnchoredInterval) where T <: AnchoredInterval
    T(Base.steprange_last(anchor(start), step, anchor(stop)))
end

function Base.length(r::StepRange{<:AnchoredInterval})
    return length(anchor(r.start):r.step:anchor(r.stop))
end

##### SET OPERATIONS #####

function Base.isempty(interval::AnchoredInterval{P,T}) where {P,T}
    return P == zero(P) && !isclosed(interval)
end

# When intersecting two `AnchoredInterval`s attempt to return an `AnchoredInterval`
function Base.intersect(a::AnchoredInterval{P,T}, b::AnchoredInterval{Q,T}) where {P,Q,T}
    interval = invoke(intersect, Tuple{AbstractInterval{T}, AbstractInterval{T}}, a, b)

    sp = isa(P, Period) ? canonicalize(typeof(P), span(interval)) : span(interval)
    if P ≤ zero(P)
        anchor = last(interval)
        new_P = -sp
    else
        anchor = first(interval)
        new_P = sp
    end

    L, R = bounds_types(interval)
    return AnchoredInterval{new_P, T, L, R}(anchor)
end

##### UTILITIES #####

function canonicalize(target_type::Type{<:Period}, p::P) where P <: Period
    Q, max_val = coarserperiod(P)
    if (value(p) % max_val == 0) && (max_val > 1)
        p = canonicalize(target_type, convert(Q, p))
    end
    return p
end

canonicalize(target_type::Type{P}, p::P) where P <: Period = p

##### TIME ZONES #####

function TimeZones.astimezone(i::AnchoredInterval{P, ZonedDateTime, L, R}, tz::TimeZone) where {P,L,R}
    return AnchoredInterval{P, ZonedDateTime, L, R}(astimezone(anchor(i), tz))
end

TimeZones.timezone(i::AnchoredInterval{P, ZonedDateTime}) where P = timezone(anchor(i))


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 5
const MAXN = 500000

function wallis(n)
    res = BigFloat(1)
    for i = 1:div(n, 2)
        res *= BigFloat(4) * BigFloat(i * i) / (BigFloat(4) * BigFloat(i * i) - BigFloat(1))
    end

    return BigFloat(2) * res
end

function main()
    n = GetDigit_n.getn_even(wallis, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(wallis, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(wallis, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, wallis(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()