text:
# Implements two view types for "converting" between arrays-of-colors
# and arrays-of-numbers (with the "first dimension" corresponding to
# color channels)
#  - ChannelView: view a color array as if it were an array of numbers
#  - ColorView: view an array of numbers as if it were an array of colors
# Examples:
#    img is a m-by-n Array{RGB{Float32}}
#    ChannelView(img) is a 3-by-m-by-n AbstractArray{Float32}
#
#    buffer is a 3-by-m-by-n Array{U8}
#    ColorView{RGB}(buffer) is an m-by-n AbstractArray{RGB{U8}}

# "First dimension" applies to colors like RGB; by default, Gray
# images don't use a whole dimension (of size 1) just to encode
# colors. But it's easy to change that behavior with the flip of a
# switch:
const squeeze1 = true # when true, don't use a dimension for the color channel of grayscale

typealias Color1{T} Colorant{T,1}
typealias Color2{T} Colorant{T,2}
typealias Color3{T} Colorant{T,3}
typealias Color4{T} Colorant{T,4}
typealias AColor{N,C,T} AlphaColor{C,T,N}
typealias ColorA{N,C,T} ColorAlpha{C,T,N}

## ChannelView

immutable ChannelView{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    parent::A

    function ChannelView{C<:Colorant}(parent::AbstractArray{C})
        n = length(channelview_indices(parent))
        n == N || throw(DimensionMismatch("for an $N-dimensional ChannelView with color type $C, input dimensionality should be $n instead of $(ndims(parent))"))
        new(parent)
    end
end

"""
    ChannelView(A)

creates a "view" of the Colorant array `A`, splitting out (if
necessary) the separate color channels of `eltype(A)` into a new first
dimension. For example, if `A` is a m-by-n RGB{U8} array,
`ChannelView(A)` will return a 3-by-m-by-n U8 array. Color spaces with
a single element (i.e., grayscale) do not add a new first dimension of
`A`.

Of relevance for types like RGB and BGR, the channels of the returned
array will be in constructor-argument order, not memory order (see
`reinterpret` if you want to use memory order).

The opposite transformation is implemented by `ColorView`.
"""
ChannelView(parent::AbstractArray) = _channelview(parent, channelview_indices(parent))
function _channelview{C<:Colorant,N}(parent::AbstractArray{C}, inds::Indices{N})
    # Creating a ChannelView in a type-stable fashion requires use of tuples to compute N+1
    ChannelView{eltype(C),N,typeof(parent)}(parent)
end

typealias Color1Array{C<:Color1,N} AbstractArray{C,N}
typealias ChannelView1{T,N,A<:Color1Array} ChannelView{T,N,A}

Base.parent(A::ChannelView) = A.parent
parenttype{T,N,A}(::Type{ChannelView{T,N,A}}) = A
@inline Base.size(A::ChannelView)    = channelview_size(parent(A))
@inline Base.indices(A::ChannelView) = channelview_indices(parent(A))

# Can be LinearFast for grayscale (1-channel images), otherwise must be LinearSlow
@pure Base.linearindexing{T<:ChannelView1}(::Type{T}) = Base.linearindexing(parenttype(T))

# colortype(A::ChannelView) = eltype(parent(A))

Base.@propagate_inbounds function Base.getindex{T,N}(A::ChannelView{T,N}, I::Vararg{Int,N})
    @boundscheck checkbounds(A, I...)
    P = parent(A)
    Ic, Ia = indexsplit(P, I)
    @inbounds ret = tuplify(P[Ia...])[Ic]
    ret
end

Base.@propagate_inbounds function Base.setindex!{T,N}(A::ChannelView{T,N}, val, I::Vararg{Int,N})
    @boundscheck checkbounds(A, I...)
    P = parent(A)
    Ic, Ia = indexsplit(P, I)
    @inbounds c = P[Ia...]
    @inbounds P[Ia...] = setchannel(c, val, Ic)
    val
end

function Base.similar{S,N}(A::ChannelView, ::Type{S}, dims::NTuple{N,Int})
    P = parent(A)
    check_ncolorchan(P, dims)
    ChannelView(similar(P, base_colorant_type(eltype(P)){S}, chanparentsize(P, dims)))
end

## ColorView

immutable ColorView{C<:Colorant,N,A<:AbstractArray} <: AbstractArray{C,N}
    parent::A

    function ColorView{T<:Number}(parent::AbstractArray{T})
        n = length(colorviewsize(C, parent))
        n == N || throw(DimensionMismatch("for an $N-dimensional ColorView with color type $C, input dimensionality should be $n instead of $(ndims(parent))"))
        checkdim1(C, size(parent))
        new(parent)
    end
end

"""
    ColorView{C}(A)

creates a "view" of the numeric array `A`, interpreting the first
dimension of `A` as if were the channels of a Colorant `C`. The first
dimension must have the proper number of elements for the constructor
of `C`. For example, if `A` is a 3-by-m-by-n U8 array,
`ColorView{RGB}(A)` will create an m-by-n array with element type
`RGB{U8}`. Color spaces with a single element (i.e., grayscale) do not
"consume" the first dimension of `A`.

Of relevance for types like RGB and BGR, the elements of `A`
are interpreted in constructor-argument order, not memory order (see
`reinterpret` if you want to use memory order).

The opposite transformation is implemented by `ChannelView`.
"""
function (::Type{ColorView{C}}){C<:Colorant,T<:Number}(parent::AbstractArray{T})
    # Creating a ColorView in a type-stable fashion requires use of tuples to compute N+1
    _colorview(base_colorant_type(C){T}, parent, colorviewsize(C, parent))
end
function _colorview{C,N}(::Type{C}, parent::AbstractArray, sz::NTuple{N,Int})
    ColorView{C,N,typeof(parent)}(parent)
end

ColorView(::AbstractArray) = error("specify the desired colorspace with ColorView{C}(parent)")

Base.parent(A::ColorView) = A.parent
parenttype{T,N,A}(::Type{ColorView{T,N,A}}) = A
@inline Base.size(A::ColorView) = colorviewsize(eltype(A), parent(A))

@pure Base.linearindexing{C<:Color1,N,A<:AbstractArray}(::Type{ColorView{C,N,A}}) = Base.linearindexing(A)
@pure Base.linearindexing{V<:ColorView}(::Type{V}) = Base.LinearSlow()

Base.@propagate_inbounds function Base.getindex{C,N}(A::ColorView{C,N}, I::Vararg{Int,N})
    P = parent(A)
    @boundscheck Base.checkbounds_indices(Bool, parentindices(C, indices(P)), I) || Base.throw_boundserror(A, I)
    @inbounds ret = C(getchannels(P, C, I)...)
    ret
end

Base.@propagate_inbounds function Base.setindex!{C,N}(A::ColorView{C,N}, val::C, I::Vararg{Int,N})
    P = parent(A)
    @boundscheck Base.checkbounds_indices(Bool, parentindices(C, indices(P)), I) || Base.throw_boundserror(A, I)
    setchannels!(P, val, I)
    val
end
Base.@propagate_inbounds function Base.setindex!{C,N}(A::ColorView{C,N}, val, I::Vararg{Int,N})
    setindex!(A, convert(C, val), I...)
end

# A grayscale ColorView can be LinearFast, so support this too
Base.@propagate_inbounds function Base.getindex{C<:Color1,N}(A::ColorView{C,N}, i::Int)
    P = parent(A)
    @boundscheck checkindex(Bool, linearindices(P), i) || Base.throw_boundserror(A, i)
    @inbounds ret = C(getchannels(P, C, i)[1])
    ret
end
Base.@propagate_inbounds function Base.setindex!{C<:Color1,N}(A::ColorView{C,N}, val::C, i::Int)
    P = parent(A)
    @boundscheck checkindex(Bool, linearindices(P), i) || Base.throw_boundserror(A, i)
    setchannels!(P, val, i)
    val
end
Base.@propagate_inbounds function Base.setindex!{C<:Color1,N}(A::ColorView{C,N}, val, i::Int)
    setindex!(A, convert(C, val), i)
end

function Base.similar{S<:Colorant,N}(A::ColorView, ::Type{S}, dims::NTuple{N,Int})
    P = parent(A)
    ColorView{S}(similar(P, celtype(eltype(S), eltype(P)), colparentsize(S, dims)))
end
function Base.similar{S<:Number,N}(A::ColorView, ::Type{S}, dims::NTuple{N,Int})
    P = parent(A)
    similar(P, S, dims)
end

## Construct a view that's conceptually equivalent to a ChannelView or ColorView,
## but which may be simpler (i.e., strip off a wrapper or use reinterpret)

"""
    channelview(A)

returns a view of `A`, splitting out (if necessary) the color channels
of `A` into a new first dimension. This is almost identical to
`ChannelView(A)`, except that if `A` is a `ColorView`, it will simply
return the parent of `A`, or will use `reinterpret` when appropriate.
Consequently, the output may not be a `ChannelView` array.
"""
channelview(A::AbstractArray) = ChannelView(A)
channelview(A::ColorView) = parent(A)
channelview{T}(A::Array{RGB{T}}) = reinterpret(T, A)
channelview{C<:AbstractRGB}(A::Array{C}) = ChannelView(A) # BGR, RGB1, etc don't satisfy conditions
channelview{C<:Color}(A::Array{C}) = reinterpret(eltype(C), A)

"""
    colorview(C, A)

returns a view of the numeric array `A`, interpreting successive
elements of `A` as if they were channels of Colorant `C`. This is
almost identical to `ColorView{C}(A)`, except that if `A` is a
`ChannelView`, it will simply return the parent of `A`, or use
`reinterpret` when appropriate. Consequently, the output may not be a
`ColorView` array.

# Example
```jl
A = rand(3, 10, 10)
img = colorview(RGB, A)
```
"""
colorview{C<:Colorant}(::Type{C}, A::AbstractArray) = ColorView{C}(A)
colorview{C<:RGB}(::Type{C}, A::Array) = reinterpret(C, A)
colorview{C<:AbstractRGB}(::Type{C}, A::Array) = ColorView{C}(A)
colorview{C<:Color}(::Type{C}, A::Array) = reinterpret(C, A)
function colorview{C<:Colorant}(::Type{C}, A::ChannelView)
    P = parent(A)
    _colorview(base_colorant_type(C), base_colorant_type(eltype(P)), P)
end
_colorview{C<:Colorant}(::Type{C}, ::Type{C}, A::AbstractArray) = A
_colorview(::Type, ::Type, A::AbstractArray) = throw(ArgumentError("changing colorspaces is not supported"))

## Tuple & indexing utilities
# color->number
@inline channelview_size{C<:Colorant}(parent::AbstractArray{C}) = (length(C), size(parent)...)
@inline channelview_indices{C<:Colorant}(parent::AbstractArray{C}) =
    _cvi(Base.OneTo(length(C)), indices(parent))
_cvi(rc, ::Tuple{}) = (rc,)
_cvi{R<:AbstractUnitRange}(rc, inds::Tuple{R,Vararg{R}}) = (convert(R, rc), inds...)
if squeeze1
    @inline channelview_size{C<:Color1}(parent::AbstractArray{C}) = size(parent)
    @inline channelview_indices{C<:Color1}(parent::AbstractArray{C}) = indices(parent)
end

function check_ncolorchan{C<:Colorant}(::AbstractArray{C}, dims)
    dims[1] == length(C) || throw(DimensionMismatch("new array has $(dims[1]) color channels, must have $(length(C))"))
end
chanparentsize{C<:Colorant}(::AbstractArray{C}, dims) = tail(dims)
@inline colparentsize{C<:Colorant}(::Type{C}, dims) = (length(C), dims...)
if squeeze1
    check_ncolorchan{C<:Color1}(::AbstractArray{C}, dims) = nothing
    chanparentsize{C<:Color1}(::AbstractArray{C}, dims) = dims
    colparentsize{C<:Color1}(::Type{C}, dims) = dims
end

@inline indexsplit{C<:Colorant}(A::AbstractArray{C}, I) = I[1], tail(I)

if squeeze1
    @inline indexsplit{C<:Color1}(A::AbstractArray{C}, I) = 1, I
end

# number->color
@inline colorviewsize{C<:Colorant}(::Type{C}, parent::AbstractArray) = tail(size(parent))
if squeeze1
    @inline colorviewsize{C<:Color1}(::Type{C}, parent::AbstractArray) = size(parent)
end

function checkdim1{C<:Colorant}(::Type{C}, dims)
    dims[1] == length(C) || throw(DimensionMismatch("dimension 1 must have size $(length(C))"))
    nothing
end
if squeeze1
    checkdim1{C<:Color1}(::Type{C}, dims) = nothing
end

parentindices(::Type, inds) = tail(inds)
if squeeze1
    parentindices{C<:Color1}(::Type{C}, inds) = inds
end

celtype{T}(::Type{Any}, ::Type{T}) = T
celtype{T1,T2}(::Type{T1}, ::Type{T2}) = T1

## Low-level color utilities

tuplify(c::Color1) = (comp1(c),)
tuplify(c::Color3) = (comp1(c), comp2(c), comp3(c))
tuplify(c::Color2) = (comp1(c), alpha(c))
tuplify(c::Color4) = (comp1(c), comp2(c), comp3(c), alpha(c))

"""
    getchannels(P, C::Type, I)

Get a tuple of all channels needed to construct a Colorant of type `C`
from an `P::AbstractArray{<:Number}`.
"""
getchannels
if squeeze1
    @inline getchannels{C<:Color1}(P, ::Type{C}, I) = (@inbounds ret = (P[I...],); ret)
    @inline getchannels{C<:Color1}(P, ::Type{C}, I::Real) = (@inbounds ret = (P[I],); ret)
else
    @inline getchannels{C<:Color1}(P, ::Type{C}, I) = (@inbounds ret = (P[1, I...],); ret)
    @inline getchannels{C<:Color1}(P, ::Type{C}, I::Real) = (@inbounds ret = (P[1, I],); ret)
end
@inline function getchannels{C<:Color2}(P, ::Type{C}, I)
    @inbounds ret = (P[1,I...], P[2,I...])
    ret
end
@inline function getchannels{C<:Color3}(P, ::Type{C}, I)
    @inbounds ret = (P[1,I...], P[2,I...],P[3,I...])
    ret
end
@inline function getchannels{C<:Color4}(P, ::Type{C}, I)
    @inbounds ret = (P[1,I...], P[2,I...], P[3,I...], P[4,I...])
    ret
end

# setchannel (similar to setfield!)
# These don't check bounds since that's already done
"""
    setchannel(c, val, idx)

Equivalent to:

    cc = copy(c)
    cc[idx] = val
    cc

for immutable colors. `idx` is interpreted in the sense of constructor
arguments, so `setchannel(c, 0.5, 1)` would set red color channel for
any `c::AbstractRGB`, even if red isn't the first field in the type.
"""
setchannel{T}(c::Colorant{T,1}, val, Ic::Int) = typeof(c)(val)

setchannel{C,T}(c::TransparentColor{C,T,2}, val, Ic::Int) =
    typeof(c)(ifelse(Ic==1,val,comp1(c)),
              ifelse(Ic==2,val,alpha(c)))

setchannel{T}(c::Colorant{T,3}, val, Ic::Int) = typeof(c)(ifelse(Ic==1,val,comp1(c)),
                                                          ifelse(Ic==2,val,comp2(c)),
                                                          ifelse(Ic==3,val,comp3(c)))
setchannel{C,T}(c::TransparentColor{C,T,4}, val, Ic::Int) =
    typeof(c)(ifelse(Ic==1,val,comp1(c)),
              ifelse(Ic==2,val,comp2(c)),
              ifelse(Ic==3,val,comp3(c)),
              ifelse(Ic==4,val,alpha(c)))

"""
    setchannels!(P, val, I)

For a color `val`, distribute its channels along `P[:, I...]` for
`P::AbstractArray{<:Number}`.
"""
setchannels!
if squeeze1
    @inline setchannels!(P, val::Color1, I) = (@inbounds P[I...] = comp1(val); val)
else
    @inline setchannels!(P, val::Color1, I) = (@inbounds P[1,I...] = comp1(val); val)
end
@inline function setchannels!(P, val::Color2, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = alpha(val)
    val
end
@inline function setchannels!(P, val::Color3, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = comp2(val)
    @inbounds P[3,I...] = comp3(val)
    val
end
@inline function setchannels!(P, val::Color4, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = comp2(val)
    @inbounds P[3,I...] = comp3(val)
    @inbounds P[4,I...] = alpha(val)
    val
end


text:
__precompile__(true)
module Roots

import Base: *

using Polynomials
import Polynomials: roots

using ForwardDiff
using Compat

export roots

export fzero,
       fzeros,
       newton, halley,
       secant_method, steffensen,
       multroot,
       D, D2

## load in files
include("fzero.jl")
include("adiff.jl")
include("SOLVE.jl")
include("derivative_free.jl")
include("newton.jl")
include("Polys/polynomials.jl")
include("Polys/agcd.jl")
include("Polys/multroot.jl")
include("Polys/real_roots.jl")

## Main functions are
## * fzero(f, ...) to find _a_ zero of f
## * fzeros(f, ...) to attempt to find all zeros of f

"""

Find zero of a function using an iterative algorithm

* `f`: a scalar function or callable object
* `x0`: an initial guess, finite valued.

Keyword arguments:

* `ftol`: tolerance for a guess `abs(f(x)) < ftol`
* `xtol`: stop if `abs(xold - xnew) <= xtol + max(1, |xnew|)*xtolrel`
* `xtolrel`: see `xtol`
* `maxeval`: maximum number of steps
* `verbose`: Boolean. Set `true` to trace algorithm
* `order`: Can specify order of algorithm. 0 is most robust, also 1, 2, 5, 8, 16.
* `kwargs...` passed on to different algorithms. There are `maxfneval` when `order` is 1,2,5,8, or 16 and `beta` for orders 2,5,8,16,

This is a polyalgorithm redirecting different algorithms based on the value of `order`.

"""
function fzero(f, x0::Real; kwargs...)
    x = float(x0)
    isinf(x) && throw(ConvergenceFailed("An initial value must be finite"))
    derivative_free(f, x; kwargs...)
end

"""
Find zero of a function within a bracket

Uses a modified bisection method for non `big` arguments

Arguments:

* `f` A scalar function or callable object
* `a` left endpont of interval
* `b` right endpont of interval
* `xtol` optional additional tolerance on sub-bracket size.

For a bracket to be valid, it must be that `f(a)*f(b) < 0`.

For `Float64` values, the answer is such that `f(prevfloat(x)) *
f(nextfloat(x)) < 0` unless a non-zero value of `xtol` is specified in
which case, it stops when the sub-bracketing produces an bracket with
length less than `max(xtol, abs(x1)*xtolrel)`.

For `Big` values, which defaults to the algorithm of Alefeld, Potra, and Shi, a
default tolerance is used for the sub-bracket length that can be enlarged with `xtol`.

If `a==-Inf` it is replaced with `nextfloat(a)`; if `b==Inf` it is replaced with `prevfloat(b)`.

Example:

    `fzero(sin, 3, 4)` # find pi
    `fzero(sin, [big(3), 4]) find pi with more digits
"""
function fzero(f, a::Real, b::Real; kwargs...)
    a,b = sort([a,b])
    a,b = promote(float(a), b)
    (a == -Inf) && (a = nextfloat(a))
    (b == Inf) && (b = prevfloat(b))
    (isinf(a) | isinf(b)) && throw(ConvergenceFailed("A bracketing interval must be bounded"))
    find_zero(f,a,b; kwargs...)
end

"""
Find a zero with bracket specified via `[a,b]`, as `fzero(sin, [3,4])`.
"""
function fzero{T <: Real}(f, bracket::Vector{T}; kwargs...)
    fzero(f, float(bracket[1]), float(bracket[2]); kwargs...)
end

"""
Find a zero within a bracket with an initial guess to *possibly* speed things along.
"""
function fzero{T <: Real}(f, x0::Real, bracket::Vector{T}; kwargs...)
    a,b = sort(map(float,bracket))
    (a == -Inf) && (a = nextfloat(a))
    (b == Inf) && (b = prevfloat(b))
    (isinf(a) | isinf(b)) && throw(ConvergenceFailed("A bracketing interval must be bounded"))
    try
        ex = a42a(f, a, b, float(x0); kwargs...)
    catch ex
        if isa(ex, StateConverged)
            return(ex.x0)
        else
            rethrow(ex)
        end
    end
end

"""
Find zero using Newton's method.
"""
fzero(f::Function, fp::Function, x0::Real; kwargs...) = newton(f, fp, float(x0); kwargs...)

"""

Find zero of a polynomial with derivative free algorithm.

Arguments:

* `p` a `Poly` object
* `x0` an initial guess

Returns:

An approximate root or an error.

See `fzeros(p)` to return all real roots.

"""
function fzero(p::Poly, x0::Real, args...; kwargs...)
    fzero(convert(Function, p), float(x0), args...; kwargs...)
end

function fzero{T <: Real}(p::Poly, bracket::Vector{T}; kwargs...)
    a, b = float(bracket[1]), float(bracket[2])
    fzero(convert(Function, p), a, b; kwargs...)
end
function fzero{T <: Real}(p::Poly, x0::Real, bracket::Vector{T}; kwargs...)
    fzero(convert(Function,p), float(x0), map(float,bracket); kwargs...)
end

## fzeros

"""

Find real zeros of a polynomial

args:

`f`: a Polynomial function of R -> R. May also be of `Poly` type.

For polynomials in Z[x] or Q[x], the `real_roots` method will use
`Poly{Rational{BigInt}}` arithmetic, which allows it to handle much
larger degree polynomials accurately. This can be called directly
using `Polynomial` objects.

"""
fzeros(p::Poly) = real_roots(p)

function fzeros(f)
    p = poly([0.0])
    try
        p = convert(Poly, f)
    catch e
        error("If f(x) is not a polynomial in x, then an interval to search over is needed")
    end
    zs = fzeros(p)
    ## Output is mixed, integer, rational, big. We tidy up
    etype = eltype(f(0.0))
    ietype = eltype(f(0))
    out = Real[]
    for z in zs
        if isa(z, Rational)
            val = z.den == 1 ? convert(ietype, z.num) : convert(Rational{ietype}, z)
            push!(out, val)
        else
            push!(out, convert(etype, z))
        end
    end
    sort(out)
end

"""

Attempt to find all zeros of `f` within an interval `[a,b]`.

Simple algorithm that splits `[a,b]` into subintervals and checks each
for a root.  For bracketing subintervals, bisection is
used. Otherwise, a derivative-free method is used. If there are a
large number of roots found relative to the number of subintervals, the
number of subintervals is increased and the process is re-run.

There are possible issues with close-by zeros and zeros which do not
cross the origin (non-simple zeros). Answers should be confirmed
graphically, if possible.

"""
function fzeros{T <: Real}(f, bracket::Vector{T}; kwargs...)
    ## check if a poly
    try
        rts = fzeros(f)
        filter(x -> bracket[1] <= x <= bracket[2], rts)
    catch e
        find_zeros(f, float(bracket[1]), float(bracket[2]); kwargs...)
    end
end
fzeros(f::Function, a::Real, b::Real; kwargs...) = fzeros(f, [a,b]; kwargs...)

"""

Factor a polynomial function.

Finds factors numerically.

For polynomial functions over the integers or rational tries -- as naively as possible -- to factor exactly over the rationals first.

Returns a dictionary with keys that are roots and values that are multiplicities

Example:
```
factor(x -> (x-1)^3*(x-2))
```

"""
function Base.factor(f::Function)
    p = poly([0.0])
    try
        p = convert(Poly, f)
    catch e
        throw(DomainError()) # `factor` only works with polynomial functions"
    end
    factor(p)
end

end

ext: jl


text:
"""
    struct TensorKeyIterator{I<:Sector, F₁<:FusionTree{I}, F₂<:FusionTree{I}}
        rowr::SectorDict{I, FusionTreeDict{F₁, UnitRange{Int}}}
        colr::SectorDict{I, FusionTreeDict{F₂, UnitRange{Int}}}
    end

An iterator over all different splitting and fusion tree pairs.
"""
struct TensorKeyIterator{I<:Sector, F₁<:FusionTree{I}, F₂<:FusionTree{I}}
    rowr::SectorDict{I, FusionTreeDict{F₁, UnitRange{Int}}}
    colr::SectorDict{I, FusionTreeDict{F₂, UnitRange{Int}}}
end

"""
    struct TensorPairIterator{I<:Sector, F₁<:FusionTree{I}, F₂<:FusionTree{I},A<:DenseMatrix}
        rowr::SectorDict{I, FusionTreeDict{F₁, UnitRange{Int}}}
        colr::SectorDict{I, FusionTreeDict{F₂, UnitRange{Int}}}
        data::SectorDict{I, A}
    end

An iterator that has not been implemented.
"""
struct TensorPairIterator{I<:Sector, F₁<:FusionTree{I}, F₂<:FusionTree{I}, A<:DenseMatrix}
    rowr::SectorDict{I, FusionTreeDict{F₁, UnitRange{Int}}}
    colr::SectorDict{I, FusionTreeDict{F₂, UnitRange{Int}}}
    data::SectorDict{I, A}
end

const TensorIterator{I<:Sector, F₁<:FusionTree{I}, F₂<:FusionTree{I}} =
        Union{TensorKeyIterator{I, F₁, F₂}, TensorPairIterator{I, F₁, F₂}}

Base.IteratorSize(::Type{<:TensorIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:TensorIterator}) = Base.HasEltype()
Base.eltype(T::Type{TensorKeyIterator{I, F₁, F₂}}) where {I, F₁, F₂} = Tuple{F₁, F₂}

function Base.length(t::TensorKeyIterator)
    l = 0
    for (rowdict, coldict) in zip(values(t.rowr), values(t.colr))
        l += length(rowdict) * length(coldict)
    end
    return l
end

function Base.iterate(it::TensorKeyIterator)
    i = 1
    i > length(it.rowr) && return nothing
    rowit, colit = it.rowr.values[i], it.colr.values[i]

    rownext = iterate(rowit)
    colnext = iterate(colit)
    while true
        if rownext === nothing
            i += 1
        elseif colnext === nothing
            i += 1
        else
            break
        end
        i > length(it.rowr) && return nothing
        rowit, colit = it.rowr.values[i], it.colr.values[i]
        rownext = iterate(rowit)
        colnext = iterate(colit)
    end
    (f1, r1), rowstate = rownext
    (f2, r2), colstate = colnext

    return (f1, f2), (f2, i, rowstate, colstate)
end

function Base.iterate(it::TensorKeyIterator, state)
    (f2, i, rowstate, colstate) = state
    rowit, colit = it.rowr.values[i], it.colr.values[i]
    rownext = iterate(rowit, rowstate)
    if rownext !== nothing
        (f1, r1), rowstate = rownext
        return (f1, f2), (f2, i, rowstate, colstate)
    end
    colnext = iterate(colit, colstate)
    if colnext !== nothing
        rownext = iterate(rowit) # should not be nothing
        @assert rownext !== nothing
        (f1, r1), rowstate = rownext
        (f2, r2), colstate = colnext
        return (f1, f2), (f2, i, rowstate, colstate)
    end
    while true
        if rownext === nothing
            i += 1
        elseif colnext === nothing
            i += 1
        else
            break
        end
        i > length(it.rowr) && return nothing
        rowit, colit = it.rowr.values[i], it.colr.values[i]
        rownext = iterate(rowit)
        colnext = iterate(colit)
    end
    (f1, r1), rowstate = rownext
    (f2, r2), colstate = colnext

    return (f1, f2), (f2, i, rowstate, colstate)
end


text:
@symbol_func function G_CO(cur_reactor::AbstractReactor)
  cur_G = safe_symbol(cur_reactor, :B_0)

  cur_G /= (2/5) * cur_reactor.pi

  cur_G /= cur_reactor.J_TF

  cur_G /= ( 1 - epsilon_b(cur_reactor) )

  cur_G
end


text:
using ImageFiltering

function lbp_origin_mapwindow(X)
    function _lbp(block)
        gc = block[2, 2]
        offsets = ((-1, -1), (0, -1), (1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1))

        rst = 0
        for i in 1:length(offsets)
            o = offsets[i]
            rst += ifelse(gc <= block[2+o[1], 2+o[2]], 1, 0) << (i-1)
        end
        return rst
    end
    mapwindow(_lbp, X, (3, 3))
end


text:
#
#	一些比较常用的函数
#

#============================================
utf-8变长编码-> 可直接索引的[]存储?
  使用方式:
	substr("阿三科技的符号i阿斯顿飞呢",3,-2)
	substr("阿三科技的符号i阿斯顿飞呢",3,7)
========================================#

function substr(str::String,start_idx::Int,end_idx::Int)
	#半开半闭区间[start,end) 没法表示索引到结尾
	#使用[start,end] 更附和julia的原则:从1开始
	if end_idx >length(str) @error "BoundsError 索引越界";return end
	char_list=[];
	tmp=""

	for item in str		push!(char_list,item)	end
	if end_idx<0 	end_idx=(length(char_list)+end_idx+1) end
	if start_idx>end_idx	return "" end   #排除错误参数
	for item in char_list[start_idx:end_idx]
		tmp=tmp*item
	end
	#@info "$(tmp)"
	return tmp
end

#=======================================
	类型转换:  Int64('苏') ->33487
		String -> Int32 -> bitstring
		Int->Float32 ->
====================================#
function str2array(str::String)
	tmp=[]
	for item in str
		push!(tmp,item)
	end
	return tmp
end
function str2int_array(str::String)
	tmp=Array{Int32,1}()  #不要忘了维度信息
	for item in str
		push!(tmp,Int32(item))
	end
	return tmp
end
function str2int_str(str::String)
	tmp=""
	for item in str
		tmp=tmp*repr(Int(item))
	end
	return tmp
end
function str2BigInt(str::String)
	#十进制 str->BigInt转换
	table=Set(['0','1','2','3','4','5','6','7','8','9'])
	for item in str
		@assert item in table "$(item)不是十进制数字"
	end
	return parse(BigInt,str)
	#return big"1234645675643523412423123456457"
end

#using PyPlot
function 画图()
	data1=[1.75,2.35,3.0]
	sina_debt=[2.449,2.668,2.659,2.834,3.061,3.178,3.743]
	usa_debt =[1.508,1.587,1.601,1.678,1.800,1.877,2.318]
	japan_debt=[-0.097,-0.110,-0.135,-0.117,-0.1,-0.006,0.424]
	de_debt=[-0.692,-0.641,-0.612,-0.531,-0.42,-0.253,0.276]
	en_debt=[0.597,0.535,0.493,0.564,0.558,0.754,1.256]
	title("2019-12-29各国债务收益率曲线")
	xlabel("时间（年）")
	ylabel("年率")
	xticks=[1,2,3,5,7,10,30]
	x=[1,2,3,5,7,10,30]
	text(20,3.5,"中国国债")
	text(20,2,"美国国债")
	text(20,0.3,"日本国债")
	text(10,-0.5,"德国国债")
	text(20,1,"英国国债")

	plot([1,2,3],data1,label="中国利息")
	plot(x,sina_debt,"or-")
	plot(x,usa_debt,"^k-")
	plot(x,japan_debt,"b^--")
	plot(x,de_debt,"yo--")
	plot(x,en_debt,"go-")

	grid() #显示网格
	legend()
	show()
end
#=====================================

=====================================#
function 鸡蛋价格(总价,重量=45)
	单价每斤=总价/重量;
	十元重量=10/单价每斤
	十元售重=十元重量-0.5
	@info "鸡蛋价 : " 单价每斤 十元重量 十元售重
end
#鸡蛋价格(175)
#=====================================

=====================================#
function 画圆()

	for 间隔=1:50
		x=[i for i=0:间隔:2000pi];
		a=sin.(x);
		b=cos.(x);
		plot(a,b,"bo--");
		title("x=[i for i=0:"*string(间隔)*":2000pi];a=sin.(x);b=cos.(x)")
		show()
	end
end
#画圆()
#========================================

			主函数

=======================================#

function test()

	画图()

	return true

	substr("阿三科技的符号i阿斯顿飞呢",3,-2)
	substr("阿三科技的符号i阿斯顿飞呢",3,7)

	asd="阿三健康的符合牛asdf郎我烦恼"
	@show str2array(asd)
	@show str2int_array(asd)
	@show str2int_str(asd)
	@show str2BigInt(str2int_str(asd))
	@show str2BigInt("123343245e")
end
#test()

#计算杨辉三角数列 1 1开始
#保加利亚单人纸牌游戏（Bulgarian solitaire）
#取出 45 张牌，然后把它们随意分成若干堆。接下来，从每一堆里各取一张牌，叠在一起形成一堆新的牌。不断这样做下去，如果某个时候桌面上正好有 9 堆牌，并且各堆牌数分别为 1, 2, 3, 4, …, 9 ，你就获胜了

max_stars_repo_name: Jonson-Sun/C-plus-plus

max_issues_repo_name: Jonson-Sun/C-plus-plus

max_forks_repo_name: Jonson-Sun/C-plus-plus


text:
"""
    Simply cuts the elements up into n_procs chunks
    and assigns them in increasing order based on the
    element IDs

    Note: This is generally a REALLY bad idea!
    (but it will get us going)
"""
abstract MetisPartitioner <: Partitioner

const metis_library = string(ENV["PETSC_DIR"],"/lib/libmetis")

const METIS_NOPTIONS = 40
typealias idx_t Int32

# One based
const METIS_OPTION_PTYPE = 1
const METIS_OPTION_OBJTYPE = 2
const METIS_OPTION_CTYPE = 3
const METIS_OPTION_IPTYPE = 4
const METIS_OPTION_RTYPE = 5
const METIS_OPTION_DBGLVL = 6
const METIS_OPTION_NITER = 7
const METIS_OPTION_NCUTS = 8
const METIS_OPTION_SEED = 9
const METIS_OPTION_NO2HOP = 10
const METIS_OPTION_MINCONN = 11
const METIS_OPTION_CONTIG = 12
const METIS_OPTION_COMPRESS = 13
const METIS_OPTION_CCORDER = 14
const METIS_OPTION_PFACTOR = 15
const METIS_OPTION_NSEPS = 16
const METIS_OPTION_UFACTOR = 17
const METIS_OPTION_NUMBERING = 18

const METIS_PTYPE_RB = 0
const METIS_PTYPE_KWAY = 1

const METIS_OBJTYPE_CUT = 0
const METIS_OBJTYPE_VOL = 1
const METIS_OBJTYPE_NODE = 2

const METIS_IPTYPE_GROW = 0
const METIS_IPTYPE_RANDOM = 1
const METIS_IPTYPE_EDGE = 2
const METIS_IPTYPE_NODE = 3
const METIS_IPTYPE_METISRB = 4

const METIS_CTYPE_RM = 0
const METIS_CTYPE_SHEM = 1

function partition!(mesh::Mesh, ::Type{MetisPartitioner})
    n_procs = MPI.Comm_size(MPI.COMM_WORLD)
    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    options = Array{idx_t}(METIS_NOPTIONS)

    ccall((:METIS_SetDefaultOptions, metis_library), Int32, (Ref{idx_t},), options)

#    options[METIS_OPTION_NUMBERING] = 1 # 1-based numbering! Handy!  Well... I ended up doing 0-based because debugging (and the answer needs to be zero based anyway because MPI is)

    # libMesh defaults
    if n_procs <= 8
        options[METIS_OPTION_PTYPE] = METIS_PTYPE_RB
    else
        options[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY
    end

#    options[METIS_OPTION_IPTYPE] = METIS_IPTYPE_GROW
#    options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT
#    options[METIS_OPTION_CTYPE] = METIS_CTYPE_RM
#    options[METIS_OPTION_NCUTS] = 4
#    options[METIS_OPTION_NITER] = 10
#    options[METIS_OPTION_MINCONN] = 0
#    options[METIS_OPTION_UFACTOR] = 1
#    options[METIS_OPTION_DBGLVL] = (idx_t)(1) | (idx_t)(16) | (idx_t)(128) | (idx_t)(256)

    ne = (idx_t)(length(mesh.elements))
    nn = (idx_t)(length(mesh.nodes))

    eptr = Array{idx_t}(ne+1)
    eind = Array{idx_t}(ne*4)

    for elem in mesh.elements
        node_start = (4*(elem.id-1))
        eptr[elem.id] = node_start

        for n in 1:length(elem.nodes)
            eind[node_start + n] = elem.nodes[n].id - 1
        end
    end

    eptr[length(eptr)] = (ne*4)

    # Note!  Weights are required!  Right now it's just 4 (the number of nodes)
    vwgt = Array{idx_t}(ne)
    fill!(vwgt, 4)

    ncommon = (idx_t)(2)
    nparts = (idx_t)(n_procs)
    objval = (idx_t)(0)

    epart = Array{idx_t}(ne)
    npart = Array{idx_t}(nn)

    ret_val = ccall((:METIS_PartMeshDual, metis_library), Int32, (Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ptr{Float64}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}),
                    ne, nn, eptr, eind, vwgt, C_NULL, ncommon, nparts, C_NULL, options, objval, epart, npart)

    # Assign the IDs
    for elem in mesh.elements
        elem.processor_id = epart[elem.id]
    end

#    ret_val = ccall((:METIS_PartMeshNodal, metis_library), Int32, (Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ptr{Float64}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}),
#                    ne, nn, eptr, eind, C_NULL, C_NULL, nparts, C_NULL, options, objval, epart, npart)

    # numflag = 0

    # xadj = Ref{Ptr{idx_t}}()
    # adjncy = Ref{Ptr{idx_t}}()

    # ccall((:METIS_MeshToDual, metis_library), Int32, (Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{Ptr{idx_t}}, Ref{Ptr{idx_t}}),
    #       ne, nn, eptr, eind, ncommon, numflag, xadj, adjncy)

    # xadj_array = unsafe_wrap(Array, xadj[], ne + 1, true)
    # adjncy_array = unsafe_wrap(Array, adjncy[], xadj_array[end], true)

    # println(xadj_array)
    # println(adjncy_array)

    # nvtxs = length(xadj_array) - 1
    # ncon = 1

    # part = Array{idx_t}(nvtxs)

    # vwgt = Array{idx_t}(nvtxs)
    # fill!(vwgt, 4)

    # ret_val = ccall((:METIS_PartGraphRecursive, metis_library), Int32, (Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}, Ptr{Float64}, Ptr{Float64}, Ref{idx_t}, Ref{idx_t}, Ref{idx_t}),
    #                 nvtxs, ncon, xadj_array, adjncy_array, vwgt, C_NULL, C_NULL, nparts, C_NULL, C_NULL, options, objval, part)

    # println(ret_val)
    # println(part)

    assignNodeProcessorIDs(mesh)
end


text:
mutable struct PINFO
    CPolar::Point2D{Polar}
    CCart::Point2D{Cartesian}
    UCart::Point2D{Cartesian}
    t::Int64
    t_parent::Int64
    color::Vector{Float32}
end

# mutable struct PINFO{A,B,C,D}
#     CPolar::A
#     CCart::B
#     UCart::B
#     t::C
#     t_parent::C
#     color::D
# end

struct PWEIGHTS # consider making it immutable
    barycentric::MArray{Tuple{3},Float64,1,3}
    ip_weights::MArray{Tuple{6},Float64,1,6}
    node_weights::MArray{Tuple{3},Float64,1,3}
end

# struct PWEIGHTS{T,M} # consider making it immutable
#     barycentric::T
#     ip_weights::M
#     node_weights::T
# end

mutable struct PVAR{M<:AbstractArray}
    T::M
    ΔT_subgrid::M
    Fxx::M
    Fzz::M
    Fxz::M
    Fzx::M
end

init_pinfo(np) = [PINFO(Point2D{Polar}(0.0,0.0),
                  Point2D{Cartesian}(0.0,0.0),
                  Point2D{Cartesian}(0.0,0.0),
                  Int32(0),
                  Int32(0),
                  zeros(3)
                  ) for _ in 1:np]

init_pweights(np) = [PWEIGHTS(
                    (@SVector zeros(3)),
                    (@SVector zeros(6)),
                    (@SVector zeros(3))
                    ) for _ in 1:np]

init_pvars(np) = PVAR(Vector{Float64}(undef,np),
                      fill(1.0,np),
                      fill(1.0,np),
                      fill(1.0,np),
                      fill(0.0,np),
                      fill(0.0,np))

@inline function getvertices!(v,x,i0)
    @inbounds for ii = 1:3
        id = @views i0[ii]
        v[ii] = x[id]
    end
end

function getvertices!(v::Array{Point2D{Polar},1},θ,r,idx)
    @inbounds for ii = 1:3
        v[ii].x = θ[ii,idx]
        v[ii].z = r[ii,idx]
    end
end

function getvertices!(v::Array{Point2D{Cartesian},1},x,z,idx)
    @inbounds for ii = 1:3
        v[ii].x = x[ii,idx]
        v[ii].z = z[ii,idx]
    end
end

function getvertices!(v::Point2D{Polar},x,z,idx)
    @inbounds for ii = 1:3
        v[ii].x = x[ii,idx]
        v[ii].z = z[ii,idx]
    end
end

@inline getcartesian(A::PINFO) = A.CCart

@inline getpolar(A::PINFO) = A.CPolar

@inline polar2cartesian(p::Point2D{Polar}) = Point2D{Cartesian}(p.z*sin(p.x), p.z*cos(p.x))

@inline polar2x(p::Point2D{Polar}) = p.z*sin(p.x)
@inline polar2z(p::Point2D{Polar}) = p.z*cos(p.x)

@inline cartesian2θ(A::Point2D{Cartesian}) = atan(A.x, A.z)
@inline cartesian2r(A::Point2D{Cartesian}) = sqrt(A.x^2 + A.z^2)

@inline cartesian2polar(p::Point2D{Cartesian}) = Point2D{Polar}(atan(p.x, p.z), sqrt(p.x^2+p.z^2))

@inline cartesian2polar(x::AbstractArray, z::AbstractArray) = (atan.(x, z), @. sqrt(x^2+z^2))

@inline polar2cartesian(θ::AbstractArray, r::AbstractArray) = (@. r*sin(θ), @. r*cos(θ))

function cartesian2polar!(particle_info::Array{PINFO,1})
    Threads.@threads for i in axes(particle_info,1)
        particle_info[i].CPolar.x = atan(particle_info[i].CCart.x, particle_info[i].CCart.z)
        particle_info[i].CPolar.z = sqrt(particle_info[i].CCart.x^2 + particle_info[i].CCart.z^2)

        # Fix θ
        if particle_info[i].CPolar.x < 0
            particle_info[i].CPolar.x += 2π
        elseif particle_info[i].CPolar.x > 2π
            particle_info[i].CPolar.x -= 2π
        end

        # Fix radius
        if particle_info[i].CPolar.z < 1.22
            particle_info[i].CPolar.z = 1.22001
            particle_info[i].CCart.x = particle_info[i].CPolar.z * sin(particle_info[i].CPolar.x)
            particle_info[i].CCart.z = particle_info[i].CPolar.z * cos(particle_info[i].CPolar.x)

        elseif particle_info[i].CPolar.z > 2.22
            particle_info[i].CPolar.z = 2.21999
            particle_info[i].CCart.x = particle_info[i].CPolar.z * sin(particle_info[i].CPolar.x)
            particle_info[i].CCart.z = particle_info[i].CPolar.z * cos(particle_info[i].CPolar.x)
        end

    end
end

function polar2cartesian!(particle_info::Array{PINFO,1})
    Threads.@threads for i in axes(particle_info,1)
        zp, xp = particle_info[i].CPolar.z, particle_info[i].CPolar.x
        particle_info[i].CCart.x = zp * sin(xp)
        particle_info[i].CCart.z = zp * cos(xp)
    end
end

@inline function distance(p1::Point2D{Polar},p2::Array{Point2D{Polar},1})
    ld = length(p2)
    d  = Vector{Float64}(undef, ld)
    @inbounds for i in 1:ld
        d[i] = √( p1.z^2 + p2[i].z^2 - 2p1.z*p2[i].z*cos(p1.x-p2[i].x))
    end
    return d
end

@inline distance(p1::Point2D{Cartesian},p2::Point2D{Cartesian}) = √((p1.x - p2.x)^2 + (p1.z - p2.z)^2)

@inline distance(p1::Point2D{Polar},p2::Point2D{Polar}) = distance(polar2cartesian(p1), polar2cartesian(p2))

@inline distance(p1::Point2D{Polar},p2::NamedTuple) = distance(polar2cartesian(p1), polar2cartesian(p2))

@inline distance(p1::Point2D{Cartesian},p2::NamedTuple) = √((p1.x - p2.x)^2 + (p1.z - p2.z)^2)

@inline distance(x1::T, z1::T, x2::T, z2::T) where {T<:Real} = √((x1-x2)^2 + (z1-z2)^2)

@inline distance(p1::NTuple{2,T}, p2::NTuple{2,T}) where {T<:Real} = √((p1[1]-p2[1])^2 + (p1[2]-p2[2])^2)

@inline getcentroid(p1::Point2D{T},p2::Point2D{T},p3::Point2D{T}) where T =
    Point2D{T}((p1.x+p2.x+p3.x)/3,(p1.z+p2.z+p3.z)/3)

@inbounds function barycentric_coordinates(bary_dummy, in, vc, pc)
    x1,x2,x3 = vc[1].x, vc[2].x, vc[3].x
    z1,z2,z3 = vc[1].z, vc[2].z, vc[3].z
    x,z      = pc.x, pc.z
    inside   = false

    bary_dummy[1] = @muladd ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))

    if 0 <= bary_dummy[1] <= 1
        bary_dummy[2] = @muladd  ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
            ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))

        if 0 <= bary_dummy[2] <= 1
            bary_dummy[3] = 1 - bary_dummy[1] - bary_dummy[2]

            0 <= bary_dummy[3] <= 1 ? in = true : in = false
        end

    end

    return bary_dummy, inside

end

@inbounds function barycentric_coordinates2( intri, vc, pc)
    x1,x2,x3 = vc[1].x, vc[2].x, vc[3].x
    z1,z2,z3 = vc[1].z, vc[2].z, vc[3].z
    x,z = pc.x, pc.z
    intri = false

    bc1, bc2, bc3 = 0.0, 0.0, 0.0

    bc1 = @muladd ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))

    if 0 <= bc1 <= 1
        bc2 = @muladd ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
            ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
        if 0 <= bc2 <= 1
            bc3 = 1 - bc1 - bc2
            intri = 0 <= bc3 <= 1 ?  true : false
        end
    end

    bary_dummy = @MVector [bc1,bc2,bc3]

    return bary_dummy, intri

end

@inline function barycentric_coordinates3(vc, pc)
    # Coordinates of triangle's vertices
    x1,x2,x3 = vc[1].x, vc[2].x, vc[3].x
    z1,z2,z3 = vc[1].z, vc[2].z, vc[3].z
    # Coordinates of query point
    x, z = pc.x, pc.z

    bc1 = @muladd ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))

    if 0.0 <= bc1 <= 1.0
        bc2 = @muladd  ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
              ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
        if 0.0 <= bc2 <= 1.0
            bc3 = 1.0 - bc1 - bc2
            0.0 <= bc3 <= 1.0 ? intri = true : intri = false
        end
    end

    return intri

end

function check_corruption!(found, particle_fields)
    np = length(particle_fields.T)
    Threads.@threads for i in 1:np
        @inbounds if isinf(particle_fields.T[i]) ||
           isinf(particle_fields.Fxx[i]) ||
           isinf(particle_fields.Fxz[i]) ||
           isinf(particle_fields.Fzx[i]) ||
           isinf(particle_fields.Fzz[i]) ||
           isnan(particle_fields.T[i])   ||
           isnan(particle_fields.Fxx[i]) ||
           isnan(particle_fields.Fxz[i]) ||
           isnan(particle_fields.Fzx[i]) ||
           isnan(particle_fields.Fzz[i])

           found[i] = false
        end
    end
end

function purgeparticles(particle_info, particle_weights, particle_fields, found)
    ikeep = findall(found)
    particle_info = [@inbounds particle_info[i] for i in ikeep]
    particle_weights = [@inbounds particle_weights[i] for i in ikeep]

    @inbounds for i in eachindex(found)
        if ikeep[i] == false
            popat!(particle_fields.T, i)
            popat!(particle_fields.Fxx, i)
            popat!(particle_fields.Fzz, i)
            popat!(particle_fields.Fxz, i)
            popat!(particle_fields.Fzx, i)
        end
    end

    return particle_info, particle_weights,particle_fields
end

function check_crossingπ!(vertices,particle)
    if (particle.x < π)  && (maximum(@SVector [vertices[1].x, vertices[3].x, vertices[3].x]) > 2π)
        for i in 1:3
            @inbounds if vertices[i].x > π
                vertices[i].x -= 2π
            end
        end
    end
end

function check_crossingπ(vertices,xp)
    xp[1] = vertices[1].x
    xp[2] = vertices[2].x
    xp[3] = vertices[3].x
    maxθ = maximum(xp)
    minθ = minimum(xp)

    check = false
    if (maxθ - minθ) > π
        for i in 1:3
            @inbounds if vertices[i].x > π
                vertices[i].x -= 2π
                check = true
            end
        end
    end
    return vertices, check
end

@inline function fillweights!(particle_weights,bary_dummy,particle,vertices,IC,ipart)
    @inbounds @simd for ii = 1:6
        if ii < 4
            # barycentric weights
            particle_weights[ipart].barycentric[ii] = bary_dummy[ii]
            # Nodal weights (for particle 2 node)
            particle_weights[ipart].node_weights[ii] = distance(particle,vertices[ii])
        end
        # Barycentric weights (for particle 2 ip)
        particle_weights[ipart].ip_weights[ii] = distance(particle,IC[ii])
    end
end

@inline function getreps(t0, nel, maxreps)
    reps = fill(Int32(0),nel)
    toremove = Int32[]
    @inbounds for (i,t) in enumerate(t0)
        reps[t] += Int32(1)
        if reps[t] > maxreps
            push!(toremove, i)
        end
    end
    return reps, toremove
end

@inline removeparticles(PC, ikeep) = [PC[i] for i in ikeep]

@inline function removeparticles_fields(PC, ikeep)
    nkeep = length(ikeep)
    Tp = Vector{Float64}(undef, nkeep)
    Fxxp = similar(Tp)
    Fzzp = similar(Tp)
    Fxzp = similar(Tp)
    Fzxp = similar(Tp)

    for (j, i) in enumerate(ikeep)
        @inbounds Tp[j] = PC.T[i]
        @inbounds Fxxp[j] = PC.Fxx[i]
        @inbounds Fzzp[j] = PC.Fzz[i]
        @inbounds Fxzp[j] = PC.Fxz[i]
        @inbounds Fzxp[j] = PC.Fzx[i]
    end
    PC.T = Tp
    PC.Fxx = Fxxp
    PC.Fzz = Fzzp
    PC.Fxz = Fxzp
    PC.Fzx = Fzxp
    PC
end

function addreject(T, F, gr, θThermal, rThermal, IC, particle_info, particle_weights, particle_fields; min_num_particles = 6)

    e2n, e2n_p1 = gr.e2n, gr.e2n_p1

    nel  = size(e2n_p1,2)
    vertices = [Point2D{Polar}(rand(),rand()) for _ in 1:3] # element vertices
    t0 = [particle_info[i].t for i in axes(particle_info,1)]

    max_ppel, min_ppel = 15, min_num_particles

    # --- Find number of particles per element
    reps, irem = getreps(t0, nel, max_ppel)
    # np = length(particle_info)
    # ikeep = setdiff(1:np, irem)
    # if !isempty(irem)
    #     particle_info = removeparticles(particle_info, ikeep)
    #     particle_weights = removeparticles(particle_weights, ikeep)
    #     particle_fields = removeparticles_fields(particle_fields, ikeep)
    # end

    # -- Create particles to be added
    iadd = findall(x->x<min_ppel, reps) # elements to be injected
    nxel = @. max(min_ppel - reps,0) # number of particles to inject x element
    nadd = sum(nxel)  # total number of new particles

    println(nadd, " particles added \n")

    # -- Allocatestuff
    vertices = [Point2D{Polar}(0.0, 0.0) for _ in 1:3] # element vertices

    if !isempty(iadd)

        @inbounds for iel in iadd
            # -- Parent element
            iparent = div(iel, 4, RoundUp)
            # Vertices of previous parent element
            getvertices!(vertices, θThermal, rThermal, iel)

            # # -- Find neighbouring particles
            # nb = neighbours0[iel]
            # iₙ = ismember(t0, nb)
            # particlesₙ = [particle_info[i].CPolar for i in iₙ]
            # Tₙ = particle_fields.T[iₙ]
            # Fxxₙ = particle_fields.Fxx[iₙ]
            # Fzzₙ = particle_fields.Fzz[iₙ]
            # Fxzₙ = particle_fields.Fxz[iₙ]
            # Fzxₙ = particle_fields.Fzx[iₙ]

            # -- Add particle stuff
            nrand = nxel[iel]
            for  _ in 1:nrand
                # -- Particle coordinates
                particle = randomintriangle(vertices)

                # -- Barycentric coordinates and check if point is inside triangle
                bary_coords, = barycentric_coordinates2(true, vertices, particle)

                # --- Weights (for particle → node)
                node_weights = distance(particle,vertices)
                nw = @MVector [node_weights[1],node_weights[2],node_weights[3]]

                # --- Weights (for particle → integration point)
                ip_weights = distance(particle, IC[iparent,:])
                iw = @MVector [ip_weights[1],ip_weights[2],ip_weights[3],
                               ip_weights[4],ip_weights[5],ip_weights[6]]

                # #  -- Interpolate fields to new particles
                # ω = distance(particle, particlesₙ).^-1
                # sumω = sum(ω)
                # # Deformation gradient
                # Fxxp = mydot(ω,Fxxₙ)/sumω
                # Fzzp = mydot(ω,Fzzₙ)/sumω
                # Fzxp = mydot(ω,Fxzₙ)/sumω
                # Fxzp = mydot(ω,Fzxₙ)/sumω
                # # Temperature
                # Tp = mydot(ω,Tₙ)/sumω

                # Deformation gradient
                idx = view(e2n,1:6,iparent)
                wF = @. 1/iw^2
                Fxxp = ip2particle([F[i][1,1] for i in idx], wF)
                Fzzp = ip2particle([F[i][2,2] for i in idx], wF)
                Fzxp = ip2particle([F[i][2,1] for i in idx], wF)
                Fxzp = ip2particle([F[i][1,2] for i in idx], wF)

                # Temperature
                idx = view(e2n_p1,:,iel)
                Tp = weightednode2particle(view(T, idx), fw.(bary_coords))

                # -- Add to Particle Structures
                push!(particle_info, PINFO(
                                           particle,
                                           polar2cartesian(particle),
                                           Point2D{Cartesian}(0.0,0.0),
                                           iel,
                                           iparent,
                                           zeros(3)
                                           )
                     )

                push!(particle_weights,PWEIGHTS(
                                                bary_coords,
                                                iw,
                                                nw
                                                )
                      )

                push!(particle_fields.T, Tp)
                push!(particle_fields.ΔT_subgrid, 0.0)
                push!(particle_fields.Fxx, Fxxp)
                push!(particle_fields.Fzz, Fzzp)
                push!(particle_fields.Fxz, Fxzp)
                push!(particle_fields.Fzx, Fzxp)

            end
        end

    end

    return particle_info, particle_weights,particle_fields

end

@inbounds function randomintriangle(vp)
    p, q = rand(), rand()
    if (p+q)>1
        p = 1-p
        q = 1-q
    end
    x  = @muladd vp[1].x + (vp[2].x-vp[1].x)*p + (vp[3].x-vp[1].x)*q
    z  = @muladd vp[1].z + (vp[2].z-vp[1].z)*p + (vp[3].z-vp[1].z)*q

    Point2D{Polar}(x,z)
end

function particles_generator(θ3, r3, IntC, e2n_p1; number_of_particles = 12)

    color           = zeros(3)
    vertices        = [Point2D{Polar}(rand(),rand()) for _ in 1:3] # element vertices
    particle_info   = PINFO[]
    particle_weights= PWEIGHTS[]
    nrand           = number_of_particles
    node_weights    = zeros(3)
    ip_weights      = similar(node_weights)
    xp = zeros(3)

    for iel in axes(e2n_p1,2)
        getvertices!(vertices, θ3, r3, iel)
        # check if element is crossing π and correct if needed
        vertices, check = check_crossingπ(vertices,xp)

        for _ in 1:nrand

            # generate radom particle inside triangle
            particle_polar  = randomintriangle(vertices)
            # check if particle was rotated
            if check == true
                if particle_polar.x > 2pi
                    particle_polar.x -= 2pi
                end
            end
            t_parent = ceil(Int,iel/4)

            # add to Particle Structures
            push!(particle_info, PINFO(
                                       particle_polar,
                                       polar2cartesian(particle_polar),
                                       Point2D{Cartesian}(0.0,0.0),
                                       iel,
                                       ceil(Int,iel/4),
                                       color
                                       )
            )

            # Barycentric coordinates and check if point is inside triangle
            bary_coords,    = barycentric_coordinates2(true, vertices, particle_polar)
            # Weights (for particle → node)
            node_weights   .= distance(particle_polar, vertices)
            nw              = @MVector [node_weights[1],node_weights[2],node_weights[3]]
            #  Weights (for particle → integration point)
            ip_weights      = distance(particle_polar, IntC[t_parent,:])
            iw              = @MVector [ip_weights[1],ip_weights[2],ip_weights[3],
                                        ip_weights[4],ip_weights[5],ip_weights[6]]

            push!(particle_weights,PWEIGHTS(
                                            bary_coords,
                                            iw,
                                            nw
                                            )
                  )
        end

    end

    return particle_info, particle_weights, init_pvars(length(particle_info))
end

@inline function bilinear_weigth(vc, pc, N)
    if N == 1 # index shuffle switcher
        x1,x2,x3 = vc[3].x, vc[1].x, vc[2].x
        z1,z2,z3 = vc[3].z, vc[1].z, vc[2].z
    elseif N == 2
        x1,x2,x3 = vc[2].x, vc[3].x, vc[1].x
        z1,z2,z3 = vc[2].z, vc[3].z, vc[1].z
    end
    x,z = pc.x, pc.z
    N =
        ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
end

function anycustom(a,b)
    out = false
    for i ∈ axes(b,1)
        if @inbounds b[i] == a
            out = true
            break
        end
    end
    return out
end

function findlostparticles(t_lost,particle_info,θThermal, rThermal)
    idx = unique(t_lost)
    idx = idx[2:end]
    vertices = [Point2D{Polar}(rand(),rand()) for _=1:3] # element vertices

    for i in idx
        intri = false
        while intri == false
            # -- i-th particle
            particle    = particle_info[i].CPolar
            # -- Check if it's in the same element as in previous tstep
            previous_el = particle_info[i].t
            # -- Vertices of previous parent element
            getvertices!(vertices, θThermal, rThermal, previous_el)
            # -- Calculate barycentric coordinates and check if point is inside triangle
            _, intri = barycentric_coordinates2(true, vertices, particle)
        end
    end
end

function check_parent(vertices,θThermal, rThermal, xp, particle, iel)
    # -- Vertices of previous parent element
    getvertices!(vertices, θThermal, rThermal, iel)
    # -- check if element is crossing π and correct if needed
    check_crossingπ!(vertices, particle)
    # -- Calculate barycentric coordinates and check if point is inside triangle
    bary_coords, intri = barycentric_coordinates2(true, vertices, particle)
    return bary_coords, intri
end

function ismember(a,b)
    v = fill(0,length(b))
    n = 0
    for i in axes(b,1)
        for j in axes(a,1)
            @inbounds if a[j] == b[i]
                n += 1
                v[n] = i
                break
            end
        end
    end
    view(v,1:n)
end

function intriangle(vc::NamedTuple, pc)
    # Coordinates of triangle's vertices
    x1,x2,x3 = vc.x[1], vc.x[2], vc.x[3]
    z1,z2,z3 = vc.z[1], vc.z[2], vc.z[3]
    # Coordinates of query point
    x, z = pc.x, pc.z

    intri = false
    bc1 = @muladd ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    if 0 <= bc1 <= 1
        bc2 = @muladd  ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
            ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
        if 0 <= bc2 <= 1
            bc3 = 1 - bc1 - bc2
            if 0 <= bc3 <= 1
                intri = true
            end
        end
    end
    intri
end

function intriangle(vc::Vector{Point2D{T}}, pc) where {T}
    # Coordinates of triangle's vertices
    x1,x2,x3 = vc[1].x, vc[2].x, vc[3].x
    z1,z2,z3 = vc[1].z, vc[2].z, vc[3].z
    # Coordinates of query point
    x, z = pc.x, pc.z

    intri = false
    bc1 = @muladd  ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
        ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    if 0 <= bc1 <= 1
        bc2 = @muladd ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
            ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
        if 0 <= bc2 <= 1
            bc3 = 1 - bc1 - bc2
            if 0 <= bc3 <= 1
                intri = true
            end
        end
    end
    intri
end

function barycentric(vc::Vector{Point2D{T}}, pc) where {T}
    # Coordinates of triangle's vertices
    x1, x2, x3 = vc[1].x, vc[2].x, vc[3].x
    z1, z2, z3 = vc[1].z, vc[2].z, vc[3].z

    bc1 = @muladd   ((z2 - z3) * (pc.x - x3) + (x3 - x2) * (pc.z - z3)) /
          ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    bc2 = @muladd  ((z3 - z1) * (pc.x - x3) + (x1 - x3) * (pc.z - z3)) /
          ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    bc3 = 1 - bc1 - bc2

    return bc1, bc2, bc3
end

function barycentric(vc::NamedTuple, pc)
    # Coordinates of triangle's vertices
    x1, x2, x3 = vc.x[1], vc.x[2], vc.x[3]
    z1, z2, z3 = vc.z[1], vc.z[2], vc.z[3]

    bc1 = @muladd ((z2 - z3) * (pc.x - x3) + (x3 - x2) * (pc.z - z3)) /
          ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    bc2 = @muladd ((z3 - z1) * (pc.x - x3) + (x1 - x3) * (pc.z - z3)) /
          ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    bc3 = 1 - bc1 - bc2

    return bc1, bc2, bc3
end

function tsearch_parallel(particle_info, particle_weights, θThermal, rThermal, coordinates, gr, IntC)
    np = length(particle_info) # num of particles
    found = fill(false, np) # true if particle found in same element as in previous t-step
    vertices = [[Point2D{Polar}(0.0, 0.0) for _ = 1:3] for _ in 1:Threads.nthreads()] # element vertices
    macro_neighbour = gr.neighbours
    θ, r = coordinates.θ, coordinates.r

    # This greedy algorithm is not enterely balanced, it is slighlty better to break
    # it into the three following separate balanced threaded loops
    @batch for ipart in 1:np
        # check if particle is in the same triangle as in previous time step i.e. found = true
        isinparent!(particle_info, θThermal, rThermal, vertices, found, ipart)
    end
    # this is the guy responsible of breaking the balanced load
    @batch for ipart in 1:np
        # if found[ipart] == false
            @inbounds found[ipart] && continue
            # check neighbour that contains new position of the particle
            isinneighbours!(particle_info, θThermal, rThermal, θ, r, vertices, ipart, found,  macro_neighbour)
        # end
    end

    @batch for ipart in 1:np
        # fill weights information
        fillparticle!(particle_weights, particle_info, vertices, IntC, θThermal, rThermal, ipart)
    end

    return particle_info, particle_weights, found
end

getvertices(θ, r, t) = (x = @views(θ[:, t]), z = @views(r[:, t]))

getvertices_e2n(gr::Grid, t::NTuple) =
    (x = (gr.θ[t[1]], gr.θ[t[2]], gr.θ[t[3]]), z = (gr.r[t[1]], gr.r[t[2]], gr.r[t[3]]))

@inbounds function isinparent!(particle_info, θThermal, rThermal, vertices, found, ipart)
    nt = Threads.threadid()
    # -- Vertices of previous parent element
    getvertices!(vertices[nt], θThermal, rThermal, particle_info[ipart].t)
    # -- Calculate barycentric coordinates and check if point is inside triangle
    found[ipart] = intriangle(vertices[nt], particle_info[ipart].CPolar)
end

@inbounds function isinneighbours!(particle_info, θThermal, rThermal, θ, r, vertices, ipart, found, macro_neighbour)
    nt = Threads.threadid()

    # check if it's still within same parent element (not included in neighbour list)
    iparent = particle_info[ipart].t_parent
    getvertices!(vertices[nt],  θ, r, iparent)
    if intriangle(vertices[nt], particle_info[ipart].CPolar)

        els_p1 = (4*(iparent-1)+1):(4*(iparent-1)+1)+3 # indices of child triangles

        for p in els_p1# indices of child triangles
            # Vertices of previous parent element
            getvertices!(vertices[nt],  θThermal, rThermal, p)
            # Calculate barycentric coordinates and check if point is inside triangle
            intri = intriangle(vertices[nt], particle_info[ipart].CPolar)
            if intri == true
                particle_info[ipart].t = p
                found[ipart] = true
                break
            end
        end

    else
        for mn in macro_neighbour[iparent]
            getvertices!(vertices[nt],  θ, r, mn)
            intri = intriangle(vertices[nt], particle_info[ipart].CPolar)

            if intri == true

                els_p1 = (4*(mn-1)+1):(4*(mn-1)+1)+3 # indices of child triangles

                for p in els_p1# indices of child triangles
                    # Vertices of previous parent element
                    getvertices!(vertices[nt],  θThermal, rThermal, p)
                    # Calculate barycentric coordinates and check if point is inside triangle
                    intri = intriangle(vertices[nt], particle_info[ipart].CPolar)
                    if intri == true
                        particle_info[ipart].t = p
                        found[ipart] = true
                        break
                    end
                end

            end
        end
    end

end

@inbounds function fillparticle!(particle_weights, particle_info, vertices, IntC, θThermal, rThermal, ipart)
    nt = Threads.threadid()
    particle_info[ipart].t_parent = ceil(Int,particle_info[ipart].t/4)
    pc = particle_info[ipart].CPolar
    # -- Check if element is crossing π and correct if needed
    getvertices!(vertices[nt], θThermal, rThermal, particle_info[ipart].t)
    # -- barycentric weights
    particle_weights[ipart].barycentric .= barycentric(vertices[nt], pc)
    # -- Nodal weights (for particle 2 node)
    particle_weights[ipart].node_weights .= distance(pc, vertices[nt])
    # particle_weights[ipart].ip_weights .= distance(pc,  view(IntC,particle_info[ipart].t_parent,:))
    @inbounds for ii = 1:6
        # -- Barycentric weights (for particle 2 ip)
        particle_weights[ipart].ip_weights[ii] =
            distance(pc, IntC[particle_info[ipart].t_parent,ii])
    end

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
include("../src/benchmark.jl")
include("../src/save.jl")

images, sim_params = main();

beam_energy_map, transmit_time_map = WaveSim.beam_energy_map_and_transmit_time_map(images, sim_params);

saveall(images, beam_energy_map, transmit_time_map, sim_params, ".")


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    Primitive{T<:Real}

`T` is the number type used to represent the primitive,  e.g., `Float64`.
Primitives are the basic elements which can be stored in bounding volume hierarchies and include surfaces and CSG objects

**Must** implement the following:
```julia
boundingbox(a::Primitive{T})::BoundingBox{T}
centroid(a::Primitive{T})::SVector{3,T}
```
"""
abstract type Primitive{T<:Real} end

"""
    Surface{T<:Real}

`T` is the number type used to represent the surface, e.g., `Float64`.
Basic `Surface`s are _not_ valid CSG objects, they function only in a stand-alone capacity.

**Must** implement the following:
```julia
surfaceintersection(surface::Surface{T}, ray::AbstractRay{T,3}) -> Union{EmptyInterval{T},Interval{T}}
normal(surface::Surface{T}) -> SVector{3,T}
interface(surface::Surface{T}) -> OpticalInterface{T}
makemesh(surface::Surface{T}) -> TriangleMesh{T}
```

In a conventional ray tracer the surface intersection function would only return the first surface the ray intersects. Because our ray tracer does CSG operations the surface intersection function intersects the ray with all leaf surfaces which are part of the CSG tree.

Each leaf surface returns one or more 1D intervals along the ray. These intervals contain the part of the ray which is inside the surface. The intervals computed at the leaves are propagated upward through the CSG tree and the CSG operations of union, intersection, and difference are applied to generate new intervals which are themselves propagated upward.

The result is a union of 1D intervals, which may be disjoint, a single interval, or empty. The union of intervals represents the parts of the ray which are inside the CSG object.

Inside is well defined for halfspaces such as cylinders and spheres which divide space into two parts, but not for Bezier or NURBS patches which generally do not enclose a volume.  For surfaces which are not halfspaces the notion of inside is defined locally by computing the angle between the incoming ray and the normal of the surface at the point of intersection. All surfaces must be defined so that the normal points to the outside of the surface.

A negative dot product between the incoming ray and the normal indicates the ray is coming from the outside of the surface and heading toward the inside. A positive dot product indicates the ray is coming from the inside of the surface and heading toward the outside.

Intervals are defined along the ray which is being intersected with the surface, so they are one dimensional. For example, assume we have a ray with origin o on the outside of a plane and an intersection with the plane at point int = o + td where t is a scalar and d is the unit direction of the ray. The inside interval will be (Intersection(t),Infinity). This interval begins at the intersection point on the plane and continues to positive infinity. The Intersection struct stores both the parametric value t and the 3D point of intersection to make various operations more efficient. But the interval operations only depend on the parametric value t.

If the origin o is on the inside of the plane then the inside interval will be (RayOrigin,Intersection(t)). Only the part of the ray from the ray origin to the intersection point is inside the plane.

It is the programmer's responsibility to return Interval results from surfaceintersection that maintain these properties.

The following must be impemented only if the surface is being used as a detector
```julia
uv(surface::Surface{T}, p::SVector{3,T}) -> SVector{2,T}
uvtopix(surface::Surface{T}, uv::SVector{2,T}, imsize::Tuple{Int,Int}) -> Tuple{Int,Int}
onsurface(surface::Surface{T}, p::SVector{3,T}) -> Bool
```
"""
abstract type Surface{T<:Real} <: Primitive{T} end
export Surface, surfaceintersection, normal, interface

"""
    ParametricSurface{T,N} <: Surface{T}

`T` is the number type used to represent the surface, e.g., `Float64`.
`N` is the dimension of the space the surface is embedded in.
`ParametricSurface`s are valid CSG objects, in some cases (where analytic intersection isn't possible) they must be wrapped in an [`AcceleratedParametricSurface`](@ref) for use.

**Must** implement the following:
```julia
uv(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> SVector{2,T}
uvrange(surface::ParametricSurface{T,N}) -> Tuple{Tuple{T,T},Tuple{T,T}}
point(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
partials(surface::ParametricSurface{T,N}, u::T, v::T) -> Tuple{SVector{N,T}, SVector{N,T}}
normal(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
inside(surface::ParametricSurface{T,N}, p: :SVector{N,T}) -> Bool
onsurface(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> Bool
surfaceintersection(surface::ParametricSurface{T,N}, AbstractRay::Ray{T,N}) -> Union{EmptyInterval{T},Interval{T},DisjointUnion{T}}
interface(surface::ParametricSurface{T,N}) -> OpticalInterface{T}
```
"""
abstract type ParametricSurface{S<:Real,N} <: Surface{S} end
export ParametricSurface, point, partials, uvrange, inside, onsurface, uv

# all subclasses must implement one of these at least...
"""
    point(surf::ParametricSurface{T}, u::T, v::T) -> SVector{3,T}
    point(surf::ParametricSurface{T}, uv::SVector{2,T}) -> SVector{3,T}

Returns the 3D point on `surf` at the given uv coordinate.
"""
point(s::ParametricSurface{T}, u::T, v::T) where {T<:Real} = point(s, SVector{2,T}(u, v))
point(s::ParametricSurface{T}, uv::SVector{2,T}) where {T<:Real} = point(s, uv[1], uv[2])
"""
    normal(surf::ParametricSurface{T}, u::T, v::T) -> SVector{3,T}
    normal(surf::ParametricSurface{T}, uv::SVector{2,T}) -> SVector{3,T}

Returns the normal to `surf` at the given uv coordinate.
"""
normal(s::ParametricSurface{T}, u::T, v::T) where {T<:Real} = normal(s, SVector{2,T}(u, v))
normal(s::ParametricSurface{T}, uv::SVector{2,T}) where {T<:Real} = normal(s, uv[1], uv[2])
"""
    partials(surf::ParametricSurface{T}, u::T, v::T) -> (SVector{3,T}, SVector{3,T})
    partials(surf::ParametricSurface{T}, uv::SVector{2,T}) -> (SVector{3,T}, SVector{3,T})

Returns a tuple of the 3D partial derivatives of `surf` with respect to u and v at the given uv coordinate.
"""
partials(s::ParametricSurface{T}, u::T, v::T) where {T<:Real} = partials(s, SVector{2,T}(u, v))
partials(s::ParametricSurface{T}, uv::SVector{2,T}) where {T<:Real} = partials(s, uv[1], uv[2])
"""
    uv(surf::ParametricSurface{T}, p::SVector{3,T}) -> SVector{2,T}
    uv(surf::ParametricSurface{T}, x::T, y::T, z::T) -> SVector{2,T}

Returns the uv coordinate on `surf` of a point, `p`, in 3D space.
If `onsurface(surf, p)` is false then the behavior is undefined, it may return an inorrect uv, an invalid uv, NaN or crash.
"""
uv(s::ParametricSurface{T,3}, x::T, y::T, z::T) where {T<:Real} = uv(s, SVector{3,T}(x, y, z))
uv(s::ParametricSurface{T,3}, p::SVector{3,T}) where {T<:Real} = uv(s, p[1], p[2], p[3])
"""
    inside(surf::ParametricSurface{T}, p::SVector{3,T}) -> Bool
    inside(surf::ParametricSurface{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _inside_ `surf`.
"""
inside(s::ParametricSurface{T,3}, x::T, y::T, z::T) where {T<:Real} = inside(s, SVector{3,T}(x, y, z))
inside(s::ParametricSurface{T,3}, p::SVector{3,T}) where {T<:Real} = inside(s, p[1], p[2], p[3])
"""
    onsurface(surf::ParametricSurface{T}, p::SVector{3,T}) -> Bool
    onsurface(surf::ParametricSurface{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _on_ `surf`.
"""
onsurface(s::ParametricSurface{T,3}, x::T, y::T, z::T) where {T<:Real} = onsurface(s, SVector{3,T}(x, y, z))
onsurface(s::ParametricSurface{T,3}, p::SVector{3,T}) where {T<:Real} = onsurface(s, p[1], p[2], p[3])
"""
    uvrange(s::ParametricSurface)
    uvrange(::Type{S}) where {S<:ParametricSurface}

Returns a tuple of the form: `((umin, umax), (vmin, vmax))` specifying the limits of the parameterisation for this surface type.
Also implemented for some `Surface`s which are not `ParametricSurface`s (e.g. `Rectangle`).
"""
uvrange(::S) where {S<:ParametricSurface} = uvrange(S)

"""
    samplesurface(surf::ParametricSurface{T,N}, samplefunction::Function, numsamples::Int = 30)

Sample a parametric surface on an even `numsamples`×`numsamples` grid in UV space with provided function
"""
function samplesurface(surf::ParametricSurface{T,N}, samplefunction::Function, numsamples::Int = 30) where {T<:Real,N}
    urange, vrange = uvrange(surf)
    rt = typeof(samplefunction(surf, urange[1], vrange[1]))
    samples = Vector{rt}(undef, (numsamples + 1)^2)
    ustep = (urange[2] - urange[1]) / numsamples
    vstep = (vrange[2] - vrange[1]) / numsamples
    for ui in 0:numsamples
        for vi in 0:numsamples
            u = urange[1] + ui * ustep
            v = vrange[1] + vi * vstep
            samples[ui * (numsamples + 1) + vi + 1] = samplefunction(surf, u, v)
        end
    end
    return samples
end
export samplesurface

"""
    triangulate(surf::ParametricSurface{S,N}, quads_per_row::Int, extensionu::Bool = false, extensionv::Bool = false, radialu::Bool = false, radialv::Bool = false)

Create an array of triangles representing the parametric surface where vertices are sampled on an even grid in UV space.
The surface can be extended by 1% in u and v separately, and specifying either u or v as being radial - i.e. detemining the radius on the surface e.g. rho for zernike - will result in that dimension being sampled using sqwrt so that area of triangles is uniform. The extension will also only apply to the maximum in this case.
"""
function triangulate(surf::ParametricSurface{T,N}, subdivisons::Int, extensionu::Bool = false, extensionv::Bool = false, radialu::Bool = false, radialv::Bool = false) where {T,N}
    triangles = newintrianglepool!(T)
    (umin, umax), (vmin, vmax) = uvrange(surf)

    if extensionu
        if !radialu
            umin -= TRIANGULATION_EXTENSION
        end
        umax += TRIANGULATION_EXTENSION
    end

    if extensionv
        if !radialv
            vmin -= TRIANGULATION_EXTENSION
        end
        vmax += TRIANGULATION_EXTENSION
    end

    # if we are using this for intersection then we need to expand the triangles very slightly to
    # avoid misses due to floating point precision
    expansion = extensionu || extensionv ? TRIANGULATION_EXPANSION : zero(T)

    for ui in 0:(subdivisons - 1)
        for vj in 0:(subdivisons - 1)
            # this can be inefficient as we evaluate each point twice
            # in practice we usually use extension in which case no point is the same anyway
            u1 = T(ui / subdivisons)
            u2 = T((ui + 1) / subdivisons)
            v1 = T(vj / subdivisons)
            v2 = T((vj + 1) / subdivisons)

            if radialu
                u1 = sqrt(max(u1, zero(T)))
                u2 = sqrt(max(u2, zero(T)))
            end
            if radialv
                v1 = sqrt(max(v1, zero(T)))
                v2 = sqrt(max(v2, zero(T)))
            end

            u1 = (umax - umin) * u1 + umin - expansion
            u2 = (umax - umin) * u2 + umin + expansion
            v1 = (vmax - vmin) * v1 + vmin - expansion
            v2 = (vmax - vmin) * v2 + vmin + expansion

            p1 = point(surf, u1, v1)
            p2 = point(surf, u2, v1)
            p3 = point(surf, u2, v2)
            p4 = point(surf, u1, v2)

            if validtri(p1, p2, p3)
                push!(triangles, Triangle(p1, p2, p3, SVector{2,T}(u1, v1), SVector{2,T}(u2, v1), SVector{2,T}(u2, v2)))
            end
            if validtri(p1, p3, p4)
                push!(triangles, Triangle(p1, p3, p4, SVector{2,T}(u1, v1), SVector{2,T}(u2, v2), SVector{2,T}(u1, v2)))
            end
        end
    end

    return triangles
end
export triangulate

"""
    makemesh(object, subdivisions::Int = 30) -> TriangleMesh

Creates a [`TriangleMesh`](@ref) from an object, either a [`ParametricSurface`](@ref), [`CSGTree`](@ref) or certain surfaces (e.g. `Circle`, `Rectangle`).
This is used for visualization purposes only.
"""
function makemesh(surface::ParametricSurface{S,N}, subdivisions::Int = 30)::TriangleMesh{S} where {S,N}
    m = TriangleMesh(triangulate(surface, subdivisions, false))
    emptytrianglepool!(S)
    return m
end
export makemesh


text:
"""
Save points by row in file .txt.
"""
function save_points_txt(filename::String,V::Points)
	io = open(filename,"w")

	for i in 1:size(V,2)
		for j in 1:size(V,1)

			if j == size(V,1)
				write(io, "$(V[j,i])")
			else
				write(io, "$(V[j,i]) ")
			end

		end
		write(io, "\n")
	end

	close(io)
end

"""
Save cells by row in file .txt.
"""
function save_cells_txt(filename::String,EV::Cells)
	io = open(filename,"w")
	dim_cell = length(EV[1])

	for i in 1:length(EV)
		for j in 1:dim_cell

			if j == dim_cell
				write(io, "$(EV[i][j])")
			else
				write(io, "$(EV[i][j]) ")
			end

		end
		write(io, "\n")
	end

	close(io)
end

"""
Save point cloud by row: "x y z r g b".
"""
function save_points_rgbs_txt(filename::String, PC::PointCloud)
	io = open(filename,"w")
	RGB = convert(Array{Int32,2},floor.(PC.rgbs*255))
	for i in 1:PC.n_points
		V = PC.coordinates
		for j in 1:PC.dimension
			write(io, "$(V[j,i]) ")
		end
		for j in 1:3
			if j == 3
				write(io, "$(RGB[j,i])")
			else
				write(io, "$(RGB[j,i]) ")
			end
		end
		write(io, "\n")
	end

	close(io)
end

#
# """
# Save extrema points of segment by row.
# """
# function save_3D_lines_txt(filename::String, lines::Array{Line,1}, affine_matrix::Matrix)
# 	io = open(filename,"w")
# 	for line in lines
# 		V,_ = Common.DrawLines(line)
# 		V1 = vcat(V,(zeros(size(V,2)))')
# 		V3D = Common.apply_matrix(affine_matrix,V1)
# 		write(io, "$(V3D[1,1]) $(V3D[2,1]) $(V3D[3,1]) $(V3D[1,2]) $(V3D[2,2]) $(V3D[3,2])\n")
# 	end
#
# 	close(io)
# end
#
# """
# Save extrema points of segment by row.
# """
# function save_2D_lines_txt(filename::String, lines::Array{Hyperplane,1})
# 	io = open(filename,"w")
# 	for line in lines
# 		V,_ = Common.DrawLines(line)
# 		write(io, "$(V[1,1]) $(V[2,1]) $(V[1,2]) $(V[2,2])\n")
# 	end
#
# 	close(io)
# end

#
# """
# """
# function save_hyperplane(filename::String, hyperplane::Hyperplane)
# 	V = hyperplane.inliers.coordinates
# 	dir = hyperplane.direction
# 	cen = hyperplane.centroid
# 	dim = length(dir)
#
# 	io = open(filename,"w")
#
# 	for j in 1:dim
# 		if j == dim
# 			write(io, "$(dir[j])")
# 		else
# 			write(io, "$(dir[j]) ")
# 		end
# 	end
#
# 	write(io, "\n")
#
# 	for j in 1:dim
# 		if j == dim
# 			write(io, "$(cen[j])")
# 		else
# 			write(io, "$(cen[j]) ")
# 		end
# 	end
#
# 	write(io, "\n")
#
# 	for i in 1:size(V,2)
# 		for j in 1:size(V,1)
#
# 			if j == size(V,1)
# 				write(io, "$(V[j,i])")
# 			else
# 				write(io, "$(V[j,i]) ")
# 			end
#
# 		end
# 		write(io, "\n")
# 	end
#
# 	close(io)
# end
#


text:
@testset "CVIs" begin
    @info "CVI Testing"

    # Set the approximation CVI tolerance for all comparisons
    tolerance = 1e-1

    # Construct the cvis
    cvis = [
        WB(),
        GD43(),
        GD53(),
        cSIL(),
        CH(),
        XB(),
        DB(),
        PS(),
        rCIP(),
    ]
    n_cvis = length(cvis)

    # Grab all the data paths for testing
    data_paths = readdir("../data", join=true)

    # Initialize the results data containers
    data, labels, n_samples = Dict(), Dict(), Dict()

    # Sanitize all the data
    p = 0.1
    @info "Subsampling data at rate: $p"
    for data_path in data_paths
        # Load the data, get a subset, and relabel in order
        local_data, local_labels = get_cvi_data(data_path)
        local_data, local_labels = get_bernoulli_subset(local_data, local_labels, p)
        local_labels = relabel_cvi_data(local_labels)

        # Store the sanitized data
        data[data_path] = local_data
        labels[data_path] = local_labels
        n_samples[data_path] = length(local_labels)
    end

    # Incremental
    @info "------- CVI Incremental -------"
    cvi_i = Dict()
    for data_path in data_paths
        @info "Data: $data_path"
        cvi_i[data_path] = deepcopy(cvis)
        for cvi in cvi_i[data_path]
            # @info "ICVI: $(typeof(cvi))"
            for ix = 1:n_samples[data_path]
                sample = data[data_path][:, ix]
                label = labels[data_path][ix]
                # param_inc!(cvi, data[:, ix], labels[ix])
                param_inc!(cvi, sample, label)
                evaluate!(cvi)
            end
        end
    end

    # Batch
    @info "------- CVI Batch -------"
    cvi_b = Dict()
    for data_path in data_paths
        @info "Data: $data_path"
        cvi_b[data_path] = deepcopy(cvis)
        for cvi in cvi_b[data_path]
            # @info "CVI: $(typeof(cvi))"
            param_batch!(cvi, data[data_path], labels[data_path])
            evaluate!(cvi)
        end
    end

    # Incremental porcelain
    @info "------- ICVI Porcelain -------"
    cvi_ip = Dict()
    # cvs_ip = Dict()
    for data_path in data_paths
        @info "Data: $data_path"
        cvi_ip[data_path] = deepcopy(cvis)
        # cvs_ip[data_path] = zeros(n_samples[data_path], n_cvis)
        # for cx = 1:n_cvis
        for cvi in cvi_ip[data_path]
            for ix = 1:n_samples[data_path]
                sample = data[data_path][:, ix]
                label = labels[data_path][ix]
                # _ = get_icvi!(cvi_ip[data_path][cx], sample, label)
                _ = get_icvi!(cvi, sample, label)
                # cvs_ip[data_path][ix, cx] = cv
            end
        end
    end

    # Batch porcelain
    @info "------- CVI Porcelain -------"
    cvi_bp = Dict()
    for data_path in data_paths
        @info "Data: $data_path"
        cvi_bp[data_path] = deepcopy(cvis)
        # cvs_b = zeros(n_cvis)
        # for cx = 1:n_cvis
        for cvi in cvi_bp[data_path]
            # cvs_b[cx] = get_cvi!(cvi_bp[cx], data, labels)
            # _ = get_cvi!(cvi_bp[data_path][cx], data[data_path], labels[data_path])
            _ = get_cvi!(cvi, data[data_path], labels[data_path])
        end
    end

    # Test that all permutations are equivalent
    for data_path in data_paths
        for cx = 1:n_cvis
            # I to B
            @test isapprox(cvi_i[data_path][cx].criterion_value,
                cvi_b[data_path][cx].criterion_value,
                atol=tolerance
            )

            # IP to BP
            @test isapprox(cvi_ip[data_path][cx].criterion_value,
                cvi_bp[data_path][cx].criterion_value,
                atol=tolerance
            )

            # I to IP
            @test isapprox(cvi_i[data_path][cx].criterion_value,
                cvi_ip[data_path][cx].criterion_value,
                atol=tolerance
            )

            # B to BP
            @test isapprox(cvi_b[data_path][cx].criterion_value,
                cvi_bp[data_path][cx].criterion_value,
                atol=tolerance
            )
        end
    end
end


text:
# # SPLINT package example

# `SPLINT` is a QCDNUM add-on for converting results computed on the
# evolution grid to cubic splines. This is convenient for
# integrating/differentiating these results.

using QCDNUM
using Printf

# ### Example evolution grid

# First, we need a QCDNUM evolution grid to work with.
# For this, we use the QCDNUM example program described in the `example` notebook.

xmin = Float64.([1.0e-4])
iwt = Int32.([1])
ng = 1
nxin = 100
iosp = 3
nx = 10

qq = Float64.([2e0, 1e4])
wt = Float64.([1e0, 1e0])
nq = 1
nqin = 60
ngq = 2
itype = 1

as0 = 0.364
r20 = 2.0

q2c = 3.0
q2b = 25.0
q0 = 2.0
iqt = 999

def = Float64.([0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,
      0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,
      0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
      0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
      0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]);

nfin = 0
x = 1.0e-3
q = 1.0e3
pdf = Array{Float64}(undef, 13)
qmz2 = 8315.25;

function func(ipdf, x)::Float64
    i = ipdf[]
    xb = x[]
    adbar = 0.1939875
    f = 0
    if (i == 0)
        ag = 1.7
        f = ag * xb^-0.1 * (1.0-xb)^5.0
    end
    if (i == 1)
        ad = 3.064320
        f = ad * xb^0.8 * (1.0-xb)^4.0
    end
    if (i == 2)
        au = 5.107200
        f = au * xb^0.8 * (1.0-xb)^3.0
    end
    if (i == 3)
        f = 0.0
    end
    if (i == 4)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0
    end
    if (i == 5)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
    end
    if (i == 6)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i == 7)
        f = 0.0
    end
    if (i == 8)
        f = 0.0
    end
    if (i == 9)
        f = 0.0
    end
    if (i == 10)
        f = 0.0
    end
    if (i == 11)
        f = 0.0
    end
    if (i == 12)
        f = 0.0
    end
    return f
end

func_c = @cfunction(func, Float64, (Ref{Int32}, Ref{Float64}))

# Let's set up and run this example evolution.

QCDNUM.qcinit(-6, " ")
nx = QCDNUM.gxmake(xmin, iwt, ng, nxin, iosp)
nq = QCDNUM.gqmake(qq, wt, ngq, nqin)
nw = QCDNUM.fillwt(itype)
QCDNUM.setord(3)
QCDNUM.setalf(as0, r20)
iqc = QCDNUM.iqfrmq(q2c)
iqb = QCDNUM.iqfrmq(q2b)
iqt = QCDNUM.iqfrmq(1e11)
QCDNUM.setcbt(0, iqc, iqb, 0)
iq0 = QCDNUM.iqfrmq(q0)
eps = QCDNUM.evolfg(itype, func_c, def, iq0)
pdf = QCDNUM.allfxq(itype, x, q, 0, 1)
asmz, a, b = QCDNUM.asfunc(qmz2)
csea = 2*pdf[3];

@printf("x, q, CharmSea = %0.4e, %0.4e, %0.4e\n", x, q, csea)
@printf("as(mz2) = %0.4e", asmz)

# ### SPLINT
#
# For the SPLINT example, we need a PDF stored and the
# corresponding `iset` and `ipdf`. We will use this to
# input the function to spline.

iset = itype
ipdf = 1; # component of PDF, -6 to 6

QCDNUM.fvalij(iset, ipdf, 1, 1, 1) # check we get values

# Now that we have run the example program,
# we can load SPLINT, and write `iset` and `ipdf` to
# the user space.

QCDNUM.ssp_spinit(100)
QCDNUM.ssp_uwrite(1, Float64(iset))
QCDNUM.ssp_uwrite(2, Float64(ipdf))

# Make spline object

iasp = QCDNUM.isp_s2make(5, 5)

# Define function to read from QCDNUM into spline

function func(ix, iq, first)::Float64

    ix = ix[] # deref ptr
    iq = iq[]

    iset = Int32(QCDNUM.dsp_uread(1))
    ipdf = Int32(QCDNUM.dsp_uread(2))

    return QCDNUM.fvalij(iset, ipdf, ix, iq, 1)
end

fun = @cfunction(func, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))

# Fill the spline and set no kinematic limit
QCDNUM.ssp_s2fill(iasp, fun, 0.0)

# We can use some helper functions to query the spline properties...

QCDNUM.isp_splinetype(iasp)

nu, u1, u2, nv, v1, v2, n = QCDNUM.ssp_splims(iasp)

# ... or copy the nodes locally and print a summary.

xarray = QCDNUM.ssp_unodes(iasp, nu, nu);
qqarray = QCDNUM.ssp_vnodes(iasp, nv, nv);

QCDNUM.ssp_nprint(iasp)

# There are routines to evaluate the function and
# its integral at desired x and qq values/ranges.

x = 0.1
q = 100.0
QCDNUM.dsp_funs2(iasp, x, q, 1)

x1 = 0.01
x2 = 0.1
q1 = 10.0
q2 = 100.0
rs = 370.0
np = 4
QCDNUM.dsp_ints2(iasp, x1, x2, q1, q2, rs, np)

# It is also possible to set user nodes, if the
# automatically chosen ones are not satisfactory.

xarr = Float64.([1e-2, 5e-2, 1e-1, 0.5])
qarr = Float64.([10, 1e2, 1e3, 5e3])
nx = length(xarr)
nq = length(qarr)

iasp = QCDNUM.isp_s2user(xarr, nx, qarr, nq)

# Let's check the update nodes:

QCDNUM.ssp_nprint(iasp)