text:
"""
    Fock(quantum_system, equations, S, symmetries)

A Fock operator consists of a `quantum_system`, from which `equations`
are variationally derived, via an overload of
`diff(quantum_system)`. `quantum_system` must also provide an overload
for [`overlap_matrix`](@ref) -> `S`.

`equations` must provide an overload for
[`energy_matrix!`](@ref). `equations` must be iterable, where each
element corresponds to the equation for one orbital, and must provide
an overload for [`hamiltonian`](@ref) and
[`energy`](@ref). Additionally, [`update!`](@ref) must be provided for
`equations`, to prepare the equation system for the next iteration.
"""
mutable struct Fock{Q<:AbstractQuantumSystem,Equations,Overlap}
    quantum_system::Q
    equations::Equations
    S::Overlap
    symmetries::Vector{Vector{Int}}
end

Fock(quantum_system::Q; kwargs...) where Q =
    Fock(quantum_system, diff(quantum_system; kwargs...),
         overlap_matrix(quantum_system), symmetries(quantum_system))

function Base.show(io::IO, ::MIME"text/plain", fock::Fock{Q,E}) where {Q,E}
    write(io, "Fock operator with\n- quantum system: ")
    show(io, "text/plain", fock.quantum_system)
    write(io, "\n- SCF equations:  ")
    for eq in fock.equations
        write(io, "\n  - ")
        show(io, "text/plain", eq)
    end
end

Base.view(::Fock{Q}, args...) where Q =
    throw(ArgumentError("`view` not implemented for `Fock{$Q}`"))

"""
    overlap_matrix(quantum_system)

Return the overlap matrix which is a suitable metric for
`quantum_system`. Usually, this is `Sᵢⱼ = ⟨i|j⟩`, where `i` and `j`
are basis functions for the orbitals. This is used to calculate
inner products between orbitals, &c.
"""
overlap_matrix(::Q) where Q =
    throw(ArgumentError("`overlap_matrix` not implemented for `$Q`"))

symmetries(::Q) where Q =
    throw(ArgumentError("`symmetries` not implemented for `$Q`"))

"""
    energy_matrix!(H::AbstractMatrix, equations[, which=:total_energy])

Calculates the energy matrix of the system of `equations`. This
overwrites the entries of `H`. _To be overloaded by the user; must
support `which=:total_energy`, `which=:double_counted_energy`, and
`which=:kinetic_energy`._
"""
energy_matrix!(H::HM, equations::E, which::Symbol=:total_energy) where {HM<:AbstractMatrix,E} =
    throw(ArgumentError("`energy_matrix!` not implemented for $E"))

"""
    energy_matrix!(H::AbstractMatrix, fock::Fock[, which=:total_energy])

Calculates the energy matrix of the quantum system of `fock`. This
overwrites the entries of `H`.
"""
function energy_matrix!(H::HM, fock::F, which::Symbol=:total_energy) where {HM<:AbstractMatrix,F<:Fock}
    H .= zero(eltype(H))
    energy_matrix!(H, fock.equations, which)
    H
end

function energy(fock::Fock, which::Symbol)
    c = coefficients(fock.quantum_system)
    nc = length(c)
    T = zeros(eltype(fock.quantum_system), nc, nc)
    energy_matrix!(T, fock, which)
    dot(c, T, c)
end

orbital_energies(fock, which::Symbol=:total) =
    map(Base.Fix2(energy, which), fock.equations)

"""
    hamiltonian(equation::Equation)

Returns the orbital Hamiltonian of `equation`. _To be overloaded by the
user._
"""
hamiltonian(equation::Equation) where Equation =
    throw(ArgumentError("`hamiltonian` not implemented for `$Equation`"))

"""
    energy(equation::Equation[, term=:total_energy])

Returns the orbital energy for the orbital governed by `equation`. _To
be overloaded by the user._
"""
energy(equation::Equation, which::Symbol=:total_energy) where Equation =
    throw(ArgumentError("`energy` not implemented for `$Equation`"))

"""
    update!(eqs; kwargs...)

Update the equation system `eqs`, for the current iteration. _To be
overloaded by the user._
"""
update!(eqs; kwargs...) =
    throw(ArgumentError("`update!` not implemented for $(typeof(eqs))"))

"""
    update!(eqs, quantum_system; kwargs...)

Update the equation system `eqs` with respect to `quantum_system`, for
the current iteration. _To be overloaded by the user._
"""
update!(eqs, quantum_system; kwargs...) =
    throw(ArgumentError("`update!` not implemented for $(typeof(eqs)), $(typeof(quantum_system))"))

update!(fock::Fock; kwargs...) = update!(fock.equations)

"""
    rotate_first_lobe!(v)

Rotate the vector `v` such that the first lobe has positive sign.
"""
function rotate_first_lobe!(v::V) where {V<:AbstractVector}
    lmul!(sign(v[1]), v)
    v
end

"""
    norm_rot!(fock, v)

Normalize and rotate the eigenvector `v` such that the first
lobe has positive sign.
"""
function norm_rot!(fock::Fock, v::V) where {V<:AbstractVector}
    normalize!(fock.quantum_system, v)
    rotate_first_lobe!(v)
end

"""
    Matrix(fock)

Compute the Fock matrix, where each element corresponds to a matrix
element between two orbitals.
"""
function LinearAlgebra.Matrix(fock::Fock)
    m = length(fock.equations)
    F = Matrix{Float64}(undef, m, m)

    orbs = orbitals(fock.quantum_system)
    tmp = Vector{eltype(orbs)}(undef, size(orbs, 1))

    for (i,eqi) in enumerate(fock.equations)
        hi = KrylovWrapper(hamiltonian(eqi))
        orbi = view(orbs, :, i)

        for (j,eqj) in enumerate(fock.equations)
            mul!(tmp, hi, orbi)
            F[i,j] = real(dot(view(orbs, :, j), tmp))
        end
    end
    F
end


text:
####################################################################################################
# highgui. High-level GUI and Media I/O
####################################################################################################

# User Interface

# createTrackbar: Creates a trackbar and attaches it to the specified window
# 1. create a TrackbarCallback function in C++ and wrap with @cxx macro
# 2. call createTrackbar() below

createTrackbar(trackbarname::String, winname::String, value::Ptr{Int32}, count::Int) =
     @cxx cv::createTrackbar(pointer(trackbarname),pointer(winname), value, count)

# createTrackbar(trackbarname::String, winname::String, value::Int, count::Int, onChange=C_NULL,
#     userdata=Ptr{Void}[0][1]) = @cxx cv::createTrackbar(pointer(trackbarname),pointer(winname), pointer([cint(value)]), count, onChange, userdata)

# trackbarname – Name of the created trackbar
# winname      – Name of the window that will be used as a parent of the created trackbar
# value        – Optional pointer to an integer variable whose value reflects the position of the slider
# count        – Maximal position of the slider. The minimal position is always 0.
# onChange     – Pointer to the function (TrackbarCallback) to be called every time the slider changes position.
#                 This function should be prototyped as void Foo(int,void*)
# userdata     – User data (void*) that is passed as is to the callback

# getTrackbarPos: Returns the trackbar position
getTrackbarPos(trackbarname::String,winname::String) = @cxx cv::getTrackbarPos(pointer(trackbarname),pointer(winname))

# namedWindow(const String& winname, int flags=WINDOW_AUTOSIZE)
namedWindow(windowName::String, flags=WINDOW_AUTOSIZE) = @cxx cv::namedWindow(pointer(windowName), flags)

# void imshow(const String& winname, InputArray mat)
imshow(windowName::String, img) = @cxx cv::imshow(pointer(windowName), img)

# waitKey
waitkey(delay) = @cxx cv::waitKey(delay)

# destroy GUI window(s)
destroyWindow(windowName::String) = @cxx cv::destroyWindow(pointer(windowName))
destroyAllWindows() = @cxx cv::destroyAllWindows()

# MoveWindow
moveWindow(winname::String, x::Int, y::Int) = @cxx cv::moveWindow(pointer(winname), x, y)
# x – The new x-coordinate of the window
# y – The new y-coordinate of the window

# ResizeWindow
resizeWindow(winname::String, width::Int, height::Int) = @cxx cv::resizeWindow(pointer(winname), width, height)

# SetMouseCallback: Sets mouse handler for the specified window
# create MouseCallback function wrapped in @cxx
setMouseCallback(winname::String, onMouse, userdata=Ptr{Void}[0]) =
    @cxx cv::setMouseCallback(pointer(winname), onMouse, userdata)

# Gets the mouse-wheel motion delta (multiple of 120)
if is_windows()
    # EVENT_MOUSEWHEEL and EVENT_MOUSEHWHEEL
    getMouseWheelDelta(flags::Int) = @cxx getMouseWheelDelta(flags)
end

# setTrackbarPos (by value)
setTrackbarPos(trackbarname::String, winname::String, pos::Int) =
   @cxx cv::setTrackbarPos(pointer(trackbarname), pointer(winname), pos)

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "HIP"
version = v"4.2.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/ROCm-Developer-Tools/HIP/archive/rocm-$(version).tar.gz",
                  "ecb929e0fc2eaaf7bbd16a1446a876a15baf72419c723734f456ee62e70b4c24"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/HIP*/

# disable broken Pre-Compiled Headers (PCH)
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/disable-pch.patch"

# force hipcc to pass paths to clang
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/hipcc-force-paths.patch"

# help hipcc to find BB C++ headers+libs+crt
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/hipcc-bb-paths.patch"

# disable tests
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/disable-tests.patch"

apk add coreutils dateutils

mkdir build && cd build
ln -s ${prefix}/bin/clang ${prefix}/tools/clang
export HIP_CLANG_PATH=${prefix}/tools
cmake -DCMAKE_INSTALL_PREFIX=${prefix}/hip \
      -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DHIP_COMPILER=clang \
      -DHIP_PLATFORM=amd \
      -DHIP_RUNTIME=rocclr \
      -DROCM_PATH=${prefix} \
      -DHSA_PATH=${prefix} \
      -DLLVM_DIR="${libdir}/cmake/llvm" \
      -DClang_DIR="${libdir}/cmake/clang" \
      -DLLD_DIR="${libdir}/cmake/ldd" \
      ..
make -j${nproc}
make install

# Cleanup
rm ${prefix}/tools/clang

# link to .hipVersion so clang can find it
ln -s ../hip/bin/.hipVersion ${bindir}/.hipVersion
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libamdhip64"], :libamdhip64, ["hip/lib"]),
    ExecutableProduct("hipcc", :hipcc, "hip/bin"),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(; name="ROCmLLVM_jll", version=v"4.2.0")),
    Dependency("hsa_rocr_jll"),
    Dependency("ROCmDeviceLibs_jll"),
    Dependency("ROCmCompilerSupport_jll", v"4.2.0"),
    Dependency("ROCmOpenCLRuntime_jll", v"4.2.0"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               julia_compat="1.7",
               preferred_gcc_version=v"8")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Test
using CBinding

const SYSROOT = Sys.isapple() ? ["-isysroot", joinpath(strip(String(read(`xcrun xcode-select --print-path`))), "Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk")] : []

@testset "Examples using CBinding.jl" begin
	for entry in readdir(@__DIR__)
		endswith(entry, ".jl") || continue
		entry == "runtests.jl" && continue

		@testset "$(first(splitext(entry)))" begin
			@test (@eval include($(entry))) isa Module
		end
	end
end


text:
module Dashboards

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

function __init__()
  Genie.Renderer.Html.register_normal_element("st__dashboard", context = Genie.Renderer.Html)
end

export dashboard

function dashboard( args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)

  wrap() do
    Genie.Renderer.Html.st__dashboard(args...; kwargs...)
  end
end

const tagname = "st-dashboard"

function dashboard(elemid, content::Union{String,Vector}; partial::Bool = true, title::String = "Stipple Dashboard", head_content::String = "",
                    class::String = "", style::String = "", channel::String = Genie.config.webchannels_default_route, kwargs...)
  content = if isa(content, Vector)
    push!(pushfirst!(content, "<$tagname>"), "</$tagname>")
  else
    string("<$tagname>", content, "</$tagname>")
  end

  kwargs = NamedTuple(delete!(Stipple.OptDict(kwargs...), :id), :id, elemid)

  Stipple.layout(Genie.Renderer.Html.div(content; kwargs...); partial=partial, title=title, class=class,
                                                              head_content=head_content, style=style, channel=channel)
end

end


text:
module Augmentor

using RecipesBase
using Requires
using Images
using AffineTransforms
using ColorTypes
using PiecewiseAffineTransforms
using Interpolations

export

    rotate_expand,
    rotate_crop,

    DisplacementField,
        uniform_displacement,
        gaussian_displacement,

    DisplacementMesh,

    ImageTransformation,
        Either,
        NoOp,
        FlipX,
        FlipY,
        Rotate90,
        Rotate180,
        Rotate270,
        Resize,
        CropRatio,
        CropSize,
        RCropSize,
        RCropRatio,
        Crop,
        Zoom,
        Scale,
        RandomDisplacement,
        SmoothedDisplacement,

    multiplier,

    ImageSource,
        DirImageSource,

    Pipeline,

    transform

include("common.jl")

include("rotate.jl")
include("crop.jl")

include("imagetransformation.jl")

include("displacementfield.jl")
include("displacementmesh.jl")

include("transformations/flipx.jl")
include("transformations/flipy.jl")
include("transformations/resize.jl")
include("transformations/cropratio.jl")
include("transformations/cropsize.jl")
include("transformations/rcropsize.jl")
include("transformations/rcropratio.jl")
include("transformations/crop.jl")
include("transformations/zoom.jl")
include("transformations/scale.jl")
include("transformations/rotate90.jl")
include("transformations/rotate180.jl")
include("transformations/rotate270.jl")
include("transformations/randomdisplacement.jl")
include("transformations/smootheddisplacement.jl")

include("imagesource.jl")
include("dirimagesource.jl")

include("pipeline.jl")
include("linearpipeline.jl")

# Lazy loading for package integrations
@require MLDataUtils include("integrations/MLDataUtils.jl")

end


text: function f() x = 0 ; while x < 10; y = x; x = x + 1; end; return y end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function xscale!(simple_plot::SimplePlot, cur_symbol::Symbol)
  ( cur_symbol == :log10 ) && ( cur_symbol = :log )
  ( cur_symbol == :lin ) && ( cur_symbol = :linear )
  ( cur_symbol == :identity ) && ( cur_symbol = :linear )

  @assert cur_symbol in [:linear, :log]
  simple_plot.layout["xaxis"]["type"] = string(cur_symbol)

  if cur_symbol == :linear
    delete!(simple_plot.layout["xaxis"], "exponentformat")
  else
    simple_plot.layout["xaxis"]["exponentformat"] = "power"
  end

  validate!(simple_plot)
end

function yscale!(simple_plot::SimplePlot, cur_symbol::Symbol)
  ( cur_symbol == :log10 ) && ( cur_symbol = :log )
  ( cur_symbol == :lin ) && ( cur_symbol = :linear )
  ( cur_symbol == :identity ) && ( cur_symbol = :linear )

  @assert cur_symbol in [:linear, :log]
  simple_plot.layout["yaxis"]["type"] = string(cur_symbol)

  if cur_symbol == :linear
    delete!(simple_plot.layout["yaxis"], "exponentformat")
  else
    simple_plot.layout["yaxis"]["exponentformat"] = "power"
  end

  validate!(simple_plot)
end

function xlims!(simple_plot::SimplePlot, varargs...)
  @assert 1 <= length(varargs) <= 2

  if length(varargs) == 1
    cur_lims = first(varargs)
    @assert isa(cur_lims, Tuple)
  else
    cur_lims = tuple(varargs...)
  end

  simple_plot.layout["xaxis"]["autorange"] = false
  simple_plot.xlims = cur_lims

  validate!(simple_plot)
end

function ylims!(simple_plot::SimplePlot, varargs...)
  @assert 1 <= length(varargs) <= 2

  if length(varargs) == 1
    cur_lims = first(varargs)
    @assert isa(cur_lims, Tuple)
  else
    cur_lims = tuple(varargs...)
  end

  simple_plot.layout["yaxis"]["autorange"] = false
  simple_plot.ylims = cur_lims

  validate!(simple_plot)
end

xscale!(cur_symbol::Symbol) = xscale!(_plot, cur_symbol)
yscale!(cur_symbol::Symbol) = yscale!(_plot, cur_symbol)

xlims!(cur_lims...) = xlims!(_plot, cur_lims...)
ylims!(cur_lims...) = ylims!(_plot, cur_lims...)

export xscale!
export yscale!

export xlims!
export ylims!


text:
# Julia wrapper for header: /usr/include/libavutil/tree.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_tree_node_alloc,
    av_tree_find,
    av_tree_insert,
    av_tree_destroy,
    av_tree_enumerate

function av_tree_node_alloc()
    ccall((:av_tree_node_alloc, libavutil), Ptr{AVTreeNode}, ())
end

function av_tree_find(root, key, cmp, next::NTuple{2, Ptr{Cvoid}})
    ccall((:av_tree_find, libavutil), Ptr{Cvoid}, (Ptr{AVTreeNode}, Ptr{Cvoid}, Ptr{Cvoid}, NTuple{2, Ptr{Cvoid}}), root, key, cmp, next)
end

function av_tree_insert(rootp, key, cmp, next)
    ccall((:av_tree_insert, libavutil), Ptr{Cvoid}, (Ptr{Ptr{AVTreeNode}}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Ptr{AVTreeNode}}), rootp, key, cmp, next)
end

function av_tree_destroy(t)
    ccall((:av_tree_destroy, libavutil), Cvoid, (Ptr{AVTreeNode},), t)
end

function av_tree_enumerate(t, opaque, cmp, enu)
    ccall((:av_tree_enumerate, libavutil), Cvoid, (Ptr{AVTreeNode}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), t, opaque, cmp, enu)
end


text:
mutable struct MaxNNode
    snakes::Vector{Snake}
    food::Food
    height::Int
    width::Int
    depth::Int
    action::Union{Direction, Nothing}
    winner::Union{Int, Nothing}
    dead::Set{Int}
end

function iterative_deepening_maxn(board::Boardstate, simulation_time::Float64, height::Int, width::Int, you_index::Int)
        simulation_time = Millisecond(floor(Int,simulation_time))
        snakes = board.snakes
        food = board.food
        root = MaxNNode(snakes, food, height, width, 0, nothing, nothing, Set{Int}())
        action_nodes = possible_successors(root, you_index)
        best_action = nothing
        depth = 1

        times = []
        start = now()
        last_time = Millisecond(0)
        values = Tuple{Vector{Float64}, MaxNNode}[]
        while now() - start + last_time < simulation_time
            start_run = now()
            empty!(values)

            for a_node in action_nodes
                if now() - start < simulation_time
                    vector = maxn(a_node, depth, next_player(you_index, length(a_node.snakes)) , length(a_node.snakes), you_index)
                    push!(values, (vector, a_node))
                end
            end
            depth += 1
            if length(values) == length(action_nodes)
                sort!(values, by= x -> x[1][you_index], rev=true)
                best_action = values[1][2].action
                action_nodes = [v[2] for v in values]
            end
            last_time = now() - start_run
            push!(times, last_time)
        end
        println("$depth")
        println("times $times")
        println("Time used" ,now() - start)
        # println([(v[1], v[2].action) for v in values])
        return best_action
end

function maxn(node::MaxNNode, depth::Int, player::Int, num_players::Int, you_index::Int)
    if depth == 0 || node.winner !== nothing || length(node.dead) == length(node.snakes)
        return score(node)
    end

    while player in node.dead
        player = next_player(player, num_players)
    end
    next_p = next_player(player, num_players)
    max_vector = fill(-Inf, 4)
    next_nodes = possible_successors(node, player)
    for successor in next_nodes
        vector = maxn(successor, depth - 1, next_p, length(node.snakes), you_index)
        if vector[player] > max_vector[player]
            max_vector = vector
        end
    end

    return max_vector
end

function possible_successors(node::MaxNNode, player::Int)
    if player in node.dead
        return []
    end
    successors = MaxNNode[]
    snake = node.snakes[player]
    for a in possible_actions(snake)
        new_snek, new_food = deepcopy(snake), deepcopy(node.food)
        move_snake!(new_snek, new_food, a)
        if !check_bounds(get_head(new_snek), node.height, node.width)
            continue
        end
        snakes = deepcopy(node.snakes)
        snakes[player] = new_snek
        dead_for_a = deepcopy(node.dead)
        winner = check_winner(node, snakes, dead_for_a)
        n = MaxNNode(snakes, new_food, node.height, node.width, node.depth + 1, a, winner, dead_for_a)
        push!(successors, n)
    end
    return successors
end

function check_winner(node::MaxNNode, snakes::Vector{Snake}, dead::Set{Int})

    for (p, snek) in enumerate(snakes)
        # starve
        if snek.health == 0
            push!(dead, p)
        end

        for (o, opponent) in enumerate(snakes)
            if o == p
                continue
            end

            if get_head(snek) in opponent.body
                # head2head
                if get_head(snek) == get_head(opponent)
                    # RAM
                    if length(snek.body) > length(opponent.body)
                        push!(dead, o)
                    # draw
                    elseif length(snek.body) == length(opponent.body)
                        push!(dead, o)
                        push!(dead, p)
                    else
                        # lost (head2head) collision with opponent
                        push!(dead, p)
                    end
                else
                    # ran into opponent
                    push!(dead, p)
                end
            end

            # collision with self
            if get_head(snek) in get_body_no_head(snek)
                push!(dead, p)
            end
            if !check_bounds(get_head(snek), node.height, node.width)
                push!(dead, p)
            end

            if length(dead) == length(snakes) - 1
                winset = setdiff(Set(collect(1:length(snakes))), dead)
                return minimum(winset)
            end
        end
    end
end

function score(node::MaxNNode)
    result = zeros(length(node.snakes))
    health = [snek.health for snek in node.snakes]
    lens = [length(snek.body) for snek in node.snakes]
    flood_values = maxN_floodfill(node.snakes, node.width, node.height)
    for i in 1:length(node.snakes)
        score = 1 * health[i] + 4 * lens[i] + 1 * flood_values[i]
        result[i] += score
    end

    if !isnothing(node.winner)
        result[node.winner] = 1000
    end

    for d in node.dead
        result[d] = -1000
    end

    return result

end

max_stars_repo_name: rhotertj/solide-schlange

max_issues_repo_name: rhotertj/solide-schlange

max_forks_repo_name: rhotertj/solide-schlange


text:
module GradientTest

import Calculus

using Test
using ForwardDiff
using ForwardDiff: Dual, Tag
using StaticArrays
using DiffTests

include(joinpath(dirname(@__FILE__), "utils.jl"))

##################
# hardcoded test #
##################

f = DiffTests.rosenbrock_1
x = [0.1, 0.2, 0.3]
v = f(x)
g = [-9.4, 15.6, 52.0]

for c in (1, 2, 3), tag in (nothing, Tag(f, eltype(x)))
    println("  ...running hardcoded test with chunk size = $c and tag = $(repr(tag))")
    cfg = ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{c}(), tag)

    @test eltype(cfg) == Dual{typeof(tag), eltype(x), c}

    @test isapprox(g, ForwardDiff.gradient(f, x, cfg))
    @test isapprox(g, ForwardDiff.gradient(f, x))

    out = similar(x)
    ForwardDiff.gradient!(out, f, x, cfg)
    @test isapprox(out, g)

    out = similar(x)
    ForwardDiff.gradient!(out, f, x)
    @test isapprox(out, g)

    out = DiffResults.GradientResult(x)
    ForwardDiff.gradient!(out, f, x, cfg)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.gradient(out), g)

    out = DiffResults.GradientResult(x)
    ForwardDiff.gradient!(out, f, x)
    @test isapprox(DiffResults.value(out), v)
end

cfgx = ForwardDiff.GradientConfig(sin, x)
@test_throws ForwardDiff.InvalidTagException ForwardDiff.gradient(f, x, cfgx)
@test ForwardDiff.gradient(f, x, cfgx, Val{false}()) == ForwardDiff.gradient(f,x)

########################
# test vs. Calculus.jl #
########################

for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    v = f(X)
    g = ForwardDiff.gradient(f, X)
    @test isapprox(g, Calculus.gradient(f, X), atol=FINITEDIFF_ERROR)
    for c in CHUNK_SIZES, tag in (nothing, Tag(f, eltype(x)))
        println("  ...testing $f with chunk size = $c and tag = $(repr(tag))")
        cfg = ForwardDiff.GradientConfig(f, X, ForwardDiff.Chunk{c}(), tag)

        out = ForwardDiff.gradient(f, X, cfg)
        @test isapprox(out, g)

        out = similar(X)
        ForwardDiff.gradient!(out, f, X, cfg)
        @test isapprox(out, g)

        out = DiffResults.GradientResult(X)
        ForwardDiff.gradient!(out, f, X, cfg)
        @test isapprox(DiffResults.value(out), v)
        @test isapprox(DiffResults.gradient(out), g)
    end
end

##########################################
# test specialized StaticArray codepaths #
##########################################

println("  ...testing specialized StaticArray codepaths")

x = rand(3, 3)

for T in (StaticArrays.SArray, StaticArrays.MArray)
    sx = T{Tuple{3,3}}(x)

    cfg = ForwardDiff.GradientConfig(nothing, x)
    scfg = ForwardDiff.GradientConfig(nothing, sx)

    actual = ForwardDiff.gradient(prod, x)
    @test ForwardDiff.gradient(prod, sx) == actual
    @test ForwardDiff.gradient(prod, sx, cfg) == actual
    @test ForwardDiff.gradient(prod, sx, scfg) == actual
    @test ForwardDiff.gradient(prod, sx, scfg) isa StaticArray
    @test ForwardDiff.gradient(prod, sx, scfg, Val{false}()) == actual
    @test ForwardDiff.gradient(prod, sx, scfg, Val{false}()) isa StaticArray

    out = similar(x)
    ForwardDiff.gradient!(out, prod, sx)
    @test out == actual

    out = similar(x)
    ForwardDiff.gradient!(out, prod, sx, cfg)
    @test out == actual

    out = similar(x)
    ForwardDiff.gradient!(out, prod, sx, scfg)
    @test out == actual

    result = DiffResults.GradientResult(x)
    result = ForwardDiff.gradient!(result, prod, x)

    result1 = DiffResults.GradientResult(x)
    result2 = DiffResults.GradientResult(x)
    result3 = DiffResults.GradientResult(x)
    result1 = ForwardDiff.gradient!(result1, prod, sx)
    result2 = ForwardDiff.gradient!(result2, prod, sx, cfg)
    result3 = ForwardDiff.gradient!(result3, prod, sx, scfg)
    @test DiffResults.value(result1) == DiffResults.value(result)
    @test DiffResults.value(result2) == DiffResults.value(result)
    @test DiffResults.value(result3) == DiffResults.value(result)
    @test DiffResults.gradient(result1) == DiffResults.gradient(result)
    @test DiffResults.gradient(result2) == DiffResults.gradient(result)
    @test DiffResults.gradient(result3) == DiffResults.gradient(result)

    sresult1 = DiffResults.GradientResult(sx)
    sresult2 = DiffResults.GradientResult(sx)
    sresult3 = DiffResults.GradientResult(sx)
    sresult1 = ForwardDiff.gradient!(sresult1, prod, sx)
    sresult2 = ForwardDiff.gradient!(sresult2, prod, sx, cfg)
    sresult3 = ForwardDiff.gradient!(sresult3, prod, sx, scfg)
    @test DiffResults.value(sresult1) == DiffResults.value(result)
    @test DiffResults.value(sresult2) == DiffResults.value(result)
    @test DiffResults.value(sresult3) == DiffResults.value(result)
    @test DiffResults.gradient(sresult1) == DiffResults.gradient(result)
    @test DiffResults.gradient(sresult2) == DiffResults.gradient(result)
    @test DiffResults.gradient(sresult3) == DiffResults.gradient(result)
end

@testset "exponential function at base zero" begin
    @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0, -0.5]), [NaN, NaN])
    @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  0.0]), [NaN, NaN])
    @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  0.5]), [Inf, NaN])
    @test isequal(ForwardDiff.gradient(t -> t[1]^t[2], [0.0,  1.5]), [0.0, 0.0])
end

# Issue 399
@testset "chunk size zero" begin
    f_const(x) = 1.0
    g_grad_const = x -> ForwardDiff.gradient(f_const, x)
    @test g_grad_const([1.0]) == [0.0]
    @test isempty(g_grad_const(zeros(Float64, 0)))
end

@testset "dimension errors for gradient" begin
    @test_throws DimensionMismatch ForwardDiff.gradient(identity, 2pi) # input
    @test_throws DimensionMismatch ForwardDiff.gradient(identity, fill(2pi, 2)) # vector_mode_gradient
    @test_throws DimensionMismatch ForwardDiff.gradient(identity, fill(2pi, 10^6)) # chunk_mode_gradient
end

@testset "ArithmeticStyle" begin
    function f(p)
        sum(collect(0.0:p[1]:p[2]))
    end
    @test ForwardDiff.gradient(f, [0.2,25.0]) == [7875.0, 0.0]
end

end # module


text:
using Elf
using Documenter

makedocs(;
    modules=[Elf],
    authors="sasuseso <sinai471530@gmail.com> and contributors",
    repo="https://github.com/sasuseso/Elf.jl/blob/{commit}{path}#L{line}",
    sitename="Elf.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sasuseso.github.io/Elf.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sasuseso/Elf.jl",
)


text:
#=
Implementation of a metric.

From scratch.

Choose one metric to include,
if you also defined another metric having the same signature.
=#

using Colors    # RGB

include("metricFunction/negLnCauchy.jl")

#=
Specialized for types from Colors module.
Often: RGB{FixedPointNumbers.N0f8}
The specialization:
1) knows the components i.e. channels r,g,b
2) knows kinds of arithmetic (fixed point versus float) for channels
=#

using FixedPointNumbers

function pointDifference(
        targetColor::RGB{FixedPointNumbers.N0f8},
        corpusColor::RGB{FixedPointNumbers.N0f8}
        )
        # TODO generic for ValueType:   RGB{T} where {T}   # ::Float32

    #=
    Sum differences over the channels.

    Colors are not iterable over channels, but have getters r,g,b

    The type is one from Colors module.
    Often: FixedPointNumbers.N0f8
    See discussion of arithmetic on Julia docs.
    May convert to float.
    Saturated and wrapping arithmetic not available yet?

    It needs to not lose magnitude through adding plus and minus
    =#

    #=
    When T is FixedPointNumbers.N0f8,
    1) difference (-)  could wrap, overflow, or saturate and still fit in N0F8
    2) abs() might be better and also fit.
    TODO resolve whether - yields type Float

    Anyway, the sum over channels does not fit in N0F8
    (without losing accuracy due to overflow)
    So we sum to a float.
    =#

    sum::Float32 = zero(Float32)

    # typeof(targetColor.r)) is Float ???
    #=
    Use abs(), we want magnitude, not a signed difference.
    Depends on whether the cauchy function is computed from a symmetric table
    =#
    for value in (  targetColor.r - corpusColor.r,
                    targetColor.g - corpusColor.g,
                    targetColor.b - corpusColor.b,
                    )
        # sum += abs(value)
        # sum += negLnCauchy(abs(value))
        sum += negLnCauchy(value)
    end

    # Ensure the sum is positive,
    @assert sum >= 0
    return sum
end

#=
The maximum value of the metric.

The maximum value when one (target) value is zero and the other (corpus) value is
the max value of the difference between values of a channel,
summed over channels.

The maximum difference in a channel is 1.0 for type FixedPointNumbers.N0f8
=#

#=
All types for channels of RGB
are basically [0.0, 1.0] ????

=#
const maxPointDifference = 3 * negLnCauchy(1.0)

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module EON

import ..ED
using ..SimLib
using ..SimLib: Maybe
using LinearAlgebra
using SharedArrays: sdata

export EigenstateOccupation, EONDataDescriptor, EONData, load_eon

### Descriptor

struct EONDataDescriptor{T} <: ED.EDDerivedDataDescriptor
    ψ::Maybe{T}
    statename::String
    derivedfrom::ED.EDDataDescriptor
end

EONDataDescriptor(statename::String, args...; kwargs...) = EONDataDescriptor(nothing, statename, args...; kwargs...)
EONDataDescriptor(state, statename::String, args...; kwargs...) = EONDataDescriptor(state, statename, EDDataDescriptor(args...; kwargs...))

### Data obj

struct EONData{T, N} <: SimLib.AbstractSimpleData
    descriptor::EONDataDescriptor{T}
    data::Array{Float64, N}
end

ED._default_folder(::EONDataDescriptor) = "eon"
ED._filename_addition(opdd::EONDataDescriptor) = "_" * opdd.statename

"""
    load_eon(statename, edd)
    load_eon(statename, model[, diagtype] [, location])
"""
load_eon(args...; kwargs...) = load(EONDataDescriptor(args...; kwargs...))

### Task

mutable struct EONTask{T} <: ED.EDTask
    statename::String
    state::T
    data
end

EigenstateOccupation(statename, state) = EONTask(statename, state, nothing)

function ED.initialize!(task::EONTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, spectral_size)
end

function ED.compute_task!(task::EONTask, evals, evecs, inds...)
    n = min(size(task.data,1), size(evecs,2))
    for (i, vec) in enumerate(eachcol(evecs))
        i <= n || break
        task.data[i, inds...] = abs2(dot(vec, task.state))
    end
end

function ED.failed_task!(task::EONTask, inds...)
    task.data[:, inds...] .= NaN64
end

function ED.assemble(task::EONTask, edd)
    EONData(EONDataDescriptor(task.state, task.statename, edd), sdata(task.data))
end

Base.summary(task::EONTask) = string(typeof(task)) * "($(task.statename))"

end #module


text:
# # Kernel Ridge Regression
#
# Building on linear regression, we can fit non-linear data sets by introducing a feature space. In a higher-dimensional feature space, we can overfit the data; ridge regression introduces regularization to avoid this. In this notebook we show how we can use KernelFunctions.jl for *kernel* ridge regression.

## Loading and setup of required packages
using KernelFunctions
using LinearAlgebra
using Distributions

## Plotting
using Plots;
default(; lw=2.0, legendfontsize=11.0, ylims=(-150, 500));

using Random: seed!
seed!(42);

##
# ## Toy data
# Here we use a one-dimensional toy problem. We generate data using the fourth-order polynomial $f(x) = (x+4)(x+1)(x-1)(x-3)$:

f_truth(x) = (x + 4) * (x + 1) * (x - 1) * (x - 3)

x_train = -5:0.5:5
x_test = -7:0.1:7

noise = rand(Uniform(-20, 20), length(x_train))
y_train = f_truth.(x_train) + noise
y_test = f_truth.(x_test)

plot(x_test, y_test; label=raw"$f(x)$")
scatter!(x_train, y_train; seriescolor=1, label="observations")

##
# ## Linear regression
# For training inputs $\mathrm{X}=(\mathbf{x}_n)_{n=1}^N$ and observations $\mathbf{y}=(y_n)_{n=1}^N$, the linear regression weights $\mathbf{w}$ using the least-squares estimator are given by
# ```math
# \mathbf{w} = (\mathrm{X}^\top \mathrm{X})^{-1} \mathrm{X}^\top \mathbf{y}
# ```
# We predict at test inputs $\mathbf{x}_*$ using
# ```math
# \hat{y}_* = \mathbf{x}_*^\top \mathbf{w}
# ```
# This is implemented by `linear_regression`:

function linear_regression(X, y, Xstar)
    weights = (X' * X) \ (X' * y)
    return Xstar * weights
end;

# A linear regression fit to the above data set:

y_pred = linear_regression(x_train, y_train, x_test)
scatter(x_train, y_train; label="observations")
plot!(x_test, y_pred; label="linear fit")

##
# ## Featurization
# We can improve the fit by including additional features, i.e. generalizing to $\tilde{\mathrm{X}} = (\phi(x_n))_{n=1}^N$, where $\phi(x)$ constructs a feature vector for each input $x$. Here we include powers of the input, $\phi(x) = (1, x, x^2, \dots, x^d)$:

function featurize_poly(x; degree=1)
    return repeat(x, 1, degree + 1) .^ (0:degree)'
end

function featurized_fit_and_plot(degree)
    X = featurize_poly(x_train; degree=degree)
    Xstar = featurize_poly(x_test; degree=degree)
    y_pred = linear_regression(X, y_train, Xstar)
    scatter(x_train, y_train; legend=false, title="fit of order $degree")
    return plot!(x_test, y_pred)
end

plot((featurized_fit_and_plot(degree) for degree in 1:4)...)

##
# Note that the fit becomes perfect when we include exactly as many orders in the features as we have in the underlying polynomial (4).
#
# However, when increasing the number of features, we can quickly overfit to noise in the data set:

featurized_fit_and_plot(20)

##
# ## Ridge regression
# To counteract this unwanted behaviour, we can introduce regularization. This leads to *ridge regression* with $L_2$ regularization of the weights ([Tikhonov regularization](https://en.wikipedia.org/wiki/Tikhonov_regularization)).
# Instead of the weights in linear regression,
# ```math
# \mathbf{w} = (\mathrm{X}^\top \mathrm{X})^{-1} \mathrm{X}^\top \mathbf{y}
# ```
# we introduce the ridge parameter $\lambda$:
# ```math
# \mathbf{w} = (\mathrm{X}^\top \mathrm{X} + \lambda \mathbb{1})^{-1} \mathrm{X}^\top \mathbf{y}
# ```
# As before, we predict at test inputs $\mathbf{x}_*$ using
# ```math
# \hat{y}_* = \mathbf{x}_*^\top \mathbf{w}
# ```
# This is implemented by `ridge_regression`:

function ridge_regression(X, y, Xstar, lambda)
    weights = (X' * X + lambda * I) \ (X' * y)
    return Xstar * weights
end

function regularized_fit_and_plot(degree, lambda)
    X = featurize_poly(x_train; degree=degree)
    Xstar = featurize_poly(x_test; degree=degree)
    y_pred = ridge_regression(X, y_train, Xstar, lambda)
    scatter(x_train, y_train; legend=false, title="\$\\lambda=$lambda\$")
    return plot!(x_test, y_pred)
end

plot((regularized_fit_and_plot(20, lambda) for lambda in (1e-3, 1e-2, 1e-1, 1))...)

##
# ## Kernel ridge regression
# Instead of constructing the feature matrix explicitly, we can use *kernels* to replace inner products of feature vectors with a kernel evaluation: $\langle \phi(x), \phi(x') \rangle = k(x, x')$ or $\tilde{\mathrm{X}} \tilde{\mathrm{X}}^\top = \mathrm{K}$, where $\mathrm{K}_{ij} = k(x_i, x_j)$.
#
# To apply this "kernel trick" to ridge regression, we can rewrite the ridge estimate for the weights
# ```math
# \mathbf{w} = (\mathrm{X}^\top \mathrm{X} + \lambda \mathbb{1})^{-1} \mathrm{X}^\top \mathbf{y}
# ```
# using the [matrix inversion lemma](https://tlienart.github.io/pub/csml/mtheory/matinvlem.html#basic_lemmas)
# as
# ```math
# \mathbf{w} = \mathrm{X}^\top (\mathrm{X} \mathrm{X}^\top + \lambda \mathbb{1})^{-1} \mathbf{y}
# ```
# where we can now replace the inner product with the kernel matrix,
# ```math
# \mathbf{w} = \mathrm{X}^\top (\mathrm{K} + \lambda \mathbb{1})^{-1} \mathbf{y}
# ```
# And the prediction yields another inner product,
# ```math
# \hat{y}_* = \mathbf{x}_*^\top \mathbf{w} = \langle \mathbf{x}_*, \mathbf{w} \rangle = \mathbf{k}_* (\mathrm{K} + \lambda \mathbb{1})^{-1} \mathbf{y}
# ```
# where $(\mathbf{k}_*)_n = k(x_*, x_n)$.
#
# This is implemented by `kernel_ridge_regression`:

function kernel_ridge_regression(k, X, y, Xstar, lambda)
    K = kernelmatrix(k, X)
    kstar = kernelmatrix(k, Xstar, X)
    return kstar * ((K + lambda * I) \ y)
end;

# Now, instead of explicitly constructing features, we can simply pass in a `PolynomialKernel` object:

function kernelized_fit_and_plot(kernel, lambda=1e-4)
    y_pred = kernel_ridge_regression(kernel, x_train, y_train, x_test, lambda)
    if kernel isa PolynomialKernel
        title = string("order ", kernel.degree)
    else
        title = string(nameof(typeof(kernel)))
    end
    scatter(x_train, y_train; label=nothing)
    return plot!(x_test, y_pred; label=nothing, title=title)
end

plot((kernelized_fit_and_plot(PolynomialKernel(; degree=degree, c=1)) for degree in 1:4)...)

##
# However, we can now also use kernels that would have an infinite-dimensional feature expansion, such as the squared exponential kernel:

kernelized_fit_and_plot(SqExponentialKernel())


text:
@testset "516.longest-palindromic-subsequence.jl" begin
    @test longest_palindrome_subseq("bbbab") == 4
    @test longest_palindrome_subseq("cbbd") == 2
    @test longest_palindrome_subseq("onhxrcdszgcncrxkjktufmgljsrhyumyxholdogjoqjewlxgxncmhjpwtrklngzfyizwdvhpxgmvzswnsyjrynvvbmsciztxlgcofetdffwxlkdteknigzjnjzfbobvpxmqonvqzyrhalqxajqpfddspnrgvbaposcezgewcsddjwvoqknmfztyfdycgaomczthegdmt") == 60
end


text:
#=
Copyright (c) 2018-2019: David Sanders and Luis Benet.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

using TaylorSeries
using TaylorIntegration: jetcoeffs!
using TaylorModels: validated_step!, TaylorModel1, TaylorModelN
import IntervalArithmetic

function validated_integ(f!, qq0::AbstractVector{T}, δq0::IntervalBox{N,T},
        t0::T, tmax::T, orderQ::Int, orderT::Int, abstol::T;
        maxsteps::Int=500, parse_eqs::Bool=true,
        check_property::Function=(t, x)->true) where {N, T<:Real}

    # Set proper parameters for jet transport
    @assert N == get_numvars()
    dof = N

    # Some variables
    R = IntervalArithmetic.Interval{T}
    q0 = IntervalBox(qq0)
    t = t0 + Taylor1(orderT)
    tI = t0 + Taylor1(orderT+1)
    δq_norm = IntervalBox(IntervalArithmetic.Interval{T}(-1,1), Val(N))
    q0box = q0 .+ δq_norm

    # Allocation of vectors
    # Output
    tv    = Vector{T}(undef, maxsteps+1)
    xv    = Vector{IntervalBox{N,T}}(undef, maxsteps+1)
    xTM1v = Matrix{TaylorModel1{TaylorN{T}, T}}(undef, dof, maxsteps+1)

    # Internals: jet transport integration
    x     = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    dx    = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    xaux  = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    x0    = Vector{TaylorN{T}}(undef, dof)
    xTMN  = Vector{TaylorModelN{N,T,T}}(undef, dof)

    # Internals: Taylor1{Interval{T}} integration
    xI    = Vector{Taylor1{IntervalArithmetic.Interval{T}}}(undef, dof)
    dxI   = Vector{Taylor1{IntervalArithmetic.Interval{T}}}(undef, dof)
    xauxI = Vector{Taylor1{IntervalArithmetic.Interval{T}}}(undef, dof)
    x0I   = Vector{IntervalArithmetic.Interval{T}}(undef, dof)

    # Set initial conditions
    zI = zero(R)
    rem = Vector{IntervalArithmetic.Interval{T}}(undef, dof)

    @inbounds for i in eachindex(x)
        qaux = normalize_taylor(qq0[i] + TaylorN(i, order=orderQ), δq0, true)
        x[i] = Taylor1(qaux, orderT)
        dx[i] = x[i]
        x0[i] = copy(qaux)
        xTMN[i] = TaylorModelN(x[i][0], zI, q0, q0box)

        xI[i] = Taylor1(q0box[i], orderT+1)
        dxI[i] = xI[i]
        x0I[i] = qaux(δq_norm)
        rem[i] = zI

        xTM1v[i, 1] = TaylorModel1(deepcopy(x[i]), zI, zI, zI)
    end

    # Output vectors
    @inbounds tv[1] = t0
    @inbounds xv[1] = IntervalBox(evaluate(xTMN, δq_norm))

    # Determine if specialized jetcoeffs! method exists (built by @taylorize)
    parse_eqs = parse_eqs && (length(methods(jetcoeffs!)) > 2)
    if parse_eqs
        try
            jetcoeffs!(Val(f!), t, x, dx)
        catch
            parse_eqs = false
        end
    end

    # Integration
    nsteps = 1
    while t0 < tmax

        # Validated step of the integration
        δt = validated_step!(f!, t, x, dx, xaux, tI, xI, dxI, xauxI,
            t0, tmax, x0, x0I, xTMN, xv, rem, δq_norm,
            q0, q0box, nsteps, orderT, abstol, parse_eqs, check_property)

        # New initial conditions and time
        nsteps += 1
        t0 += δt
        @inbounds begin
            t[0] = t0
            tI[0] = t0
            tv[nsteps] = t0
            for i in eachindex(x)
                xTM1v[i, nsteps] = TaylorModel1(deepcopy(x[i]), rem[i],
                                        IntervalArithmetic.Interval(0, 0),
                                        IntervalArithmetic.Interval(0, δt))
                aux = x[i](δt)
                x[i]  = Taylor1(aux, orderT)
                dx[i] = Taylor1(zero(aux), orderT)
            end
        end

        if nsteps > maxsteps
            info("Maximum number of integration steps reached; exiting.")
            break
        end

    end

    return view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v, :, 1:nsteps)
end


text:
using Distributed
@everywhere using StochasticDiffEq, DiffEqDevTools, Test
@everywhere using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems
@everywhere importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear, prob_sde_additive
Random.seed!(100)
dts = 1./2.^(10:-1:2) #14->7 good plot

prob = prob_sde_linear
sim  = test_convergence(dts,prob,EM(),numMonte=Int(1e4),
                        weak_timeseries_errors=true,weak_dense_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
@test abs(sim.𝒪est[:weak_L2]-1) < 0.3
@test abs(sim.𝒪est[:weak_L∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil(),numMonte=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.3

prob = prob_sde_2Dlinear
sim  = test_convergence(dts,prob,EM(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.35

prob = prob_sde_additive
sim  = test_convergence(dts,prob,EM(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.35
sim5 = test_convergence(dts,prob,SRA(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim5.𝒪est[:weak_final]-2) < 0.3
@test abs(sim5.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim5.𝒪est[:weak_l∞]-2) < 0.3
sim6 = test_convergence(dts,prob,SRA1(),numMonte=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim6.𝒪est[:weak_final]-2) < 0.3
@test abs(sim6.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim6.𝒪est[:weak_l∞]-2) < 0.3


text:
function handle_step(obj::FieldODEStepHandler, arg0::FieldODEStateInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (FieldODEStateInterpolator, jboolean), arg0, arg1)
end

function init(obj::FieldODEStepHandler, arg0::FieldODEStateAndDerivative, arg1::RealFieldElement)
    return jcall(obj, "init", void, (FieldODEStateAndDerivative, RealFieldElement), arg0, arg1)
end


text:
function evaluate_lda!(func::Functional, ρ::Array{Float64}; kwargs...)
    @assert func.family == family_lda
    n_p = Int(length(ρ) / func.n_spin)

    sizes = Dict(:E => n_p, :Vρ => 2n_p, :V2ρ2 => 3n_p, :V3ρ3 => 4n_p)
    if func.n_spin == 1
        for key in keys(sizes)
            sizes[key] = n_p
        end
    end

    dargs = Dict(kwargs)
    for key in keys(dargs)
        if !haskey(sizes, key)
            throw(ArgumentError("Unknown keyword argument: $(string(key))"))
        end
        if length(dargs[key]) != sizes[key]
            throw(DimensionMismatch("Length of keyword argument $(string(key)) == " *
                                    "$(length(dargs[key])) does not agree with expected " *
                                    "value $(sizes[key])"))
        end
    end

    ccall((:xc_lda, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                    Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          func.pointer, n_p, ρ, get(dargs, :E, C_NULL),
          get(dargs, :Vρ, C_NULL), get(dargs, :V2ρ2, C_NULL), get(dargs, :V3ρ3, C_NULL)
    )
end

function evaluate_gga!(func::Functional, ρ::Array{Float64}, σ::Array{Float64}; kwargs...)
    @assert func.family == family_gga
    n_p = Int(length(ρ) / func.n_spin)

    n_σ = 3n_p
    sizes = Dict(:E => n_p, :Vρ => 2n_p, :Vσ => 3n_p,
                 :V2ρ2 => 3n_p, :V2ρσ => 6n_p, :V2σ2 => 6n_p,
                 :V3ρ3 => 4n_p, :V3ρ2σ => 9n_p, :V3ρσ2 => 12n_p, :V3σ3 => 10n_p)
    if func.n_spin == 1
        n_σ = n_p
        for key in keys(sizes)
            sizes[key] = n_p
        end
    end

    dargs = Dict(kwargs)
    for key in keys(dargs)
        if !haskey(sizes, key)
            throw(ArgumentError("Unknown keyword argument: $(string(key))"))
        end
        if length(dargs[key]) != sizes[key]
            throw(DimensionMismatch("Length of keyword argument $(string(key)) == " *
                                    "$(length(dargs[key])) does not agree with expected " *
                                    "value $(sizes[key])"))
        end
    end

    ccall((:xc_gga, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                    Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                                    Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                                    Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                                    Ptr{Float64}),
          func.pointer, n_p, ρ, σ,
          get(dargs, :E, C_NULL),    get(dargs, :Vρ, C_NULL),    get(dargs, :Vσ, C_NULL),
          get(dargs, :V2ρ2, C_NULL), get(dargs, :V2ρσ, C_NULL),  get(dargs, :V2σ2, C_NULL),
          get(dargs, :V3ρ3, C_NULL), get(dargs, :V3ρ2σ, C_NULL), get(dargs, :V3ρσ2, C_NULL),
          get(dargs, :V3σ3, C_NULL)
    )
end


text:
using HSL

using LinearAlgebra
using SparseArrays
using Test

include("../deps/deps.jl")

if @isdefined libhsl_ma97
    include("test_ma97.jl")
end


text:
module PsPot

#
# HGH pseudopotentials
#
export PsPot_HGH
export HGH_eval_Vloc_G
export HGH_eval_proj_G
export info_PsPot_HGH
include("PsPot/PsPot_HGH.jl")

end


text:
function add_observation(obj::UpdatingMultipleLinearRegression, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "addObservation", void, (Vector{jdouble}, jdouble), arg0, arg1)
end

function add_observations(obj::UpdatingMultipleLinearRegression, arg0::Vector{Vector{jdouble}}, arg1::Vector{jdouble})
    return jcall(obj, "addObservations", void, (Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1)
end

function clear(obj::UpdatingMultipleLinearRegression)
    return jcall(obj, "clear", void, ())
end

function get_n(obj::UpdatingMultipleLinearRegression)
    return jcall(obj, "getN", jlong, ())
end

function has_intercept(obj::UpdatingMultipleLinearRegression)
    return jcall(obj, "hasIntercept", jboolean, ())
end

function regress(obj::UpdatingMultipleLinearRegression)
    return jcall(obj, "regress", RegressionResults, ())
end

function regress(obj::UpdatingMultipleLinearRegression, arg0::Vector{jint})
    return jcall(obj, "regress", RegressionResults, (Vector{jint},), arg0)
end


text:
"""
    domain = Domain(; is_in_domain = whole_domain,
        domain_params = Dict(:n₁ => 1, :n₂ => 0, :F => 0),
        domain_type = :background,
        domain_ε::Function = piecewise_constant_ε,
        domain_F::Function = piecewise_constant_F,
        is_in_subdomain = Function[],
        subdomain_params = Dict{Symbol,Float64}[],
        subdomain_type = Symbol[],
        subdomain_ε = Function[],
        subdomain_F = Function[],
        lattice = Bravais(),
        which_asymptote = :none,
        which_waveguide = 0)

- `is_in_domain` is boolean function with inputs `x`, `y`, `domain_params`

INCOMPLETE DOCUMENTATION

domain = Domain(domain; :key1 => value1, :key2 => value2, ...)

New domain structure from old with modified fields
"""
struct Domain{TSH1,TDP,TDT,TDE,TDF,TSH2,TSP,TSD,TSE,TSF}
    is_in_domain::TSH1
    domain_params::TDP
    domain_type::TDT
    domain_ε::TDE
    domain_F::TDF

    is_in_subdomain::TSH2
    subdomain_params::TSP
    subdomain_type::TSD
    subdomain_ε::TSE
    subdomain_F::TSF

    lattice::BravaisLattice
    which_asymptote::Symbol
    which_waveguide::Int

    num_subdomains::Int

    function Domain(;
        is_in_domain::Tsh1 = Universe(),
        domain_params::Tdp = Dict(:n₁ => 1, :n₂ => 0, :F => 0),
        domain_type::Tdt = GenericDomain(),
        domain_ε::Tde = piecewise_constant_ε,
        domain_F::Tdf = piecewise_constant_F,
        is_in_subdomain::Tsh2 = (),
        subdomain_params::Tsp = (),
        subdomain_type::Tsd = (),
        subdomain_ε::Tse = (),
        subdomain_F::Tsf = (),
        lattice::BravaisLattice = BravaisLattice(),
        which_asymptote::Symbol = :none,
        which_waveguide::Int = 0
        ) where {Tsh1,Tdp,Tsh2<:Tuple,Tsp,Tsd,Tdt<:AbstractDomain,Tde<:Function,Tdf<:Function,Tse,Tsf}

        num_subdomains= 1 + length(is_in_subdomain)

        @assert !(which_asymptote==:none && which_waveguide!==0) "assigned waveguide id $(which_waveguide) to non-asymptotic domain"

        return new{Tsh1,Tdp,Tdt,Tde,Tdf,Tsh2,Tsp,Tsd,Tse,Tsf}(
                    is_in_domain, domain_params, domain_type, domain_ε, domain_F,
                    is_in_subdomain, subdomain_params, subdomain_type, subdomain_ε, subdomain_F,
                    lattice, which_asymptote, which_waveguide, num_subdomains
                )
    end

    function Domain(dom::Td;
        is_in_domain = dom.is_in_domain,
        domain_params = dom.domain_params,
        domain_type = dom.domain_type,
        domain_ε = dom.domain_ε,
        domain_F = dom.domain_F,
        is_in_subdomain = dom.is_in_subdomain,
        subdomain_params = dom.subdomain_params,
        subdomain_type = dom.subdomain_type,
        subdomain_ε = dom.subdomain_ε,
        subdomain_F = dom.subdomain_F,
        lattice = dom.lattice,
        which_asymptote = dom.which_asymptote,
        which_waveguide = dom.which_waveguide
        ) where Td<:Domain

        return Domain(;is_in_domain = is_in_domain, domain_params = domain_params, domain_type = domain_type,
                domain_ε = domain_ε, domain_F = domain_F, is_in_subdomain = is_in_subdomain, subdomain_params = subdomain_params,
                subdomain_type = subdomain_type, subdomain_ε = subdomain_ε, subdomain_F = subdomain_F, lattice = lattice,
                which_asymptote = which_asymptote, which_waveguide = which_waveguide)
    end

    function Base.show(io::IO, dom::Domain)
        print(io, "Domain: \n")
        print(io, "\tDomain type: ", dom.domain_type, "\n")
        if dom.which_asymptote !== :none
            print(io, "\t\twaveguide: ", dom.which_waveguide,"\n")
            print(io, "\t\tasymptotic region: ", dom.which_asymptote, "\n")
        end
        print(io, "\tdomain shape: ", dom.is_in_domain, "\n")
        print(IOContext(io, :typeinfo => Dict), "\tdomain params: ", dom.domain_params, "\n")
        print(IOContext(io, :sub=>true, :sub1=>true), "\tdomain lattice: ", dom.lattice, "\n")
        print(IOContext(io, :typeinfo => Array{Function}), "\tbackground dielectric function: ", dom.domain_ε, "\n")
        print(IOContext(io, :typeinfo => Array{Float64}), "\tbackground pump function : ", dom.domain_F, "\n")
            print(IOContext(io, :typeinfo => Array{AbstractShape}), "\tsubdomain shapes: ", dom.is_in_subdomain, "\n")
        print(IOContext(io, :typeinfo => Array{Dict{Symbol},1}), "\tsubdomain params: ", dom.subdomain_params)
    end
end

"""
    sys = System(domains::Array{Domain})
    sys = System(domain)

collection of domains that defines System, an input for Simulation.

input can be array or scalar.

    sys = System(sys)

system object from system object
"""
struct System{TDOM,TDIS,TPAR,TE,TF}
    domains::TDOM
    ε::Array{ComplexF64,2}
    F::Array{Float64,2}
    χ::TDIS

    domain_by_region::Array{Int,1}
    params_by_region::TPAR
    ε_by_region::TE
    F_by_region::TF

    regions::Array{Int,2}
    num_prev_regions::Array{Int,1}

    waveguides::Array{Int,1}

    function System(
                domains_in::Td = (Domain(),),
                χ::Tdis = TwoLevelSystem(),
                ε::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0),
                F::Array{Float64,2} = Array{Float64}(undef, 0, 0),
                regions::Array{Int,2} = Array{Int}(undef, 0, 0)
                ) where {Td<:Tuple,Tdis<:AbstractDispersion}

        domains1 = fsort(domains_in; by=isBackground)
        domains2 = bsort(domains1; by=isDefect)
        domains3 = fsort(domains2; by=isPC)
        domains = bsort(domains3; by=isWaveguide)

        num_regions = sum(map(x->x.num_subdomains,domains))

        params_by_region = join_tuple(map(d->(d.subdomain_params...,d.domain_params),domains)...)
        ε_by_region = join_tuple(map(d->(d.subdomain_ε...,d.domain_ε),domains)...)
        F_by_region = join_tuple(map(d->(d.subdomain_F...,d.domain_F),domains)...)
        domain_by_region = Array{Int}(undef, num_regions)
        num_prev_regions = collect(cumsum_tuple(map(x->x.num_subdomains,domains)))
        idx = 1
        for i ∈ eachindex(domains)
            for j ∈ 1:domains[i].num_subdomains-1
                domain_by_region[idx] = i ; idx += 1
            end
            domain_by_region[idx] = i ; idx += 1
        end
        waveguides = sort(unique(map(d->d.which_waveguide,domains)))[2:end]
        return new{Td,Tdis,typeof(params_by_region),typeof(ε_by_region),typeof(F_by_region)}(domains, ε, F, χ, domain_by_region, params_by_region, ε_by_region, F_by_region, regions, num_prev_regions, waveguides)
    end

    System(sys::Ts) where Ts<:System = System(sys.domains)

    function Base.show(io::IO, sys::System)
        if !get(io, :sub, false)
            print(io, "System with ", length(sys.domains), " domains: \n")
        end
        domain_string = [["\tdomain ", i, " type: ", sys.domains[i].domain_type] for i ∈ eachindex(sys.domains)]
        asymptote_string = [[", ", sys.domains[i].which_asymptote] for i ∈ eachindex(sys.domains)]
        waveguide_string = [[", waveguide ", sys.domains[i].which_waveguide, "\n"] for i ∈ eachindex(sys.domains)]
        for i ∈ eachindex(sys.domains)
            print(io, domain_string[i]...)
            if sys.domains[i].which_asymptote !== :none
                print(io, asymptote_string[i]...)
                print(io, waveguide_string[i]...)
            else
                print(io, "\n")
            end
        end
        for i ∈ eachindex(sys.waveguides)
            if i==1
                print(io,"\n")
            end
            print(IOContext(io, :typeinfo => Array{Int}), "\twaveguide ", sys.waveguides[i], " domains: ", findall(sys.waveguides[i].==[sys.domains[j].which_waveguide for j ∈ eachindex(sys.domains)]), "\n")
            if i==length(sys.waveguides)
                print(io,"\n")
            end
        end
    end
end

join_tuple(x::Tuple, y...) = (x..., join_tuple(y...)...)
join_tuple(x::Tuple) = x

fsort(a::Tuple;kwargs...) =
    ssort(ssort(a...;kwargs...)...;kwargs...)
bsort(a::Tuple;kwargs...) =
    reverse(fsort(reverse(fsort(a;kwargs...));kwargs...))
ssort(a,b...;by) = by(a.domain_type) ? (a,ssort(b...;by=by)...) : (ssort(b...;by=by)...,a)
ssort(a;kwargs...)=(a,)

cumsum_tuple(a::Tuple) = reverse(_cumsum_tuple(reverse(a)...))
_cumsum_tuple(a,b...) = (sum(b),_cumsum_tuple(b...)...)
_cumsum_tuple(a) = (0,)

"""
    dis = Discretization(dis; key1 => value1, key2 => value2...)
new discretization object from old, with modified fields.
"""
struct Discretization{TCS}
    dx::NTuple{2,Float64}
    coordinate_system::TCS
    sub_pixel_num::Int
    origin::NTuple{2,Float64}

    N_tr::NTuple{2,Int}
    N::NTuple{2,Int}
    dN::NTuple{2,NTuple{2,Int}}

    x::NTuple{2,Array{Float64,2}}
    x_tr::NTuple{2,Array{Float64,2}}
    x_idx::NTuple{2,Array{Int,1}}

    X::Array{Float64,2}
    Y::Array{Float64,2}

    X_idx::Array{Int,1}

    function Discretization(dx::Tuple, coordinate_system::CS, sub_pixel_num::Int, origin::Tuple, N::Tuple, dN::NTuple{2,Tuple}) where CS

        N_tr = (N[1] - dN[1][1] - dN[1][2], N[2] - dN[2][1] - dN[2][2])

        x_tr = (Array{Float64}(undef,N_tr[1],1), Array{Float64}(undef,1,N_tr[2]))
        x = (Array{Float64}(undef,N[1],1), Array{Float64}(undef,1,N[2]))
        x_idx = (Array{Int}(undef, N_tr[1]), Array{Int}(undef, N_tr[2]))
        for j ∈ eachindex(N)
            if !isinf(dx[j])
                for i ∈ eachindex(x[j])
                    x[j][i] = origin[j] + dx[j]*(i+1/2)#(1/2:N[j]-1/2)
                end
                for i ∈ eachindex(x_tr[j])
                    x_tr[j][i] = x[j][1] + dN[j][1]*dx[j] + dx[j]*(i)#collect(0:N_tr[j]-1).*dx[j]
                end
                for i ∈ eachindex(x_idx[j])
                    x_idx[j][i] = dN[j][1] + i#(1:N_tr[j])
                end
            else
                x[j][:] .= origin[j]
                x_tr[j][:] .= x[j][1]
                x_idx[j][:] .= dN[j][1] .+ (1:N_tr[j])
            end
        end

        @assert CS ∈ (Cartesian,Polar) "unrecognized coordinate system type $CS"
        if CS == Cartesian
            X, Y = broadcast((x,y)->x,x[1],x[2]), broadcast((x,y)->y,x[1],x[2])
        else
            X = x[1].*cos.(x[2])
            Y = x[1].*sin.(x[2])
        end

        X_idx = LinearIndices(Array{Bool}(undef, N...))[x_idx...][:]
        return new{CS}(float.(dx), coordinate_system, sub_pixel_num, float.(origin), N_tr, N, dN,
            x, x_tr, x_idx, X, Y, X_idx)
    end

    Discretization(dis::Discretization{CS}; dx=dis.dx, coordinate_system=CS(), sub_pixel_num=dis.sub_pixel_num) where CS = Discretization(dx, coordinate_system, sub_pixel_num, (0.,0.), (1,1), ((0,0),(0,0)))

    function Base.show(io::IO, dis::Discretization)
        if !get(io, :sub, false)
            print(io, "Discretization: \n")
        end
        print(io, "\tN: ", dis.N, "\n",
        "\tsub-pixel number: ", dis.sub_pixel_num, "\n",
        "\tdx: (",fmt("1.4f",dis.dx[1]),",",fmt("1.4f",dis.dx[2]),")\n")
        if isCartesian(dis.coordinate_system)
            print(io, "\tcoordinates: Cartesian")
        elseif isPolar(dis.coordinate_system)
            print(io, "\tcoordinates: Polar")
        else
            print(io, "\tcoordinates: ", dis.coordinate_system)
        end
    end
end

"""
    bnd = Boundary(bnd; :key1 => value1, :key2 => value2, ...)

new boundary object from old, with modified fields.
"""
struct Boundary{TBC11,TBC12,TBC21,TBC22,TBL11,TBL12,TBL21,TBL22}
    ∂Ω::NTuple{2,NTuple{2,Float64}}
    ∂Ω_tr::NTuple{2,NTuple{2,Float64}}
    bc::Tuple{Tuple{TBC11,TBC12},Tuple{TBC21,TBC22}}
    bl::Tuple{Tuple{TBL11,TBL12},Tuple{TBL21,TBL22}}

    function Boundary(
                ∂Ω::NTuple{2,NTuple{2,Number}},
                BC::Tuple{Tuple{TBC1,TBC2},Tuple{TBC3,TBC4}} where {TBC1<:AbstractBC,TBC2<:AbstractBC,TBC3<:AbstractBC,TBC4<:AbstractBC},
                BL::Tuple{Tuple{TBL1,TBL2},Tuple{TBL3,TBL4}} where {TBL1<:AbstractBL,TBL2<:AbstractBL,TBL3<:AbstractBL,TBL4<:AbstractBL};
                warnings::Bool=true
                )
        if warnings
            (typeof(BL[1][1])<:Union{PML,cPML} && !(typeof(BC[1][1])<:DirichletBC)) ? (@warn "absorbing layer with non-dirichlet condition for dim 1, side 1, which is inconsistent") : nothing
            (typeof(BL[1][2])<:Union{PML,cPML} && !(typeof(BC[1][2])<:DirichletBC)) ? (@warn "absorbing layer with non-dirichlet condition for dim 1, side 2, which is inconsistent") : nothing
            (typeof(BL[2][1])<:Union{PML,cPML} && !(typeof(BC[2][1])<:DirichletBC)) ? (@warn "absorbing layer with non-dirichlet condition for dim 2, side 1, which is inconsistent") : nothing
            (typeof(BL[2][2])<:Union{PML,cPML} && !(typeof(BC[2][2])<:DirichletBC)) ? (@warn "absorbing layer with non-dirichlet condition for dim 2, side 2, which is inconsistent") : nothing
        end

        ∂Ω_tr1 = (∂Ω[1][1] + BL[1][1].depth, ∂Ω[1][2] - BL[1][2].depth)
        ∂Ω_tr2 = (∂Ω[2][1] + BL[2][1].depth, ∂Ω[2][2] - BL[2][2].depth)
        ∂Ω_tr = (∂Ω_tr1, ∂Ω_tr2)

        new{typeof(BC[1][1]),typeof(BC[1][2]),typeof(BC[2][1]),typeof(BC[2][2]),typeof(BL[1][1]),typeof(BL[1][2]),typeof(BL[2][1]),typeof(BL[2][2])}(map(x->map(float,x),∂Ω),map(x->map(float,x),∂Ω_tr),BC,BL)
    end

    Boundary(bnd::Boundary; ∂Ω=bnd.∂Ω, bc=bnd.bc, bl=bnd.bl) = Boundary(∂Ω, bc, bl)

    Base.show(io::IO, bnd::Boundary) = begin
        if !get(io, :sub, false)
            print(io, "Boundary: \n")
        end

        print(io, "\t\t\t-------------------------\n")

        print(io, "\tbound.\t|\t\t   ", fmt("6s",_get_b_symbol(bnd.bc[2][2])), "\t\t|\n")
        print(io, "\tcond-\t|   ",fmt("6s",_get_b_symbol(bnd.bc[1][1])), "\t\t   ", fmt("5s",_get_b_symbol(bnd.bc[1][2])), "|\n")
        print(io, "\t ition\t|\t\t   ", fmt("6s",_get_b_symbol(bnd.bc[2][1])), "\t\t|\n")

        print(io, "\t\t\t-------------------------\n")

        print(io, "\t\t\t|\t\t", fmt("+5.3f",bnd.∂Ω[2][2]), "\t\t\t|\n")
        print(io, "\t∂Ω\t\t| ", fmt("+5.3f",bnd.∂Ω[1][1]), "\t\t", fmt("+5.3f",bnd.∂Ω[1][2]), "\t|\n")
        print(io, "\t\t\t|\t\t", fmt("+5.3f",bnd.∂Ω[2][1]), "\t\t\t|\n")

        print(io, "\t\t\t-------------------------\n")

        print(io, "\t\t\t|\t\t", fmt("+5.3f",bnd.∂Ω_tr[2][2]), "\t\t\t|\n")
        print(io, "\t∂Ω_tr\t| ",fmt("+5.3f",bnd.∂Ω_tr[1][1]), "\t\t", fmt("+5.3f",bnd.∂Ω_tr[1][2]), "\t|\n")
        print(io, "\t\t\t|\t\t", fmt("+5.3f",bnd.∂Ω_tr[2][1]), "\t\t\t|\n")

        print(io, "\t\t\t-------------------------\n")

        print(io, "\tbound.\t|\t\t", fmt("8s",_get_b_symbol(bnd.bl[2][2])), "\t\t|\n")
        print(io, "\tlayer\t|",fmt("8s",_get_b_symbol(bnd.bl[1][1])), "\t\t", fmt("8s",_get_b_symbol(bnd.bl[1][2])), "|\n")
        print(io, "\t\t\t|\t\t", fmt("8s",_get_b_symbol(bnd.bl[2][1])), "\t\t|\n")

        print(io, "\t\t\t-------------------------\n")

        print(io, "\tbound.\t|\t\t", fmt("+5.3f",bnd.bl[2][2].depth), "\t\t\t|\n")
        print(io, "\tlayer\t| ",fmt("+5.3f",bnd.bl[1][1].depth), "\t\t", fmt("+5.3f",bnd.bl[1][2].depth), "\t|\n")
        print(io, "\tdepth\t|\t\t", fmt("+5.3f",bnd.bl[2][1].depth), "\t\t\t|\n")

        print(io, "\t\t\t-------------------------\n")
    end
    function _get_b_symbol(bcl::T) where T<:Union{AbstractBC,AbstractBL}
        if T<:AbstractBC
            if T<:FloquetBC
                return :p
            elseif T<:RobinBC
                return :r
            elseif T<:DirichletBC
                return :d
            elseif T<:NeumannBC
                return :n
            elseif T<:MatchedBC
                return :m
            else
                throw(ArgumentError("unrecognized boundary condition"))
            end
        else
            if T<:PML
                return :PML
            elseif T<:cPML
                return :cPML
            elseif T<:noBL
                return :none
            else
                throw(ArgumentError("unrecognized boundary layer"))
            end
        end
    end
end

"""
    channel = Channels(waveguide=0, quantum_number=0)

channel object for Scattering().

`waveguide` identifies which waveguide the channel is defined by

`quantum_number` identifies which mode of the waveguide defines the channel

    channel = Channels(chn; key1 => value1, key2 => value2, ...)

new channel object from old with modified fields.
"""
struct Channels
    waveguide::Int
    quantum_number::Int
    dispersion::Array{AbstractInterpolation,1}
    gaps::Array{Array{Float64,1},1}

    Channels(waveguide=0, quantum_number=0) = new(waveguide, quantum_number, AbstractInterpolation[], Array{Float64,1}[])

    Channels(chn::Channels; waveguide=chn.waveguide, quantum_number=chn.quantum_number) = Channels(waveguide, quantum_number)

    function Base.show(io::IO, chn::Channels)
        if get(io, :indented, false)
            print(io, "\t\twaveguide: ", chn.waveguide, "\n",
            "\t\tquantum number: ", chn.quantum_number)
        else
            print(io, "Channel: \n")
            print(io, "\twaveguide: ", chn.waveguide, "\n",
            "\tquantum number: ", chn.quantum_number, "\n")
        end
    end
end

"""
    sct = Scattering(channels::Array{Channel})

scattering object for Simulation

`channels` array of Channel objects

    sct = Scattering(sct; :key1 => value1, :key2 => value2, ...)

new scattering object from old with modified fields.
"""
struct Scattering
    channels::Array{Channels,1}
    waveguides_used::Array{Int,1}

    ε₀::Array{Array{ComplexF64,2},1}

    function Scattering(channels::Array{Channels,1}=Channels[], ε₀=[])

        channels = deepcopy(channels)
        ε₀ = deepcopy(ε₀)

        num_channels = length(channels)
        waveguides = Array{Int}(undef,num_channels)
        for i ∈ eachindex(waveguides)
            waveguides[i]=channels[i].waveguide
        end
        waveguides_used = unique(waveguides)
        ε₀ = Array{Array{ComplexF64,2},1}(undef,length(waveguides_used))
        return new(channels, waveguides_used, ε₀)
    end

    Scattering(sct::Scattering; channels=sct.channels) = Scattering(channels)

    function Base.show(io::IO, sct::Scattering)
        plot_flag=false
        if !get(io, :sub, false)
            plot_flag = true
            print(io, "Scattering with ", length(sct.channels), " channels:\n")
        end
        temp = [["\n", sct.channels[i]] for i ∈ eachindex(sct.channels)]
        for i ∈ eachindex(temp)
            if i>1
                print(io, "\n")
            end
            print(io, "\tChannel ", i, ":")
            print(IOContext(io::IO, :indented => true), temp[i]...)
        end
    end
end

"""
See also: [`System`](@ref), [`Boundary`](@ref), [`Discretization`](@ref), [`Scattering`](@ref), [`TwoLevelSystem`](@ref), [`Bravais`](@ref)

    sim = Simulation(; sys, bnd, dis, sct=Scattering(), tls=TwoLevelSystem(), lat=Bravais(bnd))

simulation object
"""
struct Simulation{TSYS,TBND,TDIS}
    sys::TSYS
    bnd::TBND
    dis::TDIS
    sct::Scattering
    lat::BravaisLattice

    function Simulation(
        bnd::Boundary,
        dis::Discretization{CS},
        sys::System,
        sct::Scattering,
        lat::BravaisLattice;
        display::Bool=false,
        warnings::Bool=true,
        smoothing::Bool=true
        ) where CS<:CoordinateSystem

        ∂Ω = bnd.∂Ω
        L = (∂Ω[1][2]-∂Ω[1][1],∂Ω[2][2]-∂Ω[2][1])

        dx = dis.dx
        @assert all(L./dx .< Inf) "lattice spacings $dx inconsistent with size $L"
        N = round.(Int,L./dx)

        dN = Array{Int}(undef,2,2)
        index1 = findmin(N)[2]
        index2 = mod1(index1+1,2)
        @assert !iszero(L[index2]) "no width along dimension $index2"
        if N[index1]≤1
            # N[index1] = 1
            # DX = [Inf, L[index2]/N[index2]]
            # dx = DX[[index1,index2]]
            # dN[:,index1] .= 0
            # dN[:,index2] .= floor.(Int, bnd.bl_depth[:,index2]/dx[index2])
            # for side ∈ 1:2
            #     bnd.BL[index2][side].depth = dN[side,index2]*dx[index2]
            # end
        else
            DX = [L[index1]/N[index1], L[index2]/round(Int,L[index2]/dx[index2])]
            dx = (DX[index1],DX[index2])
            N = round.(Int,L./dx)
            dN = ()
            for dim ∈ 1:2
                dNside = ()
                for side ∈ 1:2
                    dNside = (dNside...,floor(Int, bnd.bl[dim][side].depth/dx[dim]))
                end
                dN = (dN...,dNside)
            end

            for dim ∈ 1:2, side ∈ 1:2
                bnd.bl[dim][side].depth = dN[dim][side]*dx[dim]
            end
        end

        @assert all(sum.(dN) .< N) "boundary layers fill entire region. increase size of domain or decrease size of boundary layers."

        origin = (bnd.∂Ω[1][1],bnd.∂Ω[2][1])
        if typeof(bnd.bc[1][1])<:MatchedBC || typeof(bnd.bc[1][2])<:MatchedBC
            if typeof(bnd.bc[2][1])<:MatchedBC || typeof(bnd.bc[2][2])<:MatchedBC
                throw()
            else
                non_matched_dim = 2
            end
        else
            if typeof(bnd.bc[2][1])<:MatchedBC || typeof(bnd.bc[2][2])<:MatchedBC
                throw()
            else
                non_matched_dim = 1
            end
        end

        bb = apply_args(bnd.bc, CS(), bnd.bc[non_matched_dim][1], bnd.bc[non_matched_dim][2]; N=N, dx=dx, xmin=origin, lattice=lat)
        bnd = Boundary(∂Ω, bb, apply_args(bnd.bl; Δ=∂Ω))
        dis = Discretization(dx, dis.coordinate_system, dis.sub_pixel_num, origin, N, dN)

        # argument checking
        for side ∈ 1:2
            if warnings
                for dim ∈ 1:2
                    typeof(bnd.bc[dim][side])<:FloquetBC ? (typeof(bnd.bc[dim][mod1(side+1,2)])<:FloquetBC ? nothing : (@warn "only one side=$side of dim=$dim set to floquet")) : nothing
                end
            end
            @assert ((typeof(bnd.bc[1][side])<:FloquetBC && !isinf(lat.a)) || !(typeof(bnd.bc[1][side])<:FloquetBC)) "floquet bc in dim=1, but infinite lattice constant a"
            @assert ((typeof(bnd.bc[2][side])<:FloquetBC && !isinf(lat.b)) || !(typeof(bnd.bc[1][side])<:FloquetBC)) "floquet bc in dim=2, but infinite lattice constant b"
        end
        if warnings
            for i ∈ eachindex(sys.domains)
                if N[1] == 1 && !isinf(sys.domains[i].lattice.a)
                    @warn "one-dimensional (vertical) system, but domain $(i) has finite transverse lattice constant $(sys.domains[i].lattice.a) < ∞"
                end
                if N[2] == 1 && !isinf(sys.domains[i].lattice.b)
                    @warn "one-dimensional (horizontal) system, but domain $(i) has finite transverse lattice constant $(sys.domains[i].lattice.b) < ∞"
                end
            end
        end

        ε, F, regions = construct_εFr(bnd,dis,sys)
        smoothing ? sub_pixel_smoothing!(bnd, dis, sys, ε, F, regions; display=display) : nothing

        sys = System(sys.domains, sys.χ, ε, F, regions)

        which_waveguides = Array{Int,1}(undef,length(sys.domains))
        for j ∈ eachindex(sys.domains)
            which_waveguides[j] = sys.domains[j].which_waveguide
        end
        for i ∈ 1:length(sct.waveguides_used)
            sct.ɛ₀[i] = Array{ComplexF64}(undef,dis.N[1],dis.N[2])
            wg_inds = findall(which_waveguides .== sct.waveguides_used[i])
            wg_domains = Array{Domain}(undef,length(wg_inds))
            for j ∈ eachindex(wg_domains)
                wg_domains[j] = Domain(sys.domains[wg_inds[j]]; :which_asymptote => :none, :which_waveguide => 0)
            end
            if isempty(wg_domains)
                @warn "channels references undefined waveguide $(sct.waveguides_used[i]), scattering calculations will fail"
            else
                sct.ɛ₀[i][:] = sub_pixel_smoothing(bnd, dis, System(wg_domains), display=display)[1][:]
            end
        end
        return new{typeof(sys), typeof(bnd), typeof(dis)}(sys, bnd, dis, sct, lat)
    end

    Simulation(sim::Simulation; sys=sim.sys, bnd=sim.bnd, dis=sim.dis, sct=sim.sct, tls=sim.tls, lat=sim.lat, disp_opt=false) = deepcopy(Simulation(sys=System(sys), bnd=Boundary(bnd), dis=Discretization(dis), sct=Scattering(sct), tls=TwoLevelSystem(tls), lat=BravaisLattice(lat), disp_opt=disp_opt))

    function Base.show(io::IO, sim::Simulation)
        print(IOContext(io, :sub=>true),
        "Simulation : \n\n",
        "sys: \n", sim.sys, "\n\n",
        "bnd: \n", sim.bnd, "\n\n",
        "dis: \n", sim.dis, "\n\n",
        "sct: \n", sim.sct)
        if [sim.lat.a, sim.lat.b] !== [Inf,Inf]
            print(IOContext(io, :sub=>true, :sub1=>false), "\n\nlat: ", sim.lat)
        end
    end
end