text:
using CairoMakie
using GraphMakie
using GraphMakie.Graphs
using GraphMakie.NetworkLayout
using Makie.Colors
# using GLMakie
using Test

include("beziercurves_test.jl")

@testset "GraphMakie.jl" begin
    g = Node(wheel_digraph(10))
    f, ax, p = graphplot(g)
    f, ax, p = graphplot(g, node_attr=Attributes(visible=false))
    f, ax, p = graphplot(g, node_attr=(;visible=true))

    # try to update graph
    add_edge!(g[], 2, 4)
    g[] = g[]
    # try to update network
    p.layout = SFDP()

    # update node observables
    p.node_color = :blue
    p.node_size = 30

    # update edge observables
    p.edge_width = 5.0
    p.edge_color = :green

    # it should be also possible to pass multiple values
    f, ax, p = graphplot(g, node_color=[rand([:blue,:red,:green]) for i in 1:nv(g[])])
end

@testset "graph without edges" begin
    g = SimpleDiGraph(10)
    graphplot(g)
    g = SimpleGraph(10)
    graphplot(g)
end

@testset "selfedge without neighbors" begin
    g = SimpleGraph(10)
    add_edge!(g, 1, 1)
    graphplot(g)
end

@testset "small graphs" begin
    g = SimpleGraph(1)
    fig, ax, p = graphplot(g)

    g = complete_graph(2)
    graphplot(g)
end

@testset "single line width per edge" begin
    g = complete_graph(3)
    graphplot(g)
    graphplot(g; edge_width=10)
    graphplot(g; edge_width=[5, 10, 15])
    graphplot(g; edge_width=[5, 10, 5, 10, 5, 10])
end

@testset "Hover, click and drag Interaction" begin
    g = wheel_graph(10)
    add_edge!(g, 1, 1)
    f, ax, p = graphplot(g,
                         edge_width = [3.0 for i in 1:ne(g)],
                         edge_color = [colorant"black" for i in 1:ne(g)],
                         node_size = [20 for i in 1:nv(g)],
                         node_color = [colorant"red" for i in 1:nv(g)])
    deregister_interaction!(ax, :rectanglezoom)

    function edge_hover_action(s, idx, event, axis)
        p.edge_width[][idx]= s ? 6.0 : 3.0
        p.edge_width[] = p.edge_width[] # trigger observable
    end
    ehover = EdgeHoverHandler(edge_hover_action)
    register_interaction!(ax, :ehover, ehover)

    function node_hover_action(s, idx, event, axis)
        p.node_size[][idx] = s ? 40 : 20
        p.node_size[] = p.node_size[] # trigger observable
    end
    nhover = NodeHoverHandler(node_hover_action)
    register_interaction!(ax, :nhover, nhover)

    function node_click_action(idx, args...)
        p.node_color[][idx] = rand(RGB)
        p.node_color[] = p.node_color[]
    end
    nclick = NodeClickHandler(node_click_action)
    register_interaction!(ax, :nclick, nclick)

    function edge_click_action(idx, args...)
        p.edge_color[][idx] = rand(RGB)
        p.edge_color[] = p.edge_color[]
    end
    eclick = EdgeClickHandler(edge_click_action)
    register_interaction!(ax, :eclick, eclick)

    function node_drag_action(state, idx, event, axis)
        p[:node_pos][][idx] = event.data
        p[:node_pos][] = p[:node_pos][]
    end
    ndrag = NodeDragHandler(node_drag_action)
    register_interaction!(ax, :ndrag, ndrag)

    mutable struct EdgeDragAction
        init::Union{Nothing, Point2f} # save click position
        src::Union{Nothing, Point2f}  # save src vertex position
        dst::Union{Nothing, Point2f}  # save dst vertex position
        EdgeDragAction() = new(nothing, nothing, nothing)
    end
    function (action::EdgeDragAction)(state, idx, event, axis)
        edge = collect(edges(g))[idx]
        if state == true
            if action.src===action.dst===action.init===nothing
                action.init = event.data
                action.src = p[:node_pos][][edge.src]
                action.dst = p[:node_pos][][edge.dst]
            end
            offset = event.data - action.init
            p[:node_pos][][edge.src] = action.src + offset
            p[:node_pos][][edge.dst] = action.dst + offset
            p[:node_pos][] = p[:node_pos][] # trigger change
        elseif state == false
            action.src = action.dst = action.init =  nothing
        end
    end
    edrag = EdgeDragHandler(EdgeDragAction())
    register_interaction!(ax, :edrag, edrag)
end

@testset "align_to_direction" begin
    using GraphMakie: align_to_dir
    using LinearAlgebra: normalize

    @test align_to_dir((:left, :center)) ≈ Point(1.0, 0)
    @test align_to_dir((:right, :center)) ≈ Point(-1.0, 0)
    @test align_to_dir((:center, :center)) ≈ Point(0.0, 0)

    @test align_to_dir((:left, :top)) ≈ normalize(Point(1.0, -1.0))
    @test align_to_dir((:right, :top)) ≈ normalize(Point(-1.0, -1))
    @test align_to_dir((:center, :top)) ≈ normalize(Point(0.0, -1))

    @test align_to_dir((:left, :bottom)) ≈ normalize(Point(1.0, 1.0))
    @test align_to_dir((:right, :bottom)) ≈ normalize(Point(-1.0, 1.0))
    @test align_to_dir((:center, :bottom)) ≈ normalize(Point(0.0, 1.0))

    # g = complete_graph(9)
    # nlabels_align = vec(collect(Iterators.product((:left,:center,:right),(:top,:center,:bottom))))
    # nlabels= repr.(nlabels_align)
    # graphplot(g; nlabels, nlabels_align, nlabels_distance=20)
end

@testset "test plot accessors" begin
    g = complete_graph(10)
    nlabels = repr.(1:nv(g))
    elabels = repr.(1:ne(g))
    fig, ax, p = graphplot(g)
    @test get_edge_plot(p) isa EdgePlot
    @test get_node_plot(p)[1][] == p[:node_pos][]
    @test get_arrow_plot(p).visible[] == false
    @test get_nlabel_plot(p) === nothing
    @test get_elabel_plot(p) === nothing

    fig, ax, p = graphplot(g; nlabels)
    @test get_nlabel_plot(p)[1][] == nlabels
    @test get_elabel_plot(p) === nothing

    fig, ax, p = graphplot(g; elabels)
    @test get_nlabel_plot(p) === nothing
    @test get_elabel_plot(p)[1][] == elabels

    fig, ax, p = graphplot(g; elabels, nlabels)
    @test get_nlabel_plot(p)[1][] == nlabels
    @test get_elabel_plot(p)[1][] == elabels
end

@testset "test Pointf" begin
    using GraphMakie: Pointf

    p = Point(0.0, 0.0)
    @test typeof(Pointf(p)) == Point2f
    @test Pointf(p) == Point2f(p)

    p = Point(1, 0)
    @test typeof(Pointf(p)) == Point2f
    @test Pointf(p) == Point2f(p)

    p = Point(0.0, 1.0, 2.0)
    @test typeof(Pointf(p)) == Point3f
    @test Pointf(p) == Point3f(p)

    @test Pointf(0.0, 0.0, 0.0) isa Point3f
    @test Pointf(1.0, 1.0) isa Point2f

    @test Pointf((0.0, 0.0, 0.0 )) isa Point3f
    @test Pointf((1.0, 1.0)) isa  Point2f
end


text:
module MakieTeX
using Makie, CairoMakie, Makie.MakieLayout
using Rsvg, Cairo
using Colors, LaTeXStrings

using Makie.GeometryBasics: origin, widths
using Makie.Observables

using Poppler_jll, Perl_jll, Ghostscript_jll

include("types.jl")
include("rendering.jl")
include("recipe.jl")
include("layoutable.jl")

export TeXDocument, CachedTeX
export dvi2svg, latex2dvi, rsvg2recordsurf, svg2rsvg
export teximg, teximg!, TeXImg
export LTeX

end # document

ext: jl


text:
"""
    ObservationContinuousTime

Parametric composite type to store a single birth and death process realization,
observed continuously over time.
"""
struct ObservationContinuousTime{
  F <: AbstractFloat
}
  tot_births::Int
  tot_deaths::Int
  integrated_jump::F

  # sum of log population sizes (ancillary statistic)
  sum_log_n::Float64

  len::Int
  waiting_time::Vector{F}
  initial_population_size::Int
  increment::Vector{Int}
end

"""
    observation_continuous_time(t, x, observation_time)

Construct an object of type `ObservationContinuousTime` from a continuously
observed time series. Variable `x` is the vector of observed population sizes
while `t` is the vector of event times. `observation_time` is the total amount
of time used to observe the process, that is `t[end] ≦ observation_time`.
"""
function observation_continuous_time(
  t::Vector{F},
  x::Vector{Int},
  observation_time::F
) where {
  F <: AbstractFloat
}
  if length(x) != length(t)
    error("Arguments 't' and 'x' do not have the same length.")
  end

  if t[1] != 0
    error("The first observation time 't[1]' is not zero.")
  end

  if any(i -> i < 0, x)
    error("Vector 'x' cannot contain negative values.")
  end

  if any(i -> i < 0, t)
    error("Vector 't' cannot contain negative values.")
  end

  if !issorted(t)
    error("Vector 't' is not sorted in ascending order.")
  end

  if observation_time < t[end]
    error("'observation_time' is less than 't[end]'.")
  end

  tot_births = 0
  tot_deaths = 0
  integrated_jump = zero(F)
  len = length(t) - 1
  waiting_time = zeros(F, len)
  increment = zeros(Int, len)

  for s = 2:length(t)
    waiting_time[s - 1] = t[s] - t[s - 1]
    increment[s - 1] = x[s] - x[s - 1]

    if (abs(waiting_time[s - 1]) < floatmin(F)) && (increment[s - 1] != 0)
      error("Error at time point ", s, ". Observation time is zero but ",
            "'x[", s - 1, "]' and 'x[", s, "]' are different.")
    end

    if increment[s - 1] == -1
      tot_deaths += 1
    elseif increment[s - 1] == 1
      tot_births += 1
    else
      error("Error at time point ", s, ". Increment 'x[", s, "] - x[", s - 1,
            "]' is different from ±1.")
    end

    # add current waiting_time times the previous population size
    integrated_jump += x[s - 1] * waiting_time[s - 1]
  end

  # no events between t[end] and observation_time
  integrated_jump += x[s] * (observation_time - t[end])

  sum_log_n = sum(log, x[1:len])

  ObservationContinuousTime{F}(tot_births, tot_deaths, integrated_jump,
                               sum_log_n, len, waiting_time, x[1], increment)
end

"""
    ObservationDiscreteTimeEven

Parametric composite type to store birth and death process realizations,
observed at fixed time points with the same lag.
"""
struct ObservationDiscreteTimeEven{
  F <: AbstractFloat,
  I <: Integer
}
  # total number of samples
  n::Int

  # total number of time points
  k::Int

  # same time lag between each time points
  u::F

  # complete list of population sizes
  # column 'j' is j-th time series while row 'i' is time (i - 1) * u
  state::Matrix{I}
end

"""
    observation_discrete_time_even(u, x)

Construct an object of type `ObservationDiscreteTimeEven` from a discretely
observed time series. Variable `x` is a vector (or matrix if more than 1
observation) of observed population sizes while `u` is the fixed time lag after
which every observation was made.
"""
function observation_discrete_time_even(
  u::Real,
  x::Vector{I}
) where {
  I <: Integer
}
  if u <= 0
    error("Argument 'u' cannot be less than or equal to zero.")
  end

  if any(i -> i < 0, x)
    error("Vector 'x' cannot contain negative values")
  end

  k = length(x)
  F = typeof(float(u))

  ObservationDiscreteTimeEven{F, I}(1, k, F(u), reshape(x, k, 1))
end

function observation_discrete_time_even(
  u::Real,
  x::Matrix{I}
) where {
  I <: Integer
}
  if u <= 0
    error("Argument 'u' cannot be less than or equal to zero.")
  end

  if any(i -> i < 0, x)
    error("Matrix 'x' cannot contain negative values")
  end

  F = typeof(float(u))

  ObservationDiscreteTimeEven{F, I}(size(x, 2), size(x, 1), F(u), x)
end

"""
    ObservationDiscreteTimeUneven

Parametric composite type to store birth and death process realizations,
observed at discrete time points of varying length.
"""
struct ObservationDiscreteTimeUneven{
  F <: AbstractFloat,
  I <: Integer
}
  # complete list of waiting times
  waiting_time::Vector{F}

  # complete list of population sizes
  state::Vector{I}
end

"""
    observation_discrete_time_uneven(x, t)

Construct an object of type `ObservationDiscreteTimeUneven` from a discretely
observed time series. Variable `x` is a vector of observed population sizes
while `t` is the vector of event times.
"""
function observation_discrete_time_uneven(
  t::Vector{F},
  x::Vector{I}
) where {
  F <: AbstractFloat,
  I <: Integer
}
  if length(x) != length(t)
    error("Arguments 't' and 'x' do not have the same length.")
  end

  if t[1] != 0
    error("The first observation time 't[1]' is not zero.")
  end

  if any(i -> i < 0, x)
    error("Vector 'x' cannot contain negative values.")
  end

  if any(i -> i < 0, t)
    error("Vector 't' cannot contain negative values.")
  end

  if !issorted(t)
    error("Vector 't' is not sorted in ascending order.")
  end

  waiting_time = zeros(F, length(t) - 1)

  for s = 2:length(t)
    waiting_time[s - 1] = t[s] - t[s - 1]

    if (abs(waiting_time[s - 1]) < floatmin(F)) && (x[s - 1] != x[s])
      error("Error at time point ", s, ". Observation time is zero but ",
            "'x[", s - 1, "]' and 'x[", s, "]' are different.")
    end
  end

  ObservationDiscreteTimeUneven{F, I}(waiting_time, x)
end


text:
using Test, FiniteMesh

cd(@__DIR__)
mesh = Mesh("../mesh/Gmsh/t1.msh")

cells, nodes = read_mesh("../mesh/Gmsh/square.msh")
cellid = extract_cell(cells)
edgeNodes, edgeCells = mesh_face_connectivity_2D(cellid)
cellArea = mesh_cell_area_2D(nodes, cellid)

cells, nodes = read_mesh("../mesh/SU2/naca0012.su2")
add_group!(cells, "../mesh/SU2/naca0012.su2")

rg = [0.1, 0.3, 0.5, 0.7, 0.9]
x = hcat(rg, rg, rg)
rg = [0.1 0.3 0.5]
y = vcat(rg, rg, rg, rg, rg)

cells, points = unstructured_grid(cat(x, y, dims = 3), [0.0, 1.0], [0.0, 0.6])

_points = zeros(4, 2)
_points[2, :] = [1.0, 0.0]
_points[3, :] = [1.0, 1.0]
_points[4, :] = [0.0, 1.0]

cells = triangulate(_points)

FiniteMesh.unit_normal(rand(3), rand(3), rand(3))

cellNeighborsReg = regularize_cell_neighbor(mesh.cellid, mesh.cellNeighbors)
FiniteMesh.check_cell_neighbor(mesh.cellid, cellNeighborsReg)

cellFacesReg = regularize_cell_face(mesh.cellid, mesh.cellFaces, mesh.facePoints)
FiniteMesh.check_cell_face(mesh.cellid, cellFacesReg, mesh.facePoints)


text:
# The data folder layout
# Graph Instance (hash)
#    info.json
#    tensornetwork.json
#    SizeMax{2}/
#    ConfigsMax{2}/
#        size_10.dat
#        size_9.dat

# kwargs ∈ [optimizer, simplifier]
function foldername(basefolder::String, config::GraphProblemConfig; create, prefix="")
    # create a parameter dict
    # use the hash of dictionary as the folder name
    name = joinpath(basefolder, prefix * unique_string(config))

    # create folder
    if create && !isdir(name)
        mkpath(name)
    end
    if create
        # create info file
        js = joinpath(name, "info.toml")
        @info "writing configuration to: $(js)"
        # dump config to toml file
        to_toml(js, config)
    end
    return name
end

function save_code(folder, problem::GraphProblem)
    # write tensor network contraction pattern
    filename = joinpath(folder, "tensornetwork.json")
    @info "saving contraction tree to: $(filename)"
    GenericTensorNetworks.writejson(filename, problem.code)
end

function load_code(config::GraphProblemConfig, folder)
    # write tensor network contraction pattern
    filename = joinpath(folder, "tensornetwork.json")
    @info "loading contraction tree from: $(filename)"
    code = GenericTensorNetworks.readjson(filename)
    return problem_instance(config, code)
end

function saveconfigs(folderout, sizes, configs::AbstractVector{<:Union{ConfigEnumerator, SumProductTree}})
    for (s, c) in zip(sizes, configs)
        fname = joinpath(folderout, "size_$(s).dat")
        if c isa ConfigEnumerator
            save_configs(fname, c; format=:binary)
        else
            save_sumproduct(fname, c)
        end
    end
end

function loadconfigs(folderin, sizes; tree_storage::Bool, bitlength::Int)
    configs = []
    for s in sizes
        fname = joinpath(folderin, "size_$(s).dat")
        push!(configs, if tree_storage
                load_sumproduct(fname)
            else
                load_configs(fname; format=:binary, bitlength=bitlength)
            end
        )
    end
    return configs
end

function save_property(folder::String, property::GenericTensorNetworks.AbstractProperty, data)
    fd = joinpath(folder, "$(unique_string(property)).dat")
    @info "saving result to file/folder: $(fd)"
    if property isa SizeMax{Single} || property isa SizeMin{Single}
        writedlm(fd, data.n)
    elseif property isa SizeMax || property isa SizeMin
        writedlm(fd, getfield.(data.orders, :n))
    elseif property isa CountingAll
        writedlm(fd, data)
    elseif property isa CountingMax{Single} || property isa CountingMin{Single}
        writedlm(fd, [data.n, data.c]')
    elseif property isa CountingMax || property isa CountingMin
        K = length(data.coeffs)
        writedlm(fd, hcat(collect.(zip(data.maxorder-K+1:data.maxorder, data.coeffs))...)')
    elseif property isa GraphPolynomial
        writedlm(fd, data.coeffs)
    elseif property isa SingleConfigMax{Single} || property isa SingleConfigMin{Single}
        c = data.c.data
        writedlm(fd, [data.n, get_s(c), c...])
    elseif property isa SingleConfigMax || property isa SingleConfigMin
        writedlm(fd, vcat([[di.n, get_s(di.c.data), di.c.data...]' for di in data.orders]...))
    elseif property isa ConfigsMax || property isa ConfigsMin
        if property isa ConfigsMax{Single} || property isa ConfigsMin{Single}
            sizes = [data.n]
            bls = get_n(data.c)
            configs = [data.c]
        else
            sizes = collect(data.maxorder-length(data.coeffs)+1:data.maxorder)
            bls = get_n(data.coeffs[1])
            configs = [data.coeffs...]
        end
        !isdir(fd) && mkdir(fd)
        writedlm(joinpath(fd, "sizes.dat"), sizes)
        writedlm(joinpath(fd, "bitlength.dat"), bls)
        saveconfigs(fd, sizes, configs)
    elseif property isa ConfigsAll
        !isdir(fd) && mkdir(fd)
        writedlm(joinpath(fd, "bitlength.dat"), get_n(data))
        saveconfigs(fd, ["all"], [data])
    else
        error("unknown property: `$property`.")
    end

end
get_s(::StaticElementVector{N,S,C}) where {N,S,C} = S
get_n(::ConfigEnumerator{N,S,C}) where {N,S,C} = N
get_n(::SumProductTree{OnehotVec{N,F}}) where {N,F} = N
get_t(::ConfigsMax{K,B,T}) where {K,B,T} = T
get_t(::ConfigsMin{K,B,T}) where {K,B,T} = T
get_k(::ConfigsMax{K,B,T}) where {K,B,T} = K
get_k(::ConfigsMin{K,B,T}) where {K,B,T} = K
get_t(::ConfigsAll{T}) where {T} = T

function load_property(folder::String, property::GenericTensorNetworks.AbstractProperty; T=Float64)
    fd = joinpath(folder, "$(unique_string(property)).dat")
    if property isa SizeMax{Single} || property isa SizeMin{Single}
        return Tropical(readdlm(fd, T)[])
    elseif property isa SizeMax || property isa SizeMin
        orders = readdlm(fd, T)
        return ExtendedTropical{length(orders)}(Tropical.(vec(orders)))
    elseif property isa CountingAll
        return readdlm(fd, T)[]
    elseif property isa CountingMax{Single} || property isa CountingMin{Single}
        n, c = vec(readdlm(fd, T))
        return CountingTropical(n, c)
    elseif property isa CountingMax || property isa CountingMin
        data = readdlm(fd)
        return TruncatedPoly((data[:,2]...,), data[end,1])
    elseif property isa GraphPolynomial
        return Polynomial(vec(readdlm(fd)))
    elseif property isa SingleConfigMax{Single} || property isa SingleConfigMin{Single}
        data = readdlm(fd)
        n, s = data[1], Int(data[2])
        return CountingTropical(n, ConfigSampler(StaticElementVector(2^s, Int.(data[3:end]))))
    elseif property isa SingleConfigMax || property isa SingleConfigMin
        data = readdlm(fd)
        s = Int(data[1,2])
        orders = map(1:size(data, 1)) do i
            di = data[i,:]
            n = di[1]
            CountingTropical(n, ConfigSampler(StaticElementVector(2^s, Int.(di[3:end]))))
        end
        return ExtendedTropical{size(data, 1)}(orders)
    elseif property isa ConfigsMax || property isa ConfigsMin
        K = get_k(property)
        sizes = vec(readdlm(joinpath(fd, "sizes.dat")))[end-K+1:end]
        bitlength = Int(readdlm(joinpath(fd, "bitlength.dat"))[])
        data = loadconfigs(fd, sizes; tree_storage=get_t(property), bitlength=bitlength)
        if property isa ConfigsMax{Single} || property isa ConfigsMin{Single}
            return CountingTropical(T(sizes[]), data[])
        else
            return TruncatedPoly((data...,), T(sizes[end]))
        end
    elseif property isa ConfigsAll
        bitlength = Int(readdlm(joinpath(fd, "bitlength.dat"))[])
        data = loadconfigs(fd, ["all"]; tree_storage=get_t(property), bitlength=bitlength)
        return data[]
    else
        error("unknown property: `$property`.")
    end

end


text:
using LinearAlgebra: norm, dot, cross
using Statistics: mean
using MeshCore: IncRel, indextype, nshapes

"""
    selectnode(fens::FENodeSet; kwargs...)

Select nodes using some criterion.

# Arguments
- `v` = array of locations, one location per row
- `kwargs` = pairs of keyword argument/value

# Selection criteria:

## Distance from point

```
list = selectnode(fens.xyz, distance=1.0+0.1/2^nref, from=[0. 0.],
        inflate=tolerance);
```

## Distance from a plane

```
candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
```
The keyword `plane` defines the plane by its normal (the first two or
three numbers) and its distance from the origin (the last number). Nodes
are selected they lie on the plane,  or near the plane within the
distance `thickness` from the plane. The normal is assumed to be of unit
length, if it isn't apply as such, it will be normalized internally.

## Nearest-to a point

Find the node nearest to the location given.
```
nh = selectnode(fens, nearestto = [R+Ro/2, 0.0, 0.0] )
```
"""

function _box_outputlist!(outputlist::Vector{IT}, abox::BT, sdim::IT, v::VT) where {IT, BT, VT}
    # Helper functions
    @inline inrange(rangelo,rangehi,x) = (rangelo <= x <= rangehi)
    nn = 0
    for j in 1:length(v)
        matches = true
        for i in 1:sdim
            if !inrange(abox[2*i-1], abox[2*i], v[j][i])
                matches = false; break
            end
        end
        if matches
            nn = nn + 1; outputlist[nn] = j;
        end
    end
    return outputlist, nn
end

function _distance_outputlist!(outputlist::Vector{IT}, d, fromvalue, sdim::IT, v::VT) where {IT, BT, VT}
    # Helper functions
    nn = 0
    for j in 1:length(v)
        if norm(fromvalue-v[j]) < d
            nn = nn + 1; outputlist[nn] = j;
        end
    end
    return outputlist, nn
end

function _plane_outputlist!(outputlist::Vector{IT}, distance, normal, t, sdim::IT, v::VT) where {IT, BT, VT}
    # Helper functions
    nn = 0
    for j in 1:length(v)
        ad = dot(v[j], normal);
        if abs(distance-ad)<t
            nn = nn + 1; outputlist[nn] = j;
        end
    end
    return outputlist, nn
end

function _nearestto_outputlist!(outputlist::Vector{IT}, nearestto, sdim::IT, v::VT) where {IT, BT, VT}
    distances =  fill(0.0, length(v));
    for j in 1:length(v)
        distances[j] = norm(nearestto-v[j])
    end
    Mv,j = findmin(distances)
    return [j], 1
end

"""
    vselect(v::VT; kwargs...) where {VT}

Select locations (vertices) based on some criterion.

`VT` is an abstract array that returns the coordinates of a vertex given its
number.

## Selection criteria

### box
```
nLx = vselect(v, box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)
```

The keyword 'inflate' may be used to increase or decrease the extent of
the box (or the distance) to make sure some nodes which would be on the
boundary are either excluded or included.

### distance
```
list = vselect(v, distance=1.0+0.1/2^nref, from=[0. 0.], inflate=tolerance);
```

### plane
```
candidates = vselect(v, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
```
The keyword `plane` defines the plane by its normal (the first two or
three numbers) and its distance from the origin (the last number). Nodes
are selected they lie on the plane,  or near the plane within the
distance `thickness` from the plane. The normal is assumed to be of unit
length, if it isn't provided as such, it will be normalized internally.

### nearestto
Find the node nearest to the location given.
```
nh = vselect(v, nearestto = [R+Ro/2, 0.0, 0.0])
```

# Returns
The list of vertices that match the search criterion.
"""
function vselect(v::VT; kwargs...) where {VT}
    # Extract arguments
    box = nothing; distance = nothing; from = nothing; plane  =  nothing;
    thickness = nothing; nearestto = nothing; inflate = 0.0;
    for apair in pairs(kwargs)
        sy, val = apair
        if sy == :box
            box = val
        elseif sy == :distance
            distance = val
        elseif sy == :from
            from = val
        elseif sy == :plane
            plane = val
        elseif sy == :thickness
            thickness = val
        elseif sy == :nearestto
            nearestto = val
        elseif sy == :inflate
            inflate = val
        end
    end

    # Did we get an inflate value
    inflatevalue = 0.0;
    if inflate != nothing
        inflatevalue = Float64(inflate);
    end

    # Initialize the output list
    outputlist = zeros(Int64, length(v)); nn = 0;

    sdim = length(v[1])
    # Process the different options
    if box != nothing
        dim = Int64(round(length(box)/2.));
        @assert dim == sdim "Dimension of box not matched to dimension of array of vertices"
        abox = vec(box)
        inflatebox!(abox, inflatevalue)
        outputlist, nn = _box_outputlist!(outputlist, abox, sdim, v)
    elseif distance != nothing
        fromvalue = fill(0.0, sdim);
        if from != nothing
            fromvalue = from;
        end
        d = distance+inflatevalue;
        outputlist, nn = _distance_outputlist!(outputlist, d, fromvalue, sdim, v)
    elseif plane != nothing
        normal = plane[1:end-1];
        normal = vec(normal/norm(normal));
        thicknessvalue = 0.0;
        if thickness != nothing
            thicknessvalue = thickness;
        end
        t = thicknessvalue+inflatevalue;
        distance = plane[end];
        outputlist, nn = _plane_outputlist!(outputlist, distance, normal, t, sdim, v)
    elseif nearestto != nothing
        nearestto = vec(nearestto);
        outputlist, nn = _nearestto_outputlist!(outputlist, nearestto, sdim, v)
    end
    if (nn==0)
        outputlist = Int64[];# nothing matched
    else
        outputlist = outputlist[1:nn];
    end
    return outputlist
end

"""
    connectedv(fes::AbstractFESet)

Extract the node numbers of the nodes connected by given finite elements.

Extract the list of unique node numbers for the nodes that are connected
by the finite element set `fes`. Note that it is assumed that all the
FEs are of the same type (the same number of connected nodes by each
cell).
"""
function connectedv(ir::IncRel)
    vl = fill(zero(indextype(ir)), 0)
    for i in 1:nshapes(ir.left)
        append!(vl, ir[i])
    end
    return unique(vl);
end


text:
# Gradient regression
function get_gain(loss::S, ∑δ::AbstractVector{T}, ∑δ²::AbstractVector{T}, ∑𝑤::T, λ::T) where {S <: Union{GradientRegression,GaussianRegression}, T <: AbstractFloat}
    gain = sum((∑δ .^ 2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end

#####################
# linear
#####################
function kernel_linear_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        @inbounds δ[i] = 2 * (p[i] - t[i]) * 𝑤[i]
    end
    return
end

function kernel_linear_δ²!(δ²::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        @inbounds δ²[i] = 2 * 𝑤[i]
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Linear, δ::CuMatrix{T}, δ²::CuMatrix{T}, p::CuMatrix{T}, t::CuVector{T}, 𝑤::CuVector{T}; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_i = min(MAX_THREADS, length(t))
    threads = (thread_i)
    blocks = ceil.(Int, (length(t)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_linear_δ!(δ, p, t, 𝑤)
    @cuda blocks=blocks threads=threads kernel_linear_δ²!(δ², p, t, 𝑤)
    return
end

#####################
# Logistic
#####################
function kernel_logistic_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        @inbounds δ[i] = (sigmoid(p[i]) * (1 - t[i]) - (1 - sigmoid(p[i])) * t[i]) * 𝑤[i]
    end
    return
end

function kernel_logistic_δ²!(δ²::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        @inbounds δ²[i] = sigmoid(p[i]) * (1 - sigmoid(p[i])) * 𝑤[i]
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Logistic, δ::CuMatrix{T}, δ²::CuMatrix{T}, p::CuMatrix{T}, t::CuVector{T}, 𝑤::CuVector{T}; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_i = min(MAX_THREADS, length(t))
    threads = (thread_i)
    blocks = ceil.(Int, (length(t)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_logistic_δ!(δ, p, t, 𝑤)
    @cuda blocks=blocks threads=threads kernel_logistic_δ²!(δ², p, t, 𝑤)
    return
end

################################################################################
# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ)
################################################################################
function kernel_gauss_δ!(δ::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        δ[i,1] = (p[i,1] - t[i]) / max(Cfloat(1e-5), exp(2f0 * p[i,2])) * 𝑤[i]
        δ[i,2] = (1f0 - (p[i,1] - t[i])^2 / max(Cfloat(1e-5), exp(2f0 * p[i,2]))) * 𝑤[i]
    end
    return
end

function kernel_gauss_δ²!(δ²::CuDeviceMatrix{T}, p::CuDeviceMatrix{T}, t::CuDeviceVector{T}, 𝑤::CuDeviceVector{T}) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(t)
        δ²[i,1] = 𝑤[i] / max(Cfloat(1e-5), exp(2 * p[i,2]))
        δ²[i,2] = 2 * 𝑤[i] / max(Cfloat(1e-5), exp(2 * p[i,2])) * (p[i,1] - t[i])^2
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_grads_gpu!(loss::Gaussian, δ::CuMatrix{T}, δ²::CuMatrix{T}, p::CuMatrix{T}, t::CuVector{T}, 𝑤::CuVector{T}; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_i = min(MAX_THREADS, length(t))
    threads = (thread_i)
    blocks = ceil.(Int, (length(t)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_gauss_δ!(δ, p, t, 𝑤)
    @cuda blocks=blocks threads=threads kernel_gauss_δ²!(δ², p, t, 𝑤)
    return
end

# # GaussianRegression
# function get_gain(loss::S, ∑δ::AbstractVector{T}, ∑δ²::AbstractVector{T}, ∑𝑤::T, λ::T) where {S <: GaussianRegression, T <: AbstractFloat}
#     gain = sum((∑δ .^ 2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
#     return gain
# end


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]

#######################
print("Testing ACQR_em_starting: ")
wi = collect(range(0.1, length=3, stop=3))
Ai,Ci,Qi,Ri,Si,m1i,P1i = ACQR_em_starting(y,wi,1,true)
if sum(Ai.^2) > 1e-6
    println("OK")
end


text:
#=
bisect_right:
- Julia version: 1.0
- Author: qz
- Date: 2018-08-25
=#

using Test

using BisectPy: bisect_right

@test bisect_right([1, 2, 3, 4, 5], 3.5) == 4

@test bisect_right([1, 2, 3, 4, 5], 2) == 3

@test bisect_right([1, 2, 3, 3, 3, 5], 3) == 6


text:
@testset "GeneralParserMethods" begin

    @testset "Test input lengths" begin
        # NOTE: _pre_read... functions should call _check_seq_len
        # if _convert_to_matrix_residues is used

        ids = ["A", "B", "C"]
        seqs_equal = ["MSRSKRDNEIGDSTF", "MSRSKRDNNFGDSTF", "MSRSFYSVEIGDSTF"]
        seqs_diffs = ["MSRSKRDNEIGDSTF", "MSRSKRDNNFGDSTF", "MSRSFYSVEIGD"]
        seqs_less = ["MSRSKRDNEIGDSTF", "MSRSKRDNNFGDSTF"]

        @test MSA._check_seq_len(ids, seqs_equal) === nothing
        @test_throws ErrorException  MSA._check_seq_len(ids, seqs_diffs)
        @test_throws ErrorException  MSA._check_seq_len(ids, seqs_less)
    end

    @testset "To parse MSA & mapping" begin

        sequences = ["ADEIMSY","RCGLFTV","NQHKPW-"]
        matrixres = reshape(reinterpret(Residue, collect(1:21)), (3,7))
        msa, map = MSA._to_msa_mapping(sequences)

        @test getarray(msa) == matrixres
        @test map == ["1,2,3,4,5,6,7", "1,2,3,4,5,6,7", "1,2,3,4,5,6,"]
        @test sequencenames(msa) == ["1", "2", "3"]
        # MSA constructors adds dimension names
        # @test dimnames(msa) == [:Seq, :Col]

        msa, map = MSA._to_msa_mapping(sequences,["a/11-17","b/11-17","c/11-16"])

        @test getarray(msa) == matrixres
        @test map == ["11,12,13,14,15,16,17", "11,12,13,14,15,16,17", "11,12,13,14,15,16,"]
        @test sequencenames(msa) == ["a/11-17", "b/11-17", "c/11-16"]
        # MSA constructors adds dimension names
        # @test dimnames(msa) == [:Seq, :Col]

        @test_throws ErrorException MSA._to_msa_mapping(sequences,["a/1-7","b/1-7","c/1-7"])
        @test_throws ErrorException MSA._to_msa_mapping(["AD","EI","MSY"])
    end

    @testset "Delete full gap columns" begin

        M = reshape(reinterpret(Residue,collect(1:21)),(3,7))
        M[:,[2,4,6]] .= GAP
        msa = MultipleSequenceAlignment(M)
        named = namedmatrix(msa)
        annotated_msa = AnnotatedMultipleSequenceAlignment(M)

        d_M = deletefullgapcolumns(M)
        d_msa = deletefullgapcolumns(msa)
        d_named = deletefullgapcolumns(named)
        d_annotated_msa = deletefullgapcolumns(annotated_msa)

        @test size(d_M) == (3,4)

        for object in [d_msa, d_named, d_annotated_msa]
            @test size(object) == (3,4)
            @test getcolumnmapping(object) == [1,3,5,7]
        end

        d_msa = deletefullgapcolumns!(msa)
        d_annotated_msa = deletefullgapcolumns!(annotated_msa)

        @test d_msa == msa
        @test d_annotated_msa == annotated_msa
    end
end


text:
# Determine the measure
"""
    azdual_dict(ap::ApproximationProblem; options...)

The dual that is used to create a AZ `Z` matrix.
"""
azdual_dict(samplingstyle::DiscreteStyle, ap::ApproximationProblem; options...) =
    azdual_dict(samplingstyle, ap, discretemeasure(samplingstyle, ap; options...); options...)
azdual_dict(samplingstyle::SamplingStyle, ap::ApproximationProblem; options...) =
    azdual_dict(samplingstyle, ap, measure(samplingstyle, ap; options...); options...)

azdual_dict(ss::ProductSamplingStyle, ap::ApproximationProblem; options...) =
    TensorProductDict( elements(azdual_dict, ss, ap; options...)... )

azdual_dict(samplingstyle::SamplingStyle, ap::Platform, param, measure; options...) =
    dualdictionary(ap, param, measure; samplingstyle=samplingstyle, options...)
default_azdual_dict(::SamplingStyle, dict::Dictionary, measure; options...) =
    (@warn "azdualdict on dictionary is deprecated";BasisFunctions.default_gramdual(dict, measure; options...))


text:
# Wrapper module for scikit-learn machine learners.
module ScikitLearnWrapper

importall Orchestra.Types
importall Orchestra.Util

using PyCall
@pyimport sklearn.ensemble as ENS
@pyimport sklearn.linear_model as LM
@pyimport sklearn.lda as LDA_M
@pyimport sklearn.qda as QDA_M
@pyimport sklearn.neighbors as NN
@pyimport sklearn.svm as SVM
@pyimport sklearn.tree as TREE

export SKLLearner,
       fit!,
       transform!

# Available scikit-learn learners.
learner_dict = {
  "AdaBoostClassifier" => ENS.AdaBoostClassifier,
  "BaggingClassifier" => ENS.BaggingClassifier,
  "ExtraTreesClassifier" => ENS.ExtraTreesClassifier,
  "GradientBoostingClassifier" => ENS.GradientBoostingClassifier,
  "RandomForestClassifier" => ENS.RandomForestClassifier,
  "LDA" => LDA_M.LDA,
  "LogisticRegression" => LM.LogisticRegression,
  "PassiveAggressiveClassifier" => LM.PassiveAggressiveClassifier,
  "RidgeClassifier" => LM.RidgeClassifier,
  "RidgeClassifierCV" => LM.RidgeClassifierCV,
  "SGDClassifier" => LM.SGDClassifier,
  "KNeighborsClassifier" => NN.KNeighborsClassifier,
  "RadiusNeighborsClassifier" => NN.RadiusNeighborsClassifier,
  "NearestCentroid" => NN.NearestCentroid,
  "QDA" => QDA_M.QDA,
  "SVC" => SVM.SVC,
  "LinearSVC" => SVM.LinearSVC,
  "NuSVC" => SVM.NuSVC,
  "DecisionTreeClassifier" => TREE.DecisionTreeClassifier
}

# Wrapper for scikit-learn that provides access to most learners.
#
# Options for the specific scikit-learn learner is to be passed
# in `options[:impl_options]` dictionary.
#
# Available learners:
#
#   - "AdaBoostClassifier"
#   - "BaggingClassifier"
#   - "ExtraTreesClassifier"
#   - "GradientBoostingClassifier"
#   - "RandomForestClassifier"
#   - "LDA"
#   - "LogisticRegression"
#   - "PassiveAggressiveClassifier"
#   - "RidgeClassifier"
#   - "RidgeClassifierCV"
#   - "SGDClassifier"
#   - "KNeighborsClassifier"
#   - "RadiusNeighborsClassifier"
#   - "NearestCentroid"
#   - "QDA"
#   - "SVC"
#   - "LinearSVC"
#   - "NuSVC"
#   - "DecisionTreeClassifier"
#
type SKLLearner <: Learner
  model
  options

  function SKLLearner(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      :learner => "LinearSVC",
      # Options specific to this implementation.
      :impl_options => Dict()
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(sklw::SKLLearner, instances::Matrix, labels::Vector)
  impl_options = copy(sklw.options[:impl_options])
  learner = sklw.options[:learner]
  py_learner = learner_dict[learner]

  # Assign Orchestra-specific defaults if required
  if learner == "RadiusNeighborsClassifier"
    if get(impl_options, :outlier_label, nothing) == nothing
      impl_options[:outlier_label] = labels[rand(1:size(labels, 1))]
    end
  end

  # Train
  sklw.model = py_learner(;impl_options...)
  sklw.model[:fit](instances, labels)
end

function transform!(sklw::SKLLearner, instances::Matrix)
  return collect(sklw.model[:predict](instances))
end

end # module


text:
# $name is the parameter
function say_hello(name)
    println("Hi $name, it's great to see you!")
end
# function call
say_hello("Sir")
say_hello("Dude")

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
################################################################################
#
#   module LatPhysLuttingerTiszaPlottingPyPlot
#   -> PyPlot
#   -> LatticePhysics_Base
#   -> LatticePhysics_ReciprocalSpace
#   -> LatticePhysics_LuttingerTisza_Calculations
#   -> LinearAlgebra
#
#   --> PLOTTING OF LT BANDSTRUCTURES AND EIGENVALUES
#           - in real space (pure eigenvalues)
#           - in momentum space (along cardinal directions)
#           - in momentum space (along reciprocal lattice directions)
#           - in momentum space (along path)
#
#   --> PLOTTING OF LT GROUNDSTATE MANIFOLDS
#           - plotting of manifolds in 2D
#           - plotting of manifolds in 3D
#
################################################################################

# start of module
module LatPhysLuttingerTiszaPlottingPyPlot

# used libraries
using Reexport
using LatPhysBase
using LatPhysReciprocal
using LatPhysLuttingerTisza
using PyPlot
using LinearAlgebra

# reexport reciprocal pyplot plotting
@reexport using LatPhysReciprocalPlottingPyPlot

# import plotting
#include("LT_groundstate_plotting.jl")

# end of module
end


text:
export step!, run!, log_gen, save_gen

using Logging
using Statistics

function step!(e::Cambrian.Evolution)
    e.gen += 1
    if e.gen > 1
        e.populate(e)
    end
    e.generation(e)
    e.evaluate(e)
    if ((e.cfg["log_gen"] > 0) && mod(e.gen, e.cfg["log_gen"]) == 0)
        BERL.log_gen(e)
    end
    if ((e.cfg["save_gen"] > 0) && mod(e.gen, e.cfg["save_gen"]) == 0)
        save_gen(e)
    end
end

function run!(e::Cambrian.Evolution)
    for i in (e.gen+1):e.cfg["n_gen"]
        BERL.step!(e)
    end
end

function log_gen(e::Cambrian.Evolution)
    maxs = map(i->maximum(i.fitness), e.population)
    evals = e.gen * length(e.population)
    Cambrian.with_logger(e.log) do
        @info Formatting.format("{1} {2:04d} {3} {4:e} {5:e} {6:e} {7:e}",
                                e.id, e.gen, e.text, maximum(maxs),
                                mean(maxs), std(maxs), evals)
    end
    flush(e.log.stream)
end

function save_gen(e::Cambrian.Evolution)
    # save the entire population
    path = Formatting.format("gens/{1}/{2:04d}", e.id, e.gen)
    mkpath(path)
    sort!(e.population)
    for i in eachindex(e.population)
        f = open(Formatting.format("{1}/{2:04d}.dna", path, i), "w+")
        write(f, string(e.population[i]))
        close(f)
    end
end


text:
@testset "Promoting constructors" begin
    @testset "LorentzVector" begin
        @test Vec4(0, 1, 2, 3) == Vec4(0., 1., 2., 3.)
        @test typeof(Vec4(0, 1, 2, 3)) == LorentzVector{Float64}
        @test Vec4(0., 1, 2, 3) == Vec4(0., 1., 2., 3.)
        @test typeof(Vec4(0., 1, 2, 3)) == LorentzVector{Float64}
        @test typeof(Vec4(0f0, 1, 2, 3)) == LorentzVector{Float32}
        @test typeof(Vec4(0, 1f0, 2, 3)) == LorentzVector{Float32}
        @test typeof(Vec4(0, 1, 2f0, 3)) == LorentzVector{Float32}
        @test typeof(Vec4(0, 1, 2, 3f0)) == LorentzVector{Float32}
        @test typeof(Vec4(0f0, 1., 2f0, 3f0)) == LorentzVector{Float64}
        @test Vec4(0//1, 1, 2, 3) == Vec4(0., 1., 2., 3.)
        @test Vec4(π, π, π, π) ≈ Vec4(float(π), float(π), float(π), float(π))
    end;

    @testset "SpatialVector" begin
        @test Vec3(0, 1, 2) == Vec3(0., 1., 2.)
        @test typeof(Vec3(0, 1, 2)) == SpatialVector{Float64}
        @test Vec3(0., 1, 2) == Vec3(0., 1., 2.)
        @test typeof(Vec3(0., 1, 2)) == SpatialVector{Float64}
        @test typeof(Vec3(0f0, 1, 2)) == SpatialVector{Float32}
        @test typeof(Vec3(0, 1f0, 2)) == SpatialVector{Float32}
        @test typeof(Vec3(0, 1, 2f0)) == SpatialVector{Float32}
        @test typeof(Vec3(0f0, 1., 2f0)) == SpatialVector{Float64}
        @test Vec3(0//1, 1, 2) == Vec3(0., 1., 2.)
        @test Vec3(π, π, π) ≈ Vec3(float(π), float(π), float(π))
    end;

    @testset "Misc. constructors" begin
        x = Vec3(1, 2, 3)
        @test Vec4(0.5, x) == Vec4(0.5, 1., 2., 3.)
        u = Vec4(0, 1, 2, 3)
        @test Vec3(u) == Vec3(1., 2., 3.)
    end;

    @testset "Zero constructor" begin
        x = zero(Vec3)
        @test x === Vec3(0.0, 0.0, 0.0)
        x64 = zero(Vec3{Float64})
        @test x === x64
        x32 = zero(Vec3{Float32})
        @test x32 === Vec3{Float32}(0f0, 0f0, 0f0)
        y = zero(Vec4)
        @test y === Vec4(0.0, 0.0, 0.0, 0.0)
        y64 = zero(Vec4{Float64})
        @test y === y64
        y32 = zero(Vec4{Float32})
        @test y32 === Vec4(0f0, 0f0, 0f0, 0f0)
        @test zeros(Vec3, 2) == [x, x]
        @test zeros(Vec4, 2) == [y, y]
    end
end;


text:
module Grids

export AbstractGrid
export RectilinearGrid
export LatLonGrid
export LoadLatLonGrid
export Periodic, Bounded, short_show
export replace_grid_storage
export ΔxC, ΔyC, ΔzC, ΔxF, ΔyF, ΔzF, Ax, Ay, Az, volume

using CUDA
using Adapt

using PlanktonIndividuals.Architectures

"""
    AbstractGrid{TX, TY, TZ}
Abstract type for grids with elements of type `Float64` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{TX, TY, YZ} end

"""
    AbstractTopology
Abstract type for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic
Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded
Grid topology for bounded dimensions.
"""
struct Bounded <: AbstractTopology end

import Base: show

include("rectilinear_grid.jl")
include("lat_lon_grid.jl")
include("utils.jl")
include("areas_volumes_spacings.jl")

end


text:
module TestFilterCircuit

using TinyModia

using DifferentialEquations
using ModiaPlot
using Test

setLogMerge(false)

include("../models/Electric.jl")

FilterCircuit = Model(
    R = Resistor | Map(R=0.5u"Ω"),
    C = Capacitor | Map(C=2.0u"F", v=Var(init=0.1u"V")),
    V = ConstantVoltage | Map(V=10.0u"V"),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

filterCircuit = @instantiateModel(FilterCircuit)

simulate!(filterCircuit, Tsit5(), stopTime = 10, merge = Map(R = Map(R = 5u"Ω"), C = Map(v = 3.0u"V")),
          logParameters = true, logStates = true, requiredFinalStates = [7.424843902110655])

# Test access functions
@testset "Test variable access functions" begin
    @test ModiaPlot.getNames(filterCircuit) == ["C.C", "C.class", "C.i", "C.n.i", "C.n.v", "C.p.i", "C.p.v", "C.v",
"C.v", "R.R", "R.class", "R.i", "R.n.i", "R.n.v", "R.p.i", "R.p.v", "R.v", "V.V", "V.i", "V.n.i", "V.n.v", "V.p.i", "V.p.v", "V.v", "class", "der(C.v)", "ground.p.i", "ground.p.v", "time"]
    @test ModiaPlot.hasSignal(filterCircuit, "R.v")
    @test ModiaPlot.hasSignal(filterCircuit, "C.n.v")
    @test ModiaPlot.hasSignal(filterCircuit, "R.R")
    @test ModiaPlot.hasSignal(filterCircuit, "ground.p.i")
    @test ModiaPlot.hasSignal(filterCircuit, "R.p.vv") == false
    @test isapprox(get_lastValue(filterCircuit, "R.v") , 2.5751560978893453u"V" )
    @test isapprox(get_lastValue(filterCircuit, "C.n.v"), 0.0u"V")
    @test isapprox(get_lastValue(filterCircuit, "R.R")  , 5.0u"Ω")
    @test isapprox(get_lastValue(filterCircuit, "ground.p.i"), 0.0)
end

# Test plotting of variables, zero variables, parameters
plot(filterCircuit, [("R.v", "C.v"), ("R.R", "ground.p.i")], figure=1)

# Simulate with lower precision
filterCircuitLow = @instantiateModel(FilterCircuit, FloatType = Float32)
simulate!(filterCircuitLow, RK4(), adaptive=false, stopTime=10.0, interval=0.01,
          merge = Map(R = Map(R = 5u"Ω"), C = Map(v = 3.0u"V")),
          requiredFinalStates = Float32[7.4248414])
plot(filterCircuitLow, [("R.v", "C.v"), ("R.R", "ground.p.i")], figure=2)

end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# "classic" forward - just for testing
function Baum_forward(a,A,y)
  Nt, Ns = size(y,1), size(A,1)
  alpha = zeros(Nt,Ns)
  for j = 1:Ns
    alpha[1,j] = a[j]*y[1,j]
  end

  for t = 2:Nt
    for j = 1:Ns
      for i = 1:Ns
        alpha[t,j] = alpha[t,j] + alpha[t-1,i]*A[i,j]*y[t,j]
      end
    end
  end
  return alpha
end

function Baum_forward_matrix(a,A,y)
  Nt, Ns = size(y,1), size(A,1)
  alpha = zeros(Nt,Ns)
  alpha[1,:] .= Diagonal(y[1,:])*a
  for t = 2:Nt
    alpha[t,:] .= Diagonal(y[t,:])*A'*alpha[t-1,:]
  end
  return alpha
end

# "classic" backwards - just for testing
function Baum_backward(A,y)
  Nt, Ns = size(y,1), size(A,1)
  beta = zeros(Nt,Ns)
  beta[end,:] .= 1
  for t = Nt-1:-1:1
    for j = 1:Ns
      for k = 1:Ns
        beta[t,j] += beta[t+1,k]*A[j,k]*y[t+1,k]
      end
    end
  end
  return beta
end

function Baum_backward_matrix(A,y)
  Nt, Ns = size(y,1), size(A,1)
  beta = zeros(Nt,Ns)
  beta[Nt,:] .= ones(Ns)
  for t = Nt-1:-1:1
    beta[t,:] .= A*Diagonal(y[t+1,:])*beta[t+1,:]
  end
  return beta
end

# finite difference gradient
function gradient_fd(f, y)
  fy = f(y)
	δ = sqrt(eps())
  grad = similar(y)
  for i in eachindex(y)
		yδ = copy(y)
		yδ[i] = yδ[i]+δ
    grad[i] = ( f(yδ) - fy )./δ
	end
	return grad
end


text:
@checked struct PassRegistry
    ref::API.LLVMPassRegistryRef
end
reftype(::Type{PassRegistry}) = API.LLVMPassRegistryRef

export GlobalPassRegistry

GlobalPassRegistry() = PassRegistry(API.LLVMGetGlobalPassRegistry())


text:
@deprecate simple_hessian hessian_backback
@deprecate hessian_repeat hessian_backback
@deprecate ngradient gradient_numeric
@deprecate nhessian hessian_numeric
@deprecate NEG Base.:-
@deprecate ipush! PUSH!
@deprecate ipop! POP!


text:
using Test, PolynomialBases
using LinearAlgebra
import StaticArrays

ufunc₁(x) = sinpi(x)
ufunc₂(x) = exp(sinpi(x))
ufunc₃(x) = cospi(x^2)^5

ufunc_svector(x) = StaticArrays.SVector(ufunc₁(x), ufunc₃(x))

function tolerance(p, ufunc::typeof(ufunc₁))
    if p <= 4
        1.
    elseif p <= 6
        0.03
    elseif p <= 8
        2.e-3
    elseif p <= 10
        9.5e-4
    else
        2.e-7
    end
end
function tolerance(p, ufunc::typeof(ufunc₂))
    if p <= 4
        2.
    elseif p <= 6
        0.7
    elseif p <= 7
        0.4
    elseif p <= 8
        0.2
    elseif p <= 10
        0.06
    else
        6.e-3
    end
end
function tolerance(p, ufunc::typeof(ufunc₃))
    if p <= 2
        2.
    elseif p <= 6
        0.75
    elseif p <= 8
        0.7
    elseif p <= 10
        0.4
    else
        0.06
    end
end

# some regression tests
for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    for ufc in (ufunc₁, ufunc₂, ufunc₃), p in 3:10, T in (Float32, Float64)
        basis = try
            basis_type(p, T)
        catch m
            isa(m, DomainError) && continue
            throw(m)
        end
        u = ufc.(basis.nodes)
        u2 = compute_coefficients(ufc, basis)
        @test maximum(abs, u-u2) < 10*eps(T)
        xplot = range(-1, stop=1, length=100)
        uplot = interpolate(xplot, u, basis)
        @test norm(ufc.(xplot) - uplot, Inf) < tolerance(p, ufc)
        xplot2, uplot2 = evaluate_coefficients(u, basis, length(xplot))
        @test maximum(abs, xplot-xplot2) < 10*eps(T)
        @test maximum(abs, uplot-uplot2) < 10*eps(T)
    end
end

# compare direct interpolation and matrix
for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    for ufc in (ufunc₁, ufunc₂, ufunc₃), p in 3:10, T in (Float32, Float64)
        basis = try
            basis_type(p, T)
        catch m
            isa(m, DomainError) && continue
            throw(m)
        end
        u = ufc.(basis.nodes)
        xplot = range(-1, stop=1, length=100)
        u1 = interpolate(xplot, u, basis)
        u2 = interpolation_matrix(xplot, basis) * u
        @test u1 ≈ u2 atol=1.e-14
    end
end

# compare direct interpolation and interpolation vectors for GaussLegendre
for ufc in (ufunc₁, ufunc₂, ufunc₃), p in 1:10, T in (Float32, Float64)
    basis = GaussLegendre(p, T)
    u = ufc.(basis.nodes)
    x = T[-1, 1]
    u1 = interpolate(x, u, basis)
    u2 = similar(u1); u2[1] = dot(basis.interp_left, u); u2[2] = dot(basis.interp_right, u)
    @test u1 ≈ u2 atol=5*eps(T)
end

# change of bases
for ufc in (ufunc₁, ufunc₂, ufunc₃), p in 3:10, T in (Float32, Float64)
    basis1 = LobattoLegendre(p, T)
    basis2 = GaussLegendre(p, T)

    u1 = ufc.(basis1.nodes)
    u2 = ufc.(basis2.nodes)

    u12 = change_basis(basis1, u2, basis2)
    u21 = change_basis(basis2, u1, basis1)

    xplot = range(-1, stop=1, length=100)

    @test norm(interpolate(xplot,u1,basis1) - interpolate(xplot,u12,basis1), Inf) < tolerance(p, ufc)
    @test norm(interpolate(xplot,u2,basis2) - interpolate(xplot,u21,basis2), Inf) < tolerance(p, ufc)
end

# interpolate SVector
for T in (Float32,Float64)
    basis = GaussLegendre(5, T)
    u = compute_coefficients(ufunc_svector, basis)
    @inferred evaluate_coefficients(u, basis)
end


text:
module OCIM1
#=
This module serves to load the OCIM1 matrix and grid.
These are loaded from the public, persistant URL in FigShare.
The Julia BSON format version of the OCIM1 was created with the code
in the GitHub repository https://github.com/briochemc/OceanCirculations.
=#

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using BSON                  # For saving circulation as BSON format
using Unitful, UnitfulAstro # for units
using Reexport
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

# Create registry entry for OCIM in JLD2 format
function register_OCIM1()
    register(
        DataDep(
            "AIBECS_OCIM1",
            """
            References:
            - DeVries, T., 2014: The oceanic anthropogenic CO2 sink: Storage, air‐sea fluxes, and transports over the industrial era, Global Biogeochem. Cycles, 28, 631–647, doi:10.1002/2013GB004739.
            - DeVries, T. and F. Primeau, 2011: Dynamically and Observationally Constrained Estimates of Water-Mass Distributions and Ages in the Global Ocean. J. Phys. Oceanogr., 41, 2381–2401, https://doi.org/10.1175/JPO-D-10-05011.1
            """,
            "https://files.figshare.com/17302703/OCIM1.bson",
            "45ebe60e166dcd1b5bdbdf787713b41cb1d5bb2bf4c61a57c8497367a68354e6",
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load

Returns wet3d, grd, and T (in that order) from FigShare repository.
"""
function load()
    print("Loading OCIM1")
    register_OCIM1()
    bson_file = @datadep_str string("AIBECS_OCIM1/", "OCIM1.bson")
    BSON.@load bson_file T grid wet3D
    println(" ✔")
    @info """You are about to use OCIM1 model.
          If you use it for research, please cite:

          - DeVries, T., 2014: The oceanic anthropogenic CO2 sink: Storage, air‐sea fluxes, and transports over the industrial era, Global Biogeochem. Cycles, 28, 631–647, doi:10.1002/2013GB004739.
          - DeVries, T. and F. Primeau, 2011: Dynamically and Observationally Constrained Estimates of Water-Mass Distributions and Ages in the Global Ocean. J. Phys. Oceanogr., 41, 2381–2401, doi:10.1175/JPO-D-10-05011.1

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "DeVries_Primeau_2011" and "DeVries_2014" keys.)
          """
    return wet3D, grid, ustrip.(T)
end

end # end module

export OCIM1


text:
#####
##### Parsing utilities
#####

function Base.tryparse(::Type{PatientID}, raw::AbstractString)
    s = split(raw, ' ', keepempty=false)
    length(s) == 4 || return
    code_raw, sex_raw, dob_raw, name_raw = s
    length(sex_raw) == 1 || return
    code = edf_unknown(code_raw)
    sex = edf_unknown(first, sex_raw)
    dob = edf_unknown(raw->tryparse(Date, raw, dateformat"d-u-y"), dob_raw)
    dob === nothing && return
    name = edf_unknown(name_raw)
    return PatientID(code, sex, dob, name)
end

function Base.tryparse(::Type{RecordingID}, raw::AbstractString)
    s = split(raw, ' ', keepempty=false)
    length(s) == 5 || return
    first(s) == "Startdate" || return
    _, start_raw, admin_raw, tech_raw, equip_raw = s
    startdate = edf_unknown(raw->tryparse(Date, raw, dateformat"d-u-y"), start_raw)
    startdate === nothing && return
    admincode = edf_unknown(admin_raw)
    technician = edf_unknown(tech_raw)
    equipment = edf_unknown(equip_raw)
    return RecordingID(startdate, admincode, technician, equipment)
end

"""
    edf_unknown([f,] field::String)

Check whether the given text is simply "X", which per the EDF+ specification means that the
value is unknown, not applicable, or must be made anonymous, and if so return `missing`,
otherwise return the result of applying `f` to the given text. If unspecified, `f` defaults
to `identity`.
"""
edf_unknown(f, field::AbstractString) = field == "X" ? missing : f(field)
edf_unknown(field::AbstractString) = edf_unknown(identity, field)

"""
    set!(io::IO, [f=strip,] sigs::Vector{Signal}, name::Symbol, sz::Integer)

For each signal in `sigs`, set the field given by `name` to the result of `f` applied
to a chunk of size `sz` read from `io` as a `String`.
"""
function set!(io::IO, f, sigs::Vector{Signal}, name::Symbol, sz::Integer)
    n = length(sigs)
    T = fieldtype(Signal, name)
    @inbounds for i = 1:n
        x = f(String(Base.read(io, sz)))
        setfield!(sigs[i], name, convert(T, x))
    end
    return sigs
end

function set!(io::IO, sigs::Vector{Signal}, name::Symbol, sz::Integer)
    return set!(io, strip, sigs, name, sz)
end

parse_float(raw::AbstractString) = something(tryparse(Float32, raw), NaN32)

#####
##### Reading
#####

function read_header(io::IO)
    version = strip(String(Base.read(io, 8)))

    patient_id_raw = strip(String(Base.read(io, 80)))
    patient_id = something(tryparse(PatientID, patient_id_raw), String(patient_id_raw))

    recording_id_raw = strip(String(Base.read(io, 80)))
    recording_id = something(tryparse(RecordingID, recording_id_raw), String(recording_id_raw))

    start_raw = Base.read(io, 8)
    push!(start_raw, 0x20)  # Push a space separator
    append!(start_raw, Base.read(io, 8))  # Add the time
    # Parsing the date per the given format will validate EDF+ item 2
    start = DateTime(String(start_raw), dateformat"dd\.mm\.yy HH\.MM\.SS")
    if year(start) <= 84  # 1985 is used as a clipping date
        start += Year(2000)
    else
        start += Year(1900)
    end
    # FIXME: EDF "avoids" the Y2K problem by punting it to the year 2084, after which
    # we ignore the above entirely and use `recording_id.startdate`. We could add a
    # check here on `year(today())`, but that will be dead code for the next 60+ years.

    nb_header = parse(Int, String(Base.read(io, 8)))
    reserved = String(Base.read(io, 44))
    continuous = !startswith(reserved, "EDF+D")
    n_records = parse(Int, String(Base.read(io, 8)))
    duration = parse(Float64, String(Base.read(io, 8)))
    n_signals = parse(Int, String(Base.read(io, 4)))

    signals = [Signal() for _ = 1:n_signals]

    # TODO: EDF+ allows floating point data which does not fit within the Int16 limits.
    # See https://edfplus.info/specs/edffloat.html for details. MNE seems to implement
    # this(?)

    set!(io, signals, :label, 16)
    set!(io, signals, :transducer, 80)
    set!(io, signals, :physical_units, 8)
    set!(io, parse_float, signals, :physical_min, 8)
    set!(io, parse_float, signals, :physical_max, 8)
    set!(io, parse_float, signals, :digital_min, 8)
    set!(io, parse_float, signals, :digital_max, 8)
    set!(io, signals, :prefilter, 80)
    set!(io, x->parse(Int16, x), signals, :n_samples, 8)

    skip(io, 32 * n_signals)  # Reserved

    for signal in signals
        signal.samples = Vector{Int16}()  # Otherwise it's #undef
    end

    anno_idx = findfirst(signal->signal.label == "EDF Annotations", signals)
    if anno_idx !== nothing
        n_signals -= 1
    end

    @assert position(io) == nb_header

    h = Header(version, patient_id, recording_id, continuous, start, n_records,
               duration, n_signals, nb_header)
    return (h, signals, anno_idx)
end

function read_data!(io::IO, signals::Vector{Signal}, header::Header, ::Nothing)
    for i = 1:header.n_records, signal in signals
        append!(signal.samples, reinterpret(Int16, Base.read(io, 2 * signal.n_samples)))
    end
    @assert eof(io)
    return (signals, nothing)
end

function read_data!(io::IO, signals::Vector{Signal}, header::Header, anno_idx::Integer)
    annos = RecordAnnotation[]
    for i = 1:header.n_records
        anno = RecordAnnotation()
        anno.annotations = AnnotationsList[]
        for (j, signal) in enumerate(signals)
            n_bytes = 2 * signal.n_samples
            data = Base.read(io, n_bytes)
            if j == anno_idx
                record = IOBuffer(data)
                while !eof(record) && Base.peek(record) != 0x0
                    toplevel, offset, duration, events = read_tal(record)
                    if toplevel
                        anno.offset = offset
                        anno.event = events
                        anno.n_bytes = n_bytes
                    else
                        push!(anno.annotations, AnnotationsList(offset, duration, events))
                    end
                end
            else
                append!(signal.samples, reinterpret(Int16, data))
            end
        end
        push!(annos, anno)
    end
    deleteat!(signals, anno_idx)
    @assert eof(io)
    return (signals, annos)
end

function read_tal(io::IO)
    c = Base.read(io, UInt8)
    # Read the offset from the start time declared in the header
    @assert c === 0x2b || c === 0x2d  # + or -
    sign = c === 0x2b ? 1 : -1
    buffer = UInt8[]
    while true
        c = Base.read(io, UInt8)
        if c === 0x14 || c === 0x15
            break
        end
        push!(buffer, c)
    end
    offset = sign * parse(Float64, String(buffer))
    # Read the duration, if present
    if c === 0x15
        duration = parse(Float64, String(readuntil(io, 0x14)))
    else
        duration = nothing
    end
    c = Base.read(io, UInt8)
    record = c === 0x14  # Whether the annotation applies to the entire record
    # Read the annotation text, if present
    if c !== 0x0
        raw = readuntil(io, 0x0)
        pushfirst!(raw, c)
        events = convert(Vector{String}, split(String(raw), '\x14', keepempty=false))
    else
        events = String[]
    end
    return (record, offset, duration, events)
end

"""
    EDF.read(file::AbstractString)

Read the given file and return an `EDF.File` object containing the parsed data.
"""
function read(file::AbstractString)
    open(file, "r") do io
        header, data, anno_idx = read_header(io)
        data, annos = read_data!(io, data, header, anno_idx)
        File(header, data, annos)
    end
end


text:
module RayTracer

using Zygote, Flux, Images, Distributed, Statistics
import Base.show

using Compat
import Compat.isnothing

# Rendering Utilities
include("utils.jl")
include("light.jl")
include("materials.jl")
include("objects.jl")
include("camera.jl")
include("optimize.jl")

# Acceleration Structures
include("bvh.jl")

# Renderers
include("renderers/blinnphong.jl")
include("renderers/rasterizer.jl")
include("renderers/accelerated_raytracer.jl")

# Image Utilities
include("imutils.jl")

# Differentiable Rendering
include("gradients/zygote.jl")
include("gradients/numerical.jl")

end


text:
using LinearAlgebra, Random

"""
    generate_procrustes_data(m, n, d, σ=0.1, T=Float64)

Generate random (generalized) Procrustes data.

# Input
- `m :: Integer`: number of images.
- `n :: Integer`: number of image landmarks.
- `d :: Integer`: dimension of image (landmark points).
- `σ :: Number = 0.1`: noise level adding to each rotated image.
- `T :: DataType = Float64`: data type.

# Output
- `A :: Vector{Matrix{T}}`: `m` rotated images.
- `S :: Matrix{Matrix{T}}`: matrix of `S[i, j]` for OTSM problem.
- `Ā :: Matrix{T}`: the true/center image.
- `A_manopt :: Array{T, 3}`: 3D array of data suitable for Manopt software.
"""
function generate_procrustes_data(
    m :: Integer,
    n :: Integer,
    d :: Integer,
    σ :: Number   = 0.1,
    T :: DataType = Float64
    )
    # true center image
    Ā = randn(T, n, d)
    # randomly rotated images
    A = Vector{Matrix{T}}(undef, m)
    for k in 1:m
        # Haar measure (uniform) on O(n)
        # M with iid std normal entries
        # M = QR
        Q, R = qr!(randn(T, d, d))
        O = Matrix(Q)
        # In Julia, diagonal entries of R not necessarily positive
        for j in 1:d
            R[j, j] < 0 && (O[:, j] *= -1)
        end
        # now O is uniform on O(n), not on SO(n)
        # if det(O)=-1 (reflection), swap the first two columns
        det(O) < 0 && (O[:, [1, 2]] = O[:, [2, 1]])
        # rotate each point (row) in A
        A[k] = Ā * transpose(O) + σ * randn(T, n, d)
    end
    # set S[i, j]
    S = [A[i]'A[j] for i in 1:m, j in 1:m]
    for i in 1:m
        fill!(S[i, i], 0)
    end
    A_manopt = Array{T, 3}(undef, d, n, m)
    for k in 1:m
        A_manopt[:, :, k] = transpose(A[k])
    end
    A, S, Ā, A_manopt
end


text:
abstract type AbstractPotentialSimulationSetup{T, N} end

"""
    PotentialSimulationSetup{T, N, S} <: AbstractPotentialSimulationSetup{T, N}

Collection struct. It holds the grid, the potential, the point types, the charge density and the dielectric distribution.
"""
struct PotentialSimulationSetup{T, N, S} <: AbstractPotentialSimulationSetup{T, N}
    grid::Grid{T, N, S}
    potential::Array{T, N}
    pointtypes::Array{PointType, N}
    ρ::Array{T, N}
    ϵ::Array{T, N}
end

struct PotentialSimulationSetupRB{T, N1, N2, S} <: AbstractPotentialSimulationSetup{T, N1}
    grid::Grid{T, N1, S}
    potential::Array{T, N2}
    pointtypes::Array{PointType, N2}
    volume_weights::Array{T, N2}
    ρ::Array{T, N2}
    ϵ::Array{T, N1}
    geom_weights::NTuple{N1, AbstractGeometricalAxisWeights{T}}
    sor_const::Array{T, 1}
    bias_voltage::T
    maximum_applied_potential::T
    minimum_applied_potential::T
    depletion_handling_potential_limit::T
    bulk_is_ptype::Bool
    grid_boundary_factors::NTuple{3, NTuple{2, T}}
end

function sizeof(fssrb::PotentialSimulationSetup{T, N})::Int where {T, N}
    s::Int = sizeof(fssrb.grid)
    s += sizeof(fssrb.pointtypes)
    s += sizeof(fssrb.potential)
    s += sizeof(fssrb.ϵ)
    s += sizeof(fssrb.ρ)
    return s
end

function sizeof(fssrb::PotentialSimulationSetupRB{T, N1, N2})::Int where {T, N1, N2}
    s::Int = sizeof(fssrb.grid)
    s += sizeof(fssrb.pointtypes)
    s += sizeof(fssrb.potential)
    s += sizeof(fssrb.volume_weights)
    s += sizeof(fssrb.ϵ)
    s += sizeof(fssrb.ρ)
    for idim in 1:N1
        s += sizeof(fssrb.geom_weights[idim].weights)
    end
    s += sizeof(fssrb.sor_const)
    s += sizeof(fssrb.bias_voltage)
    s += sizeof(fssrb.maximum_applied_potential)
    s += sizeof(fssrb.minimum_applied_potential)
    s += sizeof(fssrb.depletion_handling_potential_limit)
    s += sizeof(fssrb.bulk_is_ptype)
    s += sizeof(fssrb.grid_boundary_factors)
    return s
end

function PotentialSimulationSetup(ssd::SolidStateDetector{T, S}, grid::Grid{T, N, S} = Grid(ssd), potential_array::Union{Missing, Array{T, N}} = missing)::PotentialSimulationSetup{T, N, S} where {T, N, S}
    fssrb::PotentialSimulationSetupRB{T, N, N + 1, S} = PotentialSimulationSetupRB(ssd, grid, potential_array)
    return PotentialSimulationSetup{T, N, S}( Grid(fssrb), ElectricPotentialArray(fssrb), PointTypeArray(fssrb), ChargeDensityArray(fssrb), DielektrikumDistributionArray(fssrb)  )
end

include("BoundaryConditions/BoundaryConditions.jl")

include("PotentialSimulationSetupRBCylindrical.jl")
include("PotentialSimulationSetupRBCartesian3D.jl")

include("plot_recipes.jl")


text:
### Day 7 - Problem 1
## Internet Protocol Version 7
## Authors: Thanasis Georgiou, Marina Grigoropoulou

# Looks for an ABBA identifier
function findABBA(str)
    len = length(str)
    if len < 4
        return false
    end

    for i = 2:len - 2
        if str[i] == str[i + 1] && str[i - 1] == str[i + 2] && str[i] != str[i - 1]
            return true
        end
    end
    return false
end

# Read input
f = open("input.txt")
input = readlines(f)
close(f)

# How many support TLS
count = 0

# For each IP
for ip in input
    # Split in string that are contained in brackets
    # and ones that are not
    bracketArr = matchall(r"\[\w*\]", ip)
    map!(str -> str[2:end-1], bracketArr)

    inBracket = join(bracketArr, '-')
    outBracket = replace(ip, r"\[\w*\]", "-")

    # Find ABBA identifiers in brackets
    abbaIn = findABBA(inBracket)
    abbaOut = findABBA(outBracket)

    # If it's valid increment counter
    if abbaIn == false && abbaOut == true
        count += 1
    end
end

println("Valid IPs: $count")


text:
export current_commit, tag!
export dict_list, dict_list_count

"""
    current_commit(path = projectdir()) -> commit
Return the current active commit id of the Git repository present
in `path`, which by default is the project path. If the repository
is dirty when this function is called the string will end
with `"_dirty"`.

Return `nothing` if `path` is not a Git repository.

See also [`tag!`](@ref).

## Examples
```julia
julia> current_commit()
"96df587e45b29e7a46348a3d780db1f85f41de04"

julia> current_commit(path_to_dirty_repo)
"3bf684c6a115e3dce484b7f200b66d3ced8b0832_dirty"
```
"""
function current_commit(path = projectdir())
    # Here we test if the path is a git repository.
    try
        repo = LibGit2.GitRepo(path)
    catch er
        @warn "The current project directory is not a Git repository, "*
        "returning `nothing` instead of the commit id."
        return nothing
    end
    # then we return the current commit
    repo = LibGit2.GitRepo(path)
    c = string(LibGit2.head_oid(repo))
    if LibGit2.isdirty(repo)
        @warn "The Git repository is dirty! Adding appropriate comment to "*
        "commit id..."
        return c*"_dirty"
    end
    return c
end

"""
    tag!(d::Dict, path = projectdir()) -> d
Tag `d` by adding an extra field `commit` which will have as value
the [`current_commit`](@ref) of the repository at `path` (by default
the project's path). Does not operate if a key `commit` already exists.

Notice that if `String` is not a subtype of the value type of `d` then
a new dictionary is created and returned. Otherwise the operation
is inplace (and the dictionary is returned again).

## Examples
```julia
julia> d = Dict(:x => 3, :y => 4)
Dict{Symbol,Int64} with 2 entries:
  :y => 4
  :x => 3

julia> tag!(d)
Dict{Symbol,Any} with 3 entries:
  :y      => 4
  :commit => "96df587e45b29e7a46348a3d780db1f85f41de04"
  :x      => 3
```
"""
function tag!(d::Dict{K, T}, path = projectdir()) where {K, T}
    c = current_commit(path)
    c === nothing && return d
    if haskey(d, K("commit"))
        @warn "The dictionary already has a key named `commit`. We won't "
        "overwrite it with the commit id."
        return d
    end
    if String <: T
        d[K("commit")] = c
    else
        d = Dict{K, promote_type(T, String)}(d)
        d[K("commit")] = c
    end
    return d
end

"""
    dict_list(c)
Expand the dictionary `c` into a vector of dictionaries.
Each entry has a unique combination from the product of the `Vector`
values of the dictionary while the non-`Vector` values are kept constant
for all possibilities. The keys of the entries are the same.

Whether the values of `c` are iterable or not is of no concern;
the function considers as "iterable" only subtypes of `Vector`.

Use the function [`dict_list_count`](@ref) to get the number of
dictionaries that `dict_list` will produce.

## Examples
```julia
julia> c = Dict(:a => [1, 2], :b => 4);

julia> dict_list(c)
3-element Array{Dict{Symbol,Int64},1}:
 Dict(:a=>1,:b=>4)
 Dict(:a=>2,:b=>4)

julia> c[:model] = "linear"; c[:run] = ["bi", "tri"];

julia> dict_list(c)
4-element Array{Dict{Symbol,Any},1}:
 Dict(:a=>1,:b=>4,:run=>"bi",:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:model=>"linear")

julia> c[:e] = [[1, 2], [3, 5]];

julia> dict_list(c)
8-element Array{Dict{Symbol,Any},1}:
 Dict(:a=>1,:b=>4,:run=>"bi",:e=>[1, 2],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:e=>[1, 2],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:e=>[1, 2],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:e=>[1, 2],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"bi",:e=>[3, 5],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:e=>[3, 5],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:e=>[3, 5],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:e=>[3, 5],:model=>"linear")
```
"""
function dict_list(c)
    iterable_fields = filter(k -> typeof(c[k]) <: Vector, keys(c))
    non_iterables = setdiff(keys(c), iterable_fields)

    iterable_dict = Dict(iterable_fields .=> getindex.(Ref(c), iterable_fields))
    non_iterable_dict = Dict(non_iterables .=> getindex.(Ref(c), non_iterables))

    vec(
        map(Iterators.product(values(iterable_dict)...)) do vals
            dd = Dict(keys(iterable_dict) .=> vals)
            merge(non_iterable_dict, dd)
        end
    )
end

"""
    dict_list_count(c) -> N
Return the number of dictionaries that will be created by
calling `dict_list(c)`.
"""
function dict_list_count(c)
    iterable_fields = filter(k -> typeof(c[k]) <: Vector, keys(c))
    prod(length(c[i]) for i in iterable_fields)
end

# function ntuple_list(c)
#     iterable_fields = filter(k -> typeof(c[k]) <: Vector, collect(keys(c)))
#     non_iterables = setdiff(keys(c), iterable_fields)
#
#     iterable_vals = Tuple(collect((typeof(c[i]) for i in iterable_fields)))
#     iterable_tuple = NamedTuple{tuple(iterable_fields), iterable_vals}
#
#     Dict(iterable_fields .=> getindex.(Ref(c), iterable_fields))
#     non_iterable_tuple = Dict(
#         non_iterables .=> getindex.(Ref(c), non_iterables))
#
#     vec(
#         map(Iterators.product(values(iterable_dict)...)) do vals
#             dd = Dict(keys(iterable_dict) .=> vals)
#             merge(non_iterable_dict, dd)
#         end
#     )
# end


text:
#############################################################################
# set_variables.jl
# Defines SetVariable, which is a subtype of AbstractExpr
#############################################################################

export SetVariable
export fix!, free!, elements
export get_cv, get_sv, get_v

mutable struct SetVariable <: CombiSet
  head::Symbol
  id_hash::UInt64
  elements::ValOrNothing
  baseset::AbstractArray
  cardinality::Int

  function SetVariable(baseset::AbstractArray)
    this = new(:setvariable, 0, nothing, baseset, length(baseset))
    this.id_hash = object_id(this)
    return this
  end
  SetVariable(m::Int) = SetVariable(collect(1:m))
end

id_to_variables = Dict{UInt64, SetVariable}()

# fix set variables to hold them at their current value, and free them afterwards
function fix!(x::SetVariable)
  x.elements == nothing && error("This set variable has not been assigned elements yet; cannot fix value to nothing!")
  x
end

function fix!(x::SetVariable, v)
  # TODO: check value inclusion
  x.elements = v
  fix!(x)
end

function free!(x::SetVariable)
  # TODO this won't work if :fixed appears other than at the end of x.sets
  x.sets[end] == :fixed && pop!(x.sets)
  x
end

### Output the variables of an expression
function get_cv(x::AbstractExpr)
  if isa(x, Variable) || (typeof(x) == Constant) || (typeof(x) == SetVariable)
    if isa(x, Variable)
      cv = Variable[]
      push!(cv, x)
      return cv
    else
      return Variable[]
    end
  else
    cv = Variable[]
    for i = 1:length(x.children)
      cvi = get_cv(x.children[i])
      cv = vcat(cv, cvi)
    end
    cvSet = Set(cv)
    cv = collect(cvSet)
    return cv
  end
end

function get_cv(x::Values)
  return Variable[]
end

### Output the set variables of an expression
function get_sv(x::AbstractExpr)
  if isa(x, Variable) || (typeof(x) == Constant) || (typeof(x) == SetVariable)
    if isa(x, SetVariable)
      sv = SetVariable[]
      push!(sv, x)
      return sv
    else
      return SetVariable[]
    end
  else
    sv = SetVariable[]
    for i = 1:length(x.children)
      svi = get_sv(x.children[i])
      sv = vcat(sv, svi)
    end
    svSet = Set(sv)
    sv = collect(svSet)
    return sv
  end
end

function get_sv(f::SubmodFunc)
  return f.setvariables
end

function get_sv(x::Values)
  return []
end

### Output all the variables of an expression
function get_v(x::AbstractExpr)
  if isa(x, Variable) || (typeof(x) == Constant) || (typeof(x) == SetVariable)
    if isa(x, Variable) || typeof(x) == SetVariable
      v = AbstractExpr[]
      push!(v, x)
      return v
    else
      return []
    end
  else
    v = AbstractExpr[]
    for i = 1:length(x.children)
      vi = get_v(x.children[i])
      v = vcat(v, vi)
    end
    vSet = Set(v)
    v = collect(vSet)
    return v
  end
end

function get_v(f::SubmodFunc)
  return f.setvariables
end

function get_v(x::Values)
  return []
end


text:
using Test
using LinearAlgebra
using SparseArrays
using WriteVTK
using WaveProp
using WaveProp.Mesh
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Interpolation
using WaveProp.FEM
using WaveProp.IO

## Assembly test
@testset "Mass matrix 1-2-3D Lagrange P0 P1" begin
    ## 2D
    Geometry.clear!()
    Ω, mesh2d = WaveProp.IO.gmsh_rectangle(;h=0.1);
    for ord in [1,2,]
        ## matrix assembly
        for (pu, pv) in [(1, 1), (1, 1), (0, 1), (1, 0),]
            u = LagrangeBasis{ReferenceTriangle,pu}()
            v = LagrangeBasis{ReferenceTriangle,pv}()
            @time A = assembly(mesh2d, Ω, u, v; order=ord)
            m, n = size(A)
            surface = ones(1, m) * A * ones(n)
            @test abs(surface[1] - 1) < 1.e-12
            @show m,n
            u = LagrangeBasis{ReferenceLine,pu}()
            v = LagrangeBasis{ReferenceLine,pv}()
            @time A = assembly(mesh2d, boundary(Ω), u, v; order=ord)
            m, n = size(A)
            perimeter = ones(1, m) * A * ones(n)
            @test abs(perimeter[1] - 4) < 1.e-12
            @show m,n
        end
    end
    ## 3D
    Geometry.clear!()
    Ω, mesh3d = WaveProp.IO.gmsh_box(;h=0.5);
    for ord in [1,2,]
        ## matrix assembly
        for (pu, pv) in [(0, 0), (1, 1), (0, 1), (1, 0),]
            u = LagrangeBasis{ReferenceTetrahedron,pu}()
            v = LagrangeBasis{ReferenceTetrahedron,pv}()
            @time A = assembly(mesh3d, Ω, u, v; order=ord)
            m, n = size(A)
            volume = ones(1, m) * A * ones(n)
            @test abs(volume[1] - 1) < 1.e-12
            @show m,n
            u = LagrangeBasis{ReferenceTriangle,pu}()
            v = LagrangeBasis{ReferenceTriangle,pv}()
            @time A = assembly(mesh3d, boundary(Ω), u, v; order=ord)
            m, n = size(A)
            surface = ones(1, m) * A * ones(n)
            @test abs(surface[1] - 6) < 1.e-12
            @show m,n
            #u = LagrangeBasis{ReferenceLine,pu}()
            #v = LagrangeBasis{ReferenceLine,pv}()
            #@show skeleton(boundary(Ω)) |> tags |> length
            #@time A = assembly(mesh3d, skeleton(boundary(Ω)), u, v; order=ord)
            #m, n = size(A)
            #edgesLength = ones(1, m) * A * ones(n)
            #@test abs(edgesLength[1] - 12) < 1.e-12
            #@show m,n
        end
    end
end

## Helmholtz
k = 1
a_helmholtz(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i]) - k^2 * u(x̂)[j] * v(x̂)[i]
@testset "Helmholtz matrix 1-2-3D Lagrange P1" begin
    ## 2D
    Geometry.clear!()
    Ω, mesh2d = WaveProp.IO.gmsh_rectangle(;h=0.1);
    for ord in [1,2,]
        ## matrix assembly
        for (pu, pv) in [(1, 1),]
            u = LagrangeBasis{ReferenceTriangle,pu}()
            v = LagrangeBasis{ReferenceTriangle,pv}()
            @time A = assembly(mesh2d, Ω, u, v; order=ord, f=a_helmholtz)
            m, n = size(A)
            @show m,n
        end
    end
    ## 3D
    Ω, mesh3d = WaveProp.IO.gmsh_box(;h=0.5);
    for ord in [1,2,]
        ## matrix assembly
        for (pu, pv) in [(1, 1),]
            u = LagrangeBasis{ReferenceTetrahedron,pu}()
            v = LagrangeBasis{ReferenceTetrahedron,pv}()
            @time A = assembly(mesh3d, Ω, u, v; order=ord, f=a_helmholtz)
            m, n = size(A)
            @show m,n
        end
    end
end

## Poisson test 2D
p, q = 2, 3
sol(x,y) = sin(p*π*x) * sin(q*π*y)
f(x,y) = (p^2 + q^2) * π^2 * sol(x,y)
am(u, v) = (i,j,x̂,el,x) -> u(x̂)[j] * v(x̂)[i]
a_poisson(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i])
l_poisson(v) = (i,x̂,el,x) -> f(x...) * v(x̂)[i]
Geometry.clear!()
Ω, mesh2d = WaveProp.IO.gmsh_rectangle(;h=0.05);
u = LagrangeBasis{ReferenceTriangle,1}()
v = LagrangeBasis{ReferenceTriangle,1}()
M = assembly(mesh2d, Ω, u, v; order=2, f=am);
A = assembly(mesh2d, Ω, u, v; order=2, f=a_poisson);
b = assembly(mesh2d, Ω, v; order=2, f=l_poisson);
# Dirichlet BC
vΩdofnb = local_dof_numbering(mesh2d, Ω, v);
vΓdofnb = local_dof_numbering(mesh2d, boundary(Ω), trace(v));
glob2loc = SparseVector(length(nodes(mesh2d)),dofs(vΩdofnb),collect(1:length(dofs(vΩdofnb))));
ϵ = 1e-16
for dof in dofs(vΓdofnb)
    i = glob2loc[dof]
    A[i,i] = 1 / ϵ
    b[i] = 0
end
# Solution
u = factorize(A) \ b;
# Export
vtkfile = vtk_mesh_file(mesh2d, Ω, joinpath(@__DIR__,"poisson"))
vtkfile["u", VTKPointData()] = u;
vtkfile |> vtk_save
# Test
uref = [sol(x...) for x in nodes(mesh2d)]
e = u .- uref
err = sqrt(transpose(e) * (M * e)) / sqrt(transpose(uref) * (M * uref))
@show err
# Should be quadratic in L2 norm (and it is)
@test err < 2e-3 # for h = 0.05

## Poisson test 3D
p, q, r = 2, 3, 4
sol(x,y,z) = sin(p*π*x) * sin(q*π*y) * sin(r*π*z)
f(x,y,z) = (p^2 + q^2 + r^2) * π^2 * sol(x,y,z)
am(u, v) = (i,j,x̂,el,x) -> u(x̂)[j] * v(x̂)[i]
a_poisson(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i])
l_poisson(v) = (i,x̂,el,x) -> f(x...) * v(x̂)[i]

Ω, mesh3d = WaveProp.IO.gmsh_box(;h=0.05);
u = LagrangeBasis{ReferenceTetrahedron,1}()
v = LagrangeBasis{ReferenceTetrahedron,1}()
A = assembly(mesh3d, Ω, u, v; order=2, f=a_poisson);
M = assembly(mesh3d, Ω, u, v; order=2, f=am);
b = assembly(mesh3d, Ω, v; order=2, f=l_poisson);
# Dirichlet BC
vΩdofnb = local_dof_numbering(mesh3d, Ω, v);
vΓdofnb = local_dof_numbering(mesh3d, boundary(Ω), trace(v));
glob2loc = SparseVector(length(nodes(mesh3d)),dofs(vΩdofnb),collect(1:length(dofs(vΩdofnb))));
ϵ = 1e-16
for dof in dofs(vΓdofnb)
    i = glob2loc[dof]
    A[i,i] = 1 / ϵ
    b[i] = 0
end
# Solution
u = factorize(A) \ b;
# Export
vtkfile = vtk_mesh_file(mesh3d, Ω, joinpath(@__DIR__,"poisson"))
vtkfile["u", VTKPointData()] = u;
vtkfile |> vtk_save
# Test
uref = [sol(x...) for x in nodes(mesh3d)]
e = u .- uref
err = sqrt(transpose(e) * (M * e)) / sqrt(transpose(uref) * (M * uref))
@show err
# Should be quadratic in L2 norm (and it is)
@test err < 4e-2 # for h = 0.05

# ## Modes for 2D Neumann Laplacian
# ak(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i])
# am(u, v) = (i,j,x̂,el,x) -> u(x̂)[j] * v(x̂)[i]
# Geometry.clear!()
# Ω, mesh2d = WaveProp.IO.gmsh_disk(;h=0.01)
# u = LagrangeBasis{ReferenceTriangle,1}()
# v = LagrangeBasis{ReferenceTriangle,1}()
# K = assembly(mesh2d, Ω, u, v; order=2, f=ak);
# M = assembly(mesh2d, Ω, u, v; order=2, f=am);
# F = eigen(Matrix(K), Matrix(M))
# @show F.values |> length
# for km in 1:min(100,length(F.values))
#     vtkfile = vtk_mesh_file(mesh3d, Ω, joinpath(@__DIR__,"modes_$(km)"))
#     vtkfile["u", VTKPointData()] = F.vectors[:,km];
#     vtkfile |> vtk_save
# end

# ## Modes for 3D Neumann Laplacian
# ak(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i])
# am(u, v) = (i,j,x̂,el,x) -> u(x̂)[j] * v(x̂)[i]
# Ω, mesh3d = WaveProp.IO.gmsh_sphere(;h=0.05);
# u = LagrangeBasis{ReferenceTetrahedron,1}()
# v = LagrangeBasis{ReferenceTetrahedron,1}()
# K = assembly(mesh3d, Ω, u, v; order=2, f=ak);
# M = assembly(mesh3d, Ω, u, v; order=2, f=am);
# F = eigen(Matrix(K), Matrix(M))
# @show F.values |> length
# for km in 1:min(100,length(F.values))
#     vtkfile = vtk_mesh_file(mesh3d, Ω, joinpath(@__DIR__,"modes_$(km)"))
#     vtkfile["u", VTKPointData()] = F.vectors[:,km];
#     vtkfile |> vtk_save
# end

# ## Modes for 2D Neumann Laplacian on 3D sphere
# # WARNING: I am not sure that it is what I meant it to be though
# ak(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))[1:2,1:2]) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))[1:2,1:2]) *  grad(v)(x̂)[i])
# am(u, v) = (i,j,x̂,el,x) -> u(x̂)[j] * v(x̂)[i]
# Ω, mesh3d = WaveProp.IO.gmsh_sphere(;h=0.05);
# u = LagrangeBasis{ReferenceTriangle,1}()
# v = LagrangeBasis{ReferenceTriangle,1}()
# Γ = boundary(Ω)
# K = assembly(mesh3d, Γ, u, v; order=2, f=ak);
# M = assembly(mesh3d, Γ, u, v; order=2, f=am);
# F = eigen(Matrix(K), Matrix(M))
# @show F.values |> length
# for km in 1:min(100,length(F.values))
#     vtkfile = vtk_mesh_file(mesh3d, Γ, joinpath(@__DIR__,"modes_sphere_$(km)"))
#     vΩdofnb = local_dof_numbering(mesh3d, Ω, LagrangeBasis{ReferenceTetrahedron,1}());
#     vΓdofnb = local_dof_numbering(mesh3d, Γ, v);
#     glob2loc = SparseVector(length(nodes(mesh3d)),dofs(vΩdofnb),collect(1:length(dofs(vΩdofnb))));
#     Fkm = zeros(length(nodes(mesh3d)))
#     for dof in dofs(vΓdofnb)
#         i = glob2loc[dof]
#         Fkm[i] = F.vectors[:,km][i]
#     end
#     vtkfile["u", VTKPointData()] = Fkm;
#     vtkfile |> vtk_save
# end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
struct Point3D <: FieldVector{3, Float64}
    x::Float64
    y::Float64
    z::Float64
end

import Base: +, -, ==, *, ^, <, >, /, !=
import Base: isnan, isapprox, isequal, isless, broadcastable
import Luxor: between, distance, midpoint

Base.broadcastable(x::Point3D) = Ref(x)

function +(p1::Point3D, p2::Point3D)
    Point3D((p2.x + p1.x), (p2.y + p1.y), (p2.z + p1.z))
end
+(p1::Point3D, k::Number)              = Point3D(p1.x + k, p1.y + k, p1.z + k)
+(k::Number, p2::Point3D)              = +(p2::Point3D, k::Number)

function -(p1::Point3D, p2::Point3D)
    Point3D((p1.x - p2.x), (p1.y - p2.y), (p1.z - p2.z))
end
-(p::Point3D)                          = Point3D(-p.x, -p.y, -p.z)
-(k::Number, p1::Point3D)              = Point3D(p1.x - k, p1.y - k, p1.z - k)
-(p1::Point3D, k::Number)              = Point3D(p1.x - k, p1.y - k, p1.z - k)

*(k::Number, p2::Point3D)              = Point3D(k * p2.x, k * p2.y, k * p2.z)
*(p2::Point3D, k::Number)              = Point3D(k * p2.x, k * p2.y, k * p2.z)
*(p1::Point3D, p2::Point3D)            = Point3D(p1.x * p2.x, p1.y * p2.y, p1.z * p2.z)

"""
    distance(p1::Point3D, p2::Point3D)

Return the distance between two points.

"""
function distance(p1::Point3D, p2::Point3D)
    sqrt((p2.x - p1.x)^2 + (p2.y - p1.y)^2 + (p2.z - p1.z)^2)
end

function isnan(pt::Point3D)
    isnan(pt.x) || isnan(pt.y) || isnan(pt.z)
end

"""
    between(p1::Point3D, p2::Point3D, x=0.5)
    between((p1::Point3D, p2::Point3D), x=0.5)

Find a point on a line between two 3D points.
If `x` is 0.5, the returned point should be halfway between them.

"""
function between(p1::Point3D, p2::Point3D, x=0.5)
    return p1 + (x * (p2 - p1))
end

function between(couple::NTuple{2, Point3D}, x=0.5)
    p1, p2 = couple
    return p1 + (x * (p2 - p1))
end

/(p2::Point3D, k::Number)              = Point3D(p2.x/k, p2.y/k, p2.z/k)
^(p::Point3D, e::Integer)              = Point3D(p.x^e,  p.y^e, p.z^e)
^(p::Point3D, e::Float64)              = Point3D(p.x^e,  p.y^e, p.z^e)

# conversion

function Base.convert(type::Type{Point3D}, pt::Point)
    return Point3D(pt.x, pt.y, 0.0)
end

# some refinements
# modifying points with tuples
+(p1::Point3D, shift::NTuple{3, Real}) = Point3D(p1.x + shift[1], p1.y + shift[2], p1.z + shift[3])
-(p1::Point3D, shift::NTuple{3, Real}) = Point3D(p1.x - shift[1], p1.y - shift[2], p1.z - shift[3])
*(p1::Point3D, shift::NTuple{3, Real}) = Point3D(p1.x * shift[1], p1.y * shift[2], p1.z * shift[3])
/(p1::Point3D, shift::NTuple{3, Real}) = Point3D(p1.x / shift[1], p1.y / shift[2], p1.z / shift[3])

# comparisons

isequal(p1::Point3D, p2::Point3D)         = isapprox(p1.x, p2.x, atol=0.00000001) && isapprox(p1.y, p2.y, atol=0.00000001) && isapprox(p1.z, p2.z, atol=0.00000001)
isapprox(p1::Point3D, p2::Point3D)        = isapprox(p1.x, p2.x, atol=0.00000001) && isapprox(p1.y, p2.y, atol=0.00000001) && isapprox(p1.z, p2.z, atol=0.00000001)
isless(p1::Point3D, p2::Point3D)          = (p1.x < p2.x) && (p1.y < p2.y) && (p1.z < p2.z)
!=(p1::Point3D, p2::Point3D)              = !isequal(p1, p2)
<(p1::Point3D, p2::Point3D)               = isless(p1,p2)
>(p1::Point3D, p2::Point3D)               = p2 < p1
==(p1::Point3D, p2::Point3D)              = isequal(p1, p2)

"""
    dotproduct3D(a::Point3D, b::Point3D)

Finds the dot product of a and b
"""
function dotproduct3D(a::Point3D, b::Point3D)
    return a.x * b.x + a.y * b.y + a.z * b.z
end

"""
    magnitude(a::Point3D)

Calculates magnitude of a.
"""
function magnitude(a::Point3D)
    return sqrt(a.x * a.x + a.y * a.y + a.z * a.z)
end

"""
    anglebetweenvectors(v1::Point3D, v2::Point3D)

Calclate the angle between two vectors.
"""
function anglebetweenvectors(v1::Point3D, v2::Point3D)
    intermediate = dotproduct3D(v1, v2)/(magnitude(v1) * magnitude(v2))
    # avoid domain errors
    if isnan(intermediate)
        return 0.0
    else
        return acos(min(max(-1, intermediate), 1))
    end
end

"""
    midpoint(pt1::Point3D, pt2::Point3D)

Find the midpoint between two points. See also `between()`.
"""
function midpoint(pt1::Point3D, pt2::Point3D)
    return (pt1 + pt2)/2
end

"""
    sphericaltocartesian(ρ, θ, ϕ)

Return `Point3D(x, y, z)` corresponding to `(ρ, θ, ϕ)`:

- ρ is the distance from the origin (ie radius)

- θ is the azimuthal angle (the longitude) 0 is +x, π is -x, 2π is +x

- ϕ is the polar angle (the latitude) 0 is North Pole, π is South Pole

There are two major conventions for spherical coordinate notation.

In physics books:

(ρ, θ, φ) gives the radial distance, polar angle (latitude), and azimuthal angle (longitude)

In mathematics books:

(ρ, θ , φ ) gives the radial distance, azimuthal angle (longitude), and polar angle (latitude)

So we're using the mathematics one here.
"""
function sphericaltocartesian(ρ, θ, ϕ)
    x = ρ * sin(ϕ) * cos(θ)
    y = ρ * sin(ϕ) * sin(θ)
    z = ρ * cos(ϕ)
    return Point3D(x, y, z)
end

"""
    cartesiantospherical(x, y, z)

Return `(ρ, θ, ϕ)` (radius, longitude, latitude) of the Point3D(x, y, z).
"""
function cartesiantospherical(x, y, z)
    ϕ = atan(y, x)
    ρ = sqrt(x^2 + y^2 + z^2)
    θ = acos(z/ρ)
    return (ρ, θ, ϕ)
end

# rotations

# old versions just redirect to Rotations now
"""
    rotateX(pt3D::Point3D, rad)

Return a new point resulting from rotating the point around the x axis by an angle in radians.

Rotations are anticlockwise when looking along axis from 0 to +axis.
"""
function rotateX(pt3D::Point3D, rad)
    return RotX(rad) * pt3D
end

"""
    rotateY(pt3D::Point3D, rad)

Return a new point resulting from rotating the point around the y axis by an angle in radians.

Rotations are anticlockwise when looking along axis from 0 to +axis.
"""
function rotateY(pt3D::Point3D, rad)
    return RotY(rad) * pt3D
end

"""
    rotateZ(pt3D::Point3D, rad)

Return a new point resulting from rotating the point around the z axis by an angle in radians.
"""
function rotateZ(pt3D::Point3D, rad)
    return RotZ(rad) * pt3D
end

# rotate around axes

## copies (new points)
"""
    rotateby(pt::Point3D, angleX, angleY, angleZ)
    rotateby(ptlist::Array{Point3D, 1}, angleX, angleY, angleZ)
    rotateby(point::Point3D, r::Rotation)
    rotateby(ptlist::Array{Point3D, 1}, r::Rotation)

Return a new point/list of points resulting from rotating
around the x, y, and z axes by angleX, angleY, angleZ.

The Z rotation is first, then the Y, then the X.

A 3×3 rotation matrix parameterized by the
"Tait-Bryant" XYZ Euler angle convention,
consisting of first a rotation about the Z
axis by theta3, followed by a rotation about
the Y axis by theta2, and finally a rotation
about the X axis by theta1.
"""
function rotateby(point::Point3D, angleX::Float64, angleY::Float64, angleZ::Float64)
    return RotXYZ(angleX, angleY, angleZ) * point
end

function rotateby(ptlist::Array{Point3D, 1}, angleX::Float64, angleY::Float64, angleZ::Float64)
    return rotateby.(ptlist, angleX, angleY, angleZ)
end

function rotateby(point::Point3D, r::Rotation)
    return r * point
end

function rotateby(ptlist::Array{Point3D, 1}, r::Rotation)
    return rotateby.(ptlist, r)
end

## update an array of points
"""
    rotateby!(ptlist::Array{Point3D, 1}, angleX, angleY, angleZ)

Modify a list of points by rotating each one around the
x, y, and z axes by angleX, angleY, angleZ.
"""
function rotateby!(ptlist::Array{Point3D, 1}, angleX::Float64, angleY::Float64, angleZ::Float64)
    for i in eachindex(ptlist)
        ptlist[i] = rotateby(ptlist[i], angleX, angleY, angleZ)
    end
    return ptlist
end

# rotate a point around another point

## copy (new point)

"""
    rotateby(point::Point3D, about::Point3D, angleX, angleY, angleZ)
    rotateby(point::Point3D, about::Point3D, r::Rotation)
    rotateby(ptlist::Array{Point3D, 1}, about::Point3D, r::Rotation)
"""
function rotateby(point::Point3D, about::Point3D, angleX::Float64, angleY::Float64, angleZ::Float64)
    return RotXYZ(angleX, angleY, angleZ) * (point - about) + about
end
function rotateby(point::Point3D, about::Point3D, r::Rotation)
    return r * (point - about) + about
end
function rotateby(ptlist::Array{Point3D, 1}, about::Point3D, r::Rotation)
    return rotateby.(ptlist, about, r)
end

## update an array of points
"""
    rotateby!(ptlist::Array{Point3D, 1}, existingpt::Point3D, angleX, angleY, angleZ)
    rotateby!(ptlist::Array{Point3D, 1}, existingpt::Point3D, r::Rotation)
    rotateby!(ptlist::Array{Point3D, 1}, r::Rotation=RotXYZ{Float64})

Rotate each point in the list by rotation (or angleX, angleY, angleZ) around another point (or origin).
"""
function rotateby!(ptlist::Array{Point3D, 1}, existingpt::Point3D, angleX::Float64, angleY::Float64, angleZ::Float64)
    for i in eachindex(ptlist)
        ptlist[i] = rotateby(ptlist[i], existingpt, angleX, angleY, angleZ)
    end
    return ptlist
end
function rotateby!(ptlist::Array{Point3D, 1}, existingpt::Point3D, r::Rotation)
    for i in eachindex(ptlist)
        ptlist[i] = rotateby(ptlist[i], existingpt, r)
    end
    return ptlist
end

rotateby!(ptlist::Array{Point3D, 1}, r::Rotation=RotXYZ{Float64}) = rotateby!(ptlist, Point3D(0.0, 0.0, 0.0), r)

"""
    moveby(pt::Point3D, d::Point3D)

Return a new point that's the result of moving a point `pt` by a vector `d`.
"""
function moveby(pt::Point3D, d::Point3D)
    return pt + d
end

"""
    moveby!(ptlist::Point3D, x, y, z)
    moveby!(ptlist::Point3D, pt::Point3D)

Move all points in the list by a vector.
"""
function moveby!(ptlist::Array{Point3D, 1}, pt::Point3D)
    return ptlist .+= pt
end

moveby!(ptlist::Array{Point3D, 1}, x, y, z) = moveby!(ptlist, Point3D(x, y, z))

"""
    scaleby!(ptlist::Array{Point3D, 1}, x, y, z)

Scales a list of points by multiplying by `x` in X, `y` in Y, `z` in Z.
"""
function scaleby!(ptlist::Array{Point3D, 1}, x, y, z)
   for n in 1:length(ptlist)
       v = ptlist[n]
       ptlist[n] = Point3D(v.x * x, v.y * y, v.z * z)
   end
   return ptlist
end

"""
    surfacenormal(ptlist::Array{Point3D, 1})

Finds one of these.
"""
function surfacenormal(ptlist::Array{Point3D, 1})
   normal = Point3D(0, 0, 0)
   for i in 1:length(ptlist)
      vertexCurrent = ptlist[i]
      vertexNext    = ptlist[mod1(i + 1, end)]
      x = normal.x + ( (vertexCurrent.y - vertexNext.y) * (vertexCurrent.z + vertexNext.z))
      y = normal.y + ( (vertexCurrent.z - vertexNext.z) * (vertexCurrent.x + vertexNext.x))
      z = normal.z + ( (vertexCurrent.x - vertexNext.x) * (vertexCurrent.y + vertexNext.y))
      normal = Point3D(x, y, z)
   end
   l = magnitude(normal)
   # return normalize(normal)
   return Point3D(normal.x/l, normal.y/l, normal.z/l) / l
end

function Base.unique(pts::Array{Point3D, 1})
    apts = Point3D[]
    for pt in pts
        if pt ∉ apts
            push!(apts, pt)
        end
    end
    return apts
end


text:
using Utils

using Chain

struct Point
    x::Int
    y::Int
end
Point() = Point(0, 0)

struct Area
    ll::Point
    ur::Point
end

struct Shot
    loc::Point
    vel::Point
end
Shot(p) = Shot(Point(), p)

const file = find_input(@__FILE__)
function load(file)
    content = read(file, String)
    nums = map(r -> parse(Int, content[r]), findall(r"-?\d+", content))
    Area(Point(nums[1], nums[3]), Point(nums[2], nums[4]))
end
const input = load(file)

function part1(input::Area)
    n = -input.ll.y - 1
    n*(n + 1)÷2
end

assertequal(part1(input), 3160)

@assert input.ll.x > 0
function Base.iterate(muzzle_velocity::Point, state=Shot(muzzle_velocity))
    # Won't hit anything ever again.
    ((state.vel.y < 0 && state.loc.y < input.ll.y) ||
    (state.vel.x > 0 && state.loc.x > input.ur.x)) && return nothing

    pos = Point(state.loc.x + state.vel.x, state.loc.y + state.vel.y)
    vel = Point(max(0, state.vel.x - 1), state.vel.y - 1)
    pos, Shot(pos, vel)
end

Base.occursin(p::Point, a::Area) = a.ll.x <= p.x <= a.ur.x && a.ll.y <= p.y <= a.ur.y

part2() = count(Iterators.product(1:input.ur.x, -99:99)) do tup
    any(occursin(input), Point(tup...))
end

assertequal(part2(), 1928)

max_stars_repo_name: CmdQ/AoC

max_issues_repo_name: CmdQ/AoC

max_forks_repo_name: CmdQ/AoC