text:
@testset "Testing MLK.OpenBound" begin
    for T in (FloatingPointTypes..., IntegerTypes...)
        if T <: Integer
            @test_throws ErrorException MLK.OpenBound(one(T))
        else
            if T <: AbstractFloat
                @test_throws ErrorException MLK.OpenBound(convert(T,  NaN))
                @test_throws ErrorException MLK.OpenBound(convert(T,  Inf))
                @test_throws ErrorException MLK.OpenBound(convert(T, -Inf))
            end

            for x in (1,2)
                B = MLK.OpenBound(convert(T,x))

                @test B.value  == x
                @test eltype(B) == T
            end

            for U in FloatingPointTypes
                B = MLK.OpenBound(one(T))

                B_u = convert(MLK.OpenBound{U}, B)
                @test B_u.value == one(U)
                @test eltype(B_u) == U
            end

            B = MLK.OpenBound(zero(T))

            @test checkvalue(-one(T), B) == true
            @test checkvalue(zero(T), B) == false
            @test checkvalue(one(T),  B) == false

            @test checkvalue(B, -one(T)) == false
            @test checkvalue(B, zero(T)) == false
            @test checkvalue(B, one(T))  == true

            @test string(B) == string("OpenBound(", zero(T), ")")
            @test show(devnull, B) == nothing
        end
    end
end

@testset "Testing MLK.ClosedBound" begin
    for T in (FloatingPointTypes..., IntegerTypes...)
        if T <: AbstractFloat
            @test_throws ErrorException MLK.ClosedBound(convert(T,  NaN))
            @test_throws ErrorException MLK.ClosedBound(convert(T,  Inf))
            @test_throws ErrorException MLK.ClosedBound(convert(T, -Inf))
        end

        for x in (1,2)
            B = MLK.ClosedBound(convert(T,x))

            @test B.value  == x
            @test eltype(B) == T
        end

        for U in (FloatingPointTypes..., IntegerTypes...)
            B = MLK.ClosedBound(one(T))

            B_u = convert(MLK.ClosedBound{U}, B)
            @test B_u.value == one(U)
            @test eltype(B_u) == U
        end

        B = MLK.ClosedBound(zero(T))

        @test checkvalue(-one(T), B) == true
        @test checkvalue(zero(T), B) == true
        @test checkvalue(one(T),  B) == false

        @test checkvalue(B, -one(T)) == false
        @test checkvalue(B, zero(T)) == true
        @test checkvalue(B, one(T))  == true

        @test string(B) == string("ClosedBound(", zero(T), ")")
        @test show(devnull, B) == nothing
    end
end

@testset "Testing MLK.NullBound" begin
    for T in (FloatingPointTypes..., IntegerTypes...)
        @test eltype(NullBound(T)) == T

        for U in (FloatingPointTypes..., IntegerTypes...)
            B = MLK.NullBound(T)

            B_u = convert(MLK.NullBound{U}, B)
            @test eltype(B_u) == U
        end

        B = MLK.NullBound(T)

        @test checkvalue(-one(T), B) == true
        @test checkvalue(zero(T), B) == true
        @test checkvalue(one(T),  B) == true

        @test checkvalue(B, -one(T)) == true
        @test checkvalue(B, zero(T)) == true
        @test checkvalue(B, one(T))  == true

        @test string(B) == string("NullBound(", T, ")")
        @test show(devnull, B) == nothing
    end
end

@testset "Testing MLK.Interval" begin
    for T in FloatingPointTypes
        @test_throws ErrorException MLK.Interval(MLK.ClosedBound(one(T)), MLK.OpenBound(one(T)))
        @test_throws ErrorException MLK.Interval(MLK.OpenBound(one(T)),   MLK.ClosedBound(one(T)))
        @test_throws ErrorException MLK.Interval(MLK.OpenBound(one(T)),   MLK.OpenBound(one(T)))

        a = MLK.ClosedBound(one(T))
        b = MLK.ClosedBound(one(T))
        I = MLK.Interval(a,b)
        @test I.a == a
        @test I.b == b
        @test eltype(I) == T

        I = MLK.interval(nothing, nothing)
        @test I.a == NullBound(Float64)
        @test I.b == NullBound(Float64)
        @test eltype(I) == Float64

        for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
            for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
                I = MLK.Interval(a, b)
                @test I.a == a
                @test I.b == b
                @test eltype(I) == T

                if typeof(a) <: NullBound
                    if typeof(b) <: NullBound
                        @test I == MLK.interval(T)
                        @test string(I) == string("interval(", T, ")")
                    else
                        @test I == MLK.interval(nothing, b)
                        @test string(I) == string("interval(nothing,", string(b), ")")
                    end
                else
                    if typeof(b) <: NullBound
                        @test I == MLK.interval(a,nothing)
                        @test string(I) == string("interval(", string(a), ",nothing)")
                    else
                        @test I == MLK.interval(a, b)
                        @test string(I) == string("interval(", string(a), ",", string(b), ")")
                    end
                end

                @test MLK.checkvalue(I, convert(T,-2)) == (typeof(a) <: NullBound ? true  : false)
                @test MLK.checkvalue(I, -one(T))       == (typeof(a) <: OpenBound ? false : true)
                @test MLK.checkvalue(I, zero(T))       == true
                @test MLK.checkvalue(I,  one(T))       == (typeof(b) <: OpenBound ? false : true)
                @test MLK.checkvalue(I, convert(T,2))  == (typeof(b) <: NullBound ? true  : false)
            end
        end

        a = convert(T,7.6)
        b = convert(T,23)
        c = convert(T,13)

        @test show(devnull, MLK.interval(MLK.ClosedBound(one(T)), MLK.ClosedBound(one(T)))) == nothing
    end
end

@testset "Testing MLK.interval" begin
    @test typeof(MLK.interval(nothing, nothing)) == MLK.Interval{Float64,NullBound{Float64},NullBound{Float64}}
    for T in FloatingPointTypes
        for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
            for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
                null_a = typeof(a) <: NullBound
                null_b = typeof(b) <: NullBound

                I = interval(null_a ? nothing : a, null_b ? nothing : b)

                if null_a && null_b
                    @test typeof(I) == MLK.Interval{Float64,NullBound{Float64},NullBound{Float64}}
                else
                    @test I.a == a
                    @test I.b == b
                    @test eltype(T) == T
                end
            end
        end
    end
end

@testset "Testing MLK.checkvalue" begin
    for T in FloatingPointTypes
        l = convert(T,-9)
        a = convert(T,-3)
        b = convert(T,3)
        u = convert(T,9)
        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
            T_a = typeof(A)
            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
                T_b = typeof(B)
                I = HP.interval(A,B)

                for x in range(l, stop = a, length = 30)
                    if T_a <: NullBound || T_a <: ClosedBound && x == a
                        @test MLK.checkvalue(I,x) == true
                    else
                        @test MLK.checkvalue(I,x) == false
                    end
                end

                for x in range(a, stop = b, length = 30)
                    if x == a
                        @test MLK.checkvalue(I,x) == T_a <: OpenBound ? false : true
                    elseif x == b
                        @test MLK.checkvalue(I,x) == T_b <: OpenBound ? false : true
                    else
                        @test MLK.checkvalue(I,x) == true
                    end
                end

                for x in range(b, stop = u, length = 30)
                    if T_b <: NullBound || T_b <: ClosedBound && x == b
                        @test MLK.checkvalue(I,x) == true
                    else
                        @test MLK.checkvalue(I,x) == false
                    end
                end
            end
        end

    end
end

@testset "Testing MLK.lowerboundtheta" begin
    for T in FloatingPointTypes
        for A in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
            T_a = typeof(A)
            for B in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
                T_b = typeof(B)
                I = MLK.Interval(A, B)

                if T_a <: NullBound
                    if T_b <: NullBound
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    elseif T_b <: ClosedBound
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    else
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    end
                elseif T_a <: ClosedBound
                    if T_b <: NullBound
                        @test MLK.lowerboundtheta(I) == I.a.value
                    elseif T_b <: ClosedBound
                        @test MLK.lowerboundtheta(I) == I.a.value
                    else
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    end
                else
                    if T_b <: NullBound
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    elseif T_b <: ClosedBound
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    else
                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
                    end
                end
            end
        end
    end
end

@testset "Testing MLK.upperboundtheta" begin
    for T in FloatingPointTypes
        a = convert(T,-5)
        b = convert(T,5)
        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
            T_a = typeof(A)
            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
                T_b = typeof(B)
                I = MLK.Interval(A, B)

                if T_a <: NullBound
                    if T_b <: NullBound
                        @test MLK.upperboundtheta(I) == convert(T,Inf)
                    elseif T_b <: ClosedBound
                        @test MLK.upperboundtheta(I) == b
                    else
                        @test MLK.upperboundtheta(I) == convert(T,Inf)
                    end
                elseif T_a <: ClosedBound
                    if T_b <: NullBound
                        @test MLK.upperboundtheta(I) == convert(T,Inf)
                    elseif T_b <: ClosedBound
                        @test MLK.upperboundtheta(I) == b
                    else
                        @test MLK.upperboundtheta(I) == log(b - a)
                    end
                else
                    if T_b <: NullBound
                        @test MLK.upperboundtheta(I) == convert(T,Inf)
                    elseif T_b <: ClosedBound
                        @test MLK.upperboundtheta(I) == log(b - a)
                    else
                        @test MLK.upperboundtheta(I) == convert(T,Inf)
                    end
                end
            end
        end
    end
end

@testset "Testing HP.theta" begin
    for T in FloatingPointTypes
        l = convert(T,-9)
        a = convert(T,-3)
        b = convert(T,3)
        u = convert(T,9)
        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
            T_a = typeof(A)
            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
                T_b = typeof(B)
                I = MLK.Interval(A, B)

                for x in range(l, stop = u, length = 60)
                    if T_a <: NullBound
                        if T_b <: NullBound
                            @test HP.theta(I,x) == x
                        elseif T_b <: ClosedBound
                            if x <= b
                                @test HP.theta(I,x) == x
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        else
                            if x < b
                                @test HP.theta(I,x) == log(b-x)
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        end
                    elseif T_a <: ClosedBound
                        if T_b <: NullBound
                            if a <= x
                                @test HP.theta(I,x) == x
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        elseif T_b <: ClosedBound
                            if a <= x <= b
                                @test HP.theta(I,x) == x
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        else
                            if a <= x < b
                                @test HP.theta(I,x) == log(b-x)
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        end
                    else
                        if T_b <: NullBound
                            if a < x
                                @test HP.theta(I,x) == log(x-a)
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        elseif T_b <: ClosedBound
                            if a < x <= b
                                @test HP.theta(I,x) == log(x-a)
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        else
                            if a < x < b
                                @test HP.theta(I,x) == log(x-a) - log(b-x)
                            else
                                @test_throws DomainError HP.theta(I,x)
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "Testing HP.eta" begin
    for T in FloatingPointTypes
        a = -one(T)
        b = one(T)
        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
                I = HP.interval(A,B)
                l_theta = max(HP.lowerboundtheta(I),convert(T,-10))
                u_theta = min(HP.upperboundtheta(I),convert(T,10))

                for x in range(convert(T, -10), stop = convert(T, 10), length = 201)
                    if l_theta <= x <= u_theta
                        v = HP.eta(I,x)
                        @test isapprox(HP.theta(I,v), x)
                    else
                        @test_throws DomainError HP.eta(I,x)
                    end
                end
            end
        end
    end
end

@testset "Testing MLK.HyperParameter" begin
    for T in (FloatingPointTypes..., IntegerTypes...)
        I = interval(nothing, ClosedBound(zero(T)))

        @test_throws ErrorException HyperParameter(one(T), I)

        P = HyperParameter(one(T), interval(nothing, ClosedBound(one(T))))
        @test getindex(P.value) == one(T)
        @test eltype(P) == T

        @test getvalue(P) == one(T)

        @test checkvalue(P, convert(T,2)) == false
        @test checkvalue(P, zero(T)) == true

        setvalue!(P, zero(T))
        @test getindex(P.value) == zero(T)
        @test getvalue(P) == zero(T)

        if T in FloatingPointTypes
            for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
                for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
                    I = HP.Interval(a, b)
                    P = HyperParameter(zero(T),I)

                    @test isapprox(HP.gettheta(P), HP.theta(I, zero(T)))

                    HP.settheta!(P, HP.theta(I, convert(T,0.5)))
                    @test isapprox(HP.gettheta(P), HP.theta(I, convert(T,0.5)))
                    @test isapprox(HP.getvalue(P), convert(T,0.5))
                end
            end
        end

        P1 = HyperParameter(zero(T), interval(T))
        P2 = HyperParameter(one(T),  interval(T))
        for op in (isless, ==, +, -, *, /)
            @test op(P1, one(T)) == op(getvalue(P1), one(T))
            @test op(one(T), P1) == op(one(T), getvalue(P1))
            @test op(P1, P2)     == op(getvalue(P1), getvalue(P2))
        end

        @test show(devnull, P) == nothing
    end
end


text:
module GitHub

##########
# import #
##########

import HttpCommon,
       HttpServer,
       JSON,
       MbedTLS,
       Requests

#############
# Utilities #
#############

# include -------

include("utils/GitHubType.jl")
include("utils/auth.jl")
include("utils/requests.jl")

# export -------

export # auth.jl
       authenticate

export # requests.jl
       rate_limit

##################################
# Owners (organizations + users) #
##################################

# include -------

include("owners/owners.jl")

# export -------

export # owners.jl
       Owner,
       owner,
       orgs,
       followers,
       following,
       repos

################
# Repositories #
################

# include -------

include("repositories/repositories.jl")
include("repositories/contents.jl")
include("repositories/commits.jl")
include("repositories/branches.jl")
include("repositories/statuses.jl")

# export -------

export # repositories.jl
       Repo,
       repo,
       create_fork,
       forks,
       contributors,
       collaborators,
       iscollaborator,
       add_collaborator,
       remove_collaborator,
       stats

export # contents.jl
       Content,
       file,
       directory,
       create_file,
       update_file,
       delete_file,
       readme,
       permalink

export # commits.jl
       Commit,
       commit,
       commits

export # branches.jl
       Branch,
       branch,
       branches

export # statuses.jl
       Status,
       create_status,
       statuses,
       status

##########
# Issues #
##########

# include -------

include("issues/pull_requests.jl")
include("issues/issues.jl")
include("issues/comments.jl")

# export -------

export # pull_requests.jl
       PullRequest,
       pull_requests,
       pull_request

export # issues.jl
       Issue,
       issue,
       issues,
       create_issue,
       edit_issue

export # comments.jl
       Comment,
       comment,
       comments,
       create_comment,
       edit_comment,
       delete_comment

############
# Activity #
############

# include -------

include("activity/events.jl")
include("activity/activity.jl")

# export -------

export # activity.jl
       star,
       unstar,
       stargazers,
       starred,
       watchers,
       watched,
       watch,
       unwatch

export # events/events.jl
       WebhookEvent

export # events/listeners.jl
       EventListener,
       CommentListener

end # module GitHub

ext: jl


text:
export Tree
export findall, setchildren!, topdown, topdown_while, toxml, tosexpr, todict

mutable struct Tree
    name
    children::Vector{Tree}
    parent
end

function Tree(name, children::Vector{Tree})
    t = Tree(name, children, nothing)
    for i = length(children):-1:1
        c = children[i]
        delete!(c)
        c.parent = t
    end
    t
end
Tree(name, children::Tree...) = Tree(name, [children...])
Tree(name) = Tree(name, Tree[])

isroot(tree::Tree) = isa(tree.parent, Void)

Base.getindex(tree::Tree, key::Int) = tree.children[key]
function Base.getindex(tree::Tree, key::String)
    for c in tree.children
        c.name == key && return c
    end
    throw("$key is not found.")
end

function Base.setindex!(tree::Tree, value::Tree, key::Int)
    tree.children[key] = value
    value.parent = tree
end
Base.isempty(tree::Tree) = isempty(tree.children)
Base.length(tree::Tree) = length(tree.children)
Base.size(tree::Tree, i::Int) = size(tree.children, i)
Base.endof(tree::Tree) = endof(tree.children)

function Base.push!(tree::Tree, children::Tree...)
    push!(tree.children, children...)
    for i = length(children):-1:1
        c = children[i]
        delete!(c)
        c.parent = tree
    end
end
Base.append!(tree::Tree, children) = push!(tree, children...)
function Base.prepend!(tree::Tree, children::Vector)
    prepend!(tree.children, children)
    for i = length(children):-1:1
        c = children[i]
        delete!(c)
        c.parent = tree
    end
end
function Base.insert!(tree::Tree, i::Int, child::Tree)
    insert!(tree.children, i, child)
    child.parent = tree
end

function Base.deleteat!(tree::Tree, i::Int)
    tree[i].parent = nothing
    deleteat!(tree.children, i)
end
function Base.delete!(tree::Tree)
    tree.parent == nothing && return
    i = findfirst(x -> x == tree, tree.parent.children)
    deleteat!(tree.parent, i)
end
function Base.delete!(tree::Tree, key::String)
    inds = find(c -> c.name == key, tree.children)
    for i = length(inds):-1:1
        deleteat!(tree, inds[i])
    end
end
function removeat!(tree::Tree, i::Int)
    children = Tree[]
    tree[i].parent = nothing
    for k = 1:length(tree)
        if k == i
            append!(children, tree[k].children)
        else
            push!(children, tree[k])
        end
    end
    setchildren!(tree, children)
end
function Base.empty!(tree::Tree)
    foreach(c -> c.parent = nothing, tree.children)
    empty!(tree.children)
end

function findall(f::Function, tree::Tree)
    nodes = Tree[]
    function traverse(node::Tree)
        f(node) && push!(nodes,node)
        for c in node.children
            traverse(c)
        end
    end
    traverse(tree)
    nodes
end
findall(tree::Tree, name::String) = findall(x -> x.name == name, tree)

function replace!(oldt::Tree, newt::Tree)
    p = oldt.parent
    i = findfirst(x -> x == oldt, p.children)
    deleteat!(p, i)
    insert!(p, i, newt)
end

function setchildren!(tree::Tree, children::Vector{Tree})
    foreach(c -> c.parent = nothing, tree.children)
    foreach(c -> c.parent = tree, children)
    tree.children = children
end
setchildren!(tree::Tree, children::Tree...) = setchildren!(tree, [children...])

function topdown(f, tree::Tree)
    f(tree)
    for c in tree.children
        topdown(f, c)
    end
end

function topdown_while(f::Function, tree::Tree)
    cond = f(tree)
    @assert isa(cond,Bool)
    cond || return
    for c in tree.children
        topdown_while(f, c)
    end
end

function bottomup(f::Function, tree::Tree)
    for c in tree.children
        bottomup(f, c)
    end
    f(tree)
end

function tosexpr(tree::Tree)
    strs = String["(", tree.name]
    for c in tree.children
        push!(strs, tosexpr(c))
    end
    push!(strs, ")")
    join(strs)
end

function toxml(tree::Tree)
    function escape(str::String)
        str = replace(str, "&", "&")
        str = replace(str, ">", ">")
        str = replace(str, "<", "<")
        str = replace(str, "'", "'")
        str = replace(str, "\"", """)
        str
    end

    @assert !isempty(tree)
    strs = String[]
    tree.parent == nothing && push!(strs,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>")

    push!(strs, "<$(tree.name)>")
    for c in tree.children
        if isempty(c)
            push!(strs, escape(c.name))
        else
            !isempty(strs) && strs[end][end] == '>' && push!(strs,"\n")
            push!(strs, toxml(c))
        end
    end
    closetag = split(tree.name, " ")[1]
    !isempty(strs) && !isempty(strs[end]) && strs[end][end] == '>' && push!(strs,"\n")
    push!(strs, "</$closetag>")
    join(strs)
end

function todict(tree::Tree)
    n = count(isempty, tree.children)
    @assert n == 0 || n == length(tree)
    if n == 0
        Dict(tree.name => map(todict, tree.children))
    else
        tree.name
    end
end

function Base.string(tree::Tree, delim="")
    strs = String[]
    topdown(tree) do node
        isempty(node) && push!(strs,node.name)
    end
    join(strs, delim)
end

function Base.parse(::Type{Tree}, sexpr::String)
    sexpr = Vector{Char}(sexpr)
    function f(i::Int)
        chars = Char[]
        children = Tree[]
        while i <= length(sexpr)
            c = sexpr[i]
            if c == '('
                child, i = f(i+1)
                push!(children, child)
            elseif c == ')'
                tree = Tree(join(chars), children)
                return tree, i+1
            else
                c != ' ' && c != '\n' && push!(chars,sexpr[i])
                i += 1
            end
        end
        throw("Invalid S-expression.")
    end
    i = findfirst(c -> c == '(', sexpr)
    i == 0 && throw("Invalid S-expression.")
    tree, _ = f(i+1)
    tree
end

#=
struct TopdownIterator
    tree::Tree
end
topdown = TopdownIterator

Base.start(iter::TopdownIterator) = Tree[iter.tree]
Base.done(iter::TopdownIterator, state::Vector{Tree}) = isempty(state)
function Base.next(iter::TopdownIterator, state::Vector{Tree})
    tree = pop!(state)
    for i = length(tree):-1:1
        push!(state, tree[i])
    end
    (tree, state)
end
=#


text:
module FractionWrapper

using JavaCall

const AbstractFormat = @jimport org.hipparchus.fraction.AbstractFormat
const AttributedCharacterIterator = @jimport java.text.AttributedCharacterIterator
const BigDecimal = @jimport java.math.BigDecimal
const BigFraction = @jimport org.hipparchus.fraction.BigFraction
const BigFractionField = @jimport org.hipparchus.fraction.BigFractionField
const BigFractionFormat = @jimport org.hipparchus.fraction.BigFractionFormat
const BigInteger = @jimport java.math.BigInteger
const Class = @jimport java.lang.Class
const Currency = @jimport java.util.Currency
const FieldPosition = @jimport java.text.FieldPosition
const Format = @jimport java.text.Format
const Fraction = @jimport org.hipparchus.fraction.Fraction
const FractionField = @jimport org.hipparchus.fraction.FractionField
const FractionFormat = @jimport org.hipparchus.fraction.FractionFormat
const JString = @jimport java.lang.String
const Locale = @jimport java.util.Locale
const Number = @jimport java.lang.Number
const NumberFormat = @jimport java.text.NumberFormat
const NumberFormat_Style = @jimport java.text.NumberFormat$Style
const Object = @jimport java.lang.Object
const ParsePosition = @jimport java.text.ParsePosition
const ProperBigFractionFormat = @jimport org.hipparchus.fraction.ProperBigFractionFormat
const ProperFractionFormat = @jimport org.hipparchus.fraction.ProperFractionFormat
const RoundingMode = @jimport java.math.RoundingMode
const StringBuffer = @jimport java.lang.StringBuffer

export AbstractFormat
export BigFraction
export BigFractionField
export BigFractionFormat
export Fraction
export FractionField
export FractionFormat
export ProperBigFractionFormat
export ProperFractionFormat
export abs
export add
export big_decimal_value
export byte_value
export clone
export compare_to
export divide
export double_value
export equals
export float_value
export format
export format_big_fraction
export format_fraction
export format_to_character_iterator
export get_available_locales
export get_class
export get_compact_number_instance
export get_currency
export get_currency_instance
export get_denominator
export get_denominator_as_int
export get_denominator_as_long
export get_denominator_format
export get_field
export get_improper_instance
export get_instance
export get_integer_instance
export get_maximum_fraction_digits
export get_maximum_integer_digits
export get_minimum_fraction_digits
export get_minimum_integer_digits
export get_number_instance
export get_numerator
export get_numerator_as_int
export get_numerator_as_long
export get_numerator_format
export get_one
export get_percent_instance
export get_proper_instance
export get_reduced_fraction
export get_rounding_mode
export get_runtime_class
export get_whole_format
export get_zero
export hash_code
export int_value
export is_grouping_used
export is_parse_integer_only
export long_value
export multiply
export negate
export notify
export notify_all
export parse
export parse_object
export percentage_value
export pow
export reciprocal
export reduce
export set_currency
export set_grouping_used
export set_maximum_fraction_digits
export set_maximum_integer_digits
export set_minimum_fraction_digits
export set_minimum_integer_digits
export set_parse_integer_only
export set_rounding_mode
export short_value
export signum
export subtract
export to_string
export wait

include("abstract_format.jl")
include("big_fraction.jl")
include("big_fraction_field.jl")
include("big_fraction_format.jl")
include("fraction.jl")
include("fraction_field.jl")
include("fraction_format.jl")
include("proper_big_fraction_format.jl")
include("proper_fraction_format.jl")

end


text:
include("InvariantDistributions/InvariantDistributions.jl")

abstract type InvariantMeasure end
abstract type AbstractRectangularInvariantMeasure <: InvariantMeasure end
abstract type AbstractTriangulationInvariantMeasure <: InvariantMeasure end

abstract type InvariantMeasureEstimator end

export
    InvariantMeasure,
    InvariantMeasureEstimator,
    AbstractRectangularInvariantMeasure,
    AbstractTriangulationInvariantMeasure


text:
"""
    function cudaAdapt_(code)

Function that adapts vector declarations and math functions to CUDA.
"""
function cudaAdapt_(code::Expr)

    #Subs all math
    l = [i for i in names(Base.Math) if isdefined(CUDA,i) && i != :(^)]
    for i in l
        code = subs_(code,i,Meta.parse(string("CUDA.","$i")))
    end
    #Subs power
    code = postwalk(x->@capture(x, a_^b_) && typeof(b) != Int ? :(Float32($a)^$b) : x, code)

    #Subs vectors
    for i in [:zeros,:ones]
        code = subs_(code,i,Meta.parse(string("CUDA.","$i")))
    end
    code = subs_(code,:Array,Meta.parse(string("CUDA.CuArray")))
    code = subs_(code,INT,INTCUDA)
    code = subs_(code,FLOAT,FLOATCUDA)

    return code
end


text:
# Common XDR routines (for the XTC and TRR formats)

module XDR

export
		read_xdr_num, write_xdr_num, read_xdr_int, write_xdr_int,
		read_xdr_float, size_of_xdr_int, size_of_xdr_ints!, decode_xdr_bits!,
		decode_xdr_ints!

@inline read_xdr_num(io::IO, T::Type) = ntoh(read(io, T))

@inline write_xdr_num(io::IO, x::Number, T::Type) =
		write(io, hton(T(x)))

@inline read_xdr_int(io::IO) = read_xdr_num(io, Int32)

@inline write_xdr_int(io::IO, x::Integer) = write_xdr_num(io, x, Int32)

@inline read_xdr_float(io::IO) = read_xdr_num(io, Float32)

function size_of_xdr_int(size::T) where {T<:Integer}
	local num::T, num_of_bits::T
	num = 1
	num_of_bits = 0
	while (size >= num && num_of_bits < 32)
		num_of_bits += 1
		num <<= 1
	end
	num_of_bits
end

function size_of_xdr_ints!(bytes::AbstractVector{T},
		sizes::AbstractVector{T}) where {T<:Integer}
	local num_of_bytes::T, num_of_bits::T, tmp::T, bytecnt::T, num::T,
			last_byte::T
	num_of_bytes = 1
	num_of_bits = 0
	bytes[1] = 1
	for size in sizes
		tmp = 0
		bytecnt = 1
		while bytecnt <= num_of_bytes
			tmp = bytes[bytecnt] * size + tmp
			bytes[bytecnt] = tmp & 0xff
			tmp >>= 8
			bytecnt += 1
		end
		while tmp != 0
			bytes[bytecnt] = tmp & 0xff
			bytecnt += 1
			tmp >>= 8
		end
		num_of_bytes = bytecnt - 1
	end
	num = 1
	last_byte = bytes[num_of_bytes]
	num_of_bytes -= 1
	while last_byte >= num
		num_of_bits += 1
		num *= 2
	end
	num_of_bits + num_of_bytes * 8
end

function decode_xdr_bits!(state::AbstractVector{T}, buf::AbstractVector{UInt8},
		num_of_bits::T) where {T<:Integer}
	local mask::T, cnt::T, lastbits::T, lastbyte::T, num::T

	mask = (1 << num_of_bits) - 1
	cnt = state[1]
	lastbits = state[2]
	lastbyte = state[3]
	num = 0

	while num_of_bits >= 8
		lastbyte = (lastbyte << 8) | buf[cnt+1]
		cnt += 1
		num |= (lastbyte >> lastbits) << (num_of_bits - 8)
		num_of_bits -= 8
	end
	if num_of_bits > 0
		if lastbits < num_of_bits
			lastbits += 8
			lastbyte = (lastbyte << 8) | buf[cnt+1]
			cnt += 1
		end
		lastbits -= num_of_bits
		num |= (lastbyte >> lastbits) & ((1 << num_of_bits) - 1)
	end
	num &= mask
	state[1] = cnt
	state[2] = lastbits
	state[3] = lastbyte
	num
end

function decode_xdr_ints!(nums::AbstractVector{T}, bytes::AbstractVector{T},
		state::AbstractVector{T}, buf::AbstractVector{UInt8}, num_of_ints::T,
		num_of_bits::T, sizes::AbstractVector{T}, pos::T) where {T<:Integer}
	@inbounds begin
		local num_of_bytes::T, i::T, j::T, p::T
		bytes[2] = 0
		bytes[3] = 0
		bytes[4] = 0
		num_of_bytes = 0
		while num_of_bits > 8
			bytes[num_of_bytes+1] = decode_xdr_bits!(state, buf, T(8))
			num_of_bytes += 1
			num_of_bits -= T(8)
		end
		if num_of_bits > 0
			bytes[num_of_bytes+1] = decode_xdr_bits!(state, buf, num_of_bits)
			num_of_bytes += 1
		end
		i = num_of_ints - 1
		while i > 0
			num::Int = 0
			j = num_of_bytes - 1
			while j >= 0
				num = (num << 8) | bytes[j+1]
				p = unsafe_trunc(Int, num / sizes[i+1])
				bytes[j+1] = p
				num -= p * sizes[i+1]
				j -= 1
			end
			nums[i+pos] = num
			i -= 1
		end
		nums[pos] = bytes[1] | (bytes[2] << 8) | (bytes[3] << 16) |
				(bytes[4] << 24)
	end
	nothing
end

end # module


text:
# start the following on the command line:
# julia -p n  # starts n processes on the local machine
# julia -p 8  # starts REPL with 8 workers
workers()
# 8-element Array{Int64,1}:
#   2
#   3
#   4
#   5
#   ⋮
#   8
#   9

# iterate over workers:
for pid in workers()
	# do something with pid
end
myid() # 1
addprocs(5)
# 5-element Array{Any,1}:
#  10
#  11
#  12
#  13
#  14
nprocs() # 14
rmprocs(3) # worker with id 3 is removed

# cluster of computers:
# julia --machinefile machines driver.jl
# Run processes specified in driver.jl on hosts listed in machines

# primitive operations with processes:
r1 = remotecall(2, x -> x^2, 1000) #> RemoteRef(2,1,20)
fetch(r1) #> 1000000

remotecall_fetch(5, sin, 2pi) # -2.4492935982947064e-16

r2 = @spawnat 4 sqrt(2) # lets worker 4 calculate sqrt(2)
fetch(r2)  #> 1.4142135623730951
r = [@spawnat w sqrt(5) for w in workers()]
# or in a freshly started REPL:
# r = [@spawnat i sqrt(5) for i=1:length(workers())]
# 6-element Array{Any,1}:
#  RemoteRef(1,1,8)
#  RemoteRef(2,1,9)
#  RemoteRef(3,1,10)
#  RemoteRef(4,1,11)
#  RemoteRef(5,1,12)
#  RemoteRef(6,1,13)
r3 = @spawn sqrt(5) #> RemoteRef(5,1,26)
fetch(r3) #> 2.23606797749979

# using @everywhere to make functions available to all workers:
@everywhere w = 8
@everywhere println(myid()) # 1
        # From worker 2:  2
        # From worker 3:  3
        # From worker 4:  4
        # From worker 7:  7
        # From worker 5:  5
        # From worker 6:  6
        # From worker 8:  8
        # From worker 9:  9

@everywhere include("defs.jl")
@everywhere function fib(n)
	if (n < 2) then
	return n
	else return fib(n-1) + fib(n-2)
	end
end
require("functions")

# broadcast data to all workers:
d = "Julia"
for pid in workers()
   remotecall(pid, x -> (global d; d = x; nothing), d)
end
#= Output:
1   From worker 2:  2

    From worker 4:  4
    From worker 5:  5
    From worker 6:  6
=#

max_stars_repo_path: Module 1/Chapter08/parallel.jl

max_issues_repo_path: Module 1/Chapter08/parallel.jl

max_forks_repo_path: Module 1/Chapter08/parallel.jl


text:
module Trees
using ForceImport
@force using Base

export Tree, EmptyTree, TreeNode,
       isterminal, insert_child!,
       tree, lisp_tree_structure, parenthesis_to_brackets,
       start, next, done, eltype,
       leafs, leaf_data, data,
       iterate, eltype, IteratorSize

abstract type Tree{T} end

mutable struct EmptyTree{T} <: Tree{T}
end

mutable struct TreeNode{T} <: Tree{T}
  data    ::T
  parent  ::Tree{T}
  children::Vector{TreeNode{T}}
end

isterminal(tree::TreeNode{T}) where T = isempty(tree.children)

TreeNode(data::T, parent::Tree{T}) where T = TreeNode(data, parent, Vector{TreeNode{T}}())

# TreeNode{T}(data::T, parent::Tree{T}) = TreeNode(data, parent, Vector{TreeNode{T}}())
TreeNode(data, T=typeof(data)) = TreeNode(data, EmptyTree{T}())

data(t::TreeNode) = t.data

show(io::IO, tree::EmptyTree{T}) where T = print(io::IO, "[]")

function show(io::IO, tree::TreeNode{T}) where T
  print(io, "[", tree.data)
  for child in tree.children
    print(io, child)
  end
  print(io, "]")
end

function insert_child!(tree::TreeNode{T}, data::T) where T
  push!(tree.children, TreeNode(data, tree))
end

function insert_child!(tree::TreeNode, child::TreeNode)
    child.parent = tree
    push!(tree.children, child)
end

function tree(str::AbstractString)
  node = TreeNode("")
  str = replace(str, " "=>"")[2:end-1]

  for c in str
    if c == '['
      insert_child!(node, "")
      node = node.children[end]
    elseif c == ']'
      node = node.parent
    else
      node.data = string(node.data, c)
    end
  end
  node
end

function tree(lst::Union{Tuple, AbstractVector}, T=typeof(lst[1]))
    node = TreeNode(lst[1], T)
    if length(lst) > 1
        for x in lst[2:end]
            insert_child!(node, tree(x, T))
        end
    end
    node
end

function lisp_tree_structure(str::String)
  # str = replace(str, "  ", " ")
  node = TreeNode("")
  str = replace(str, "  ", " ")[2:end-1]
  prev = '['
  start_reading_terminal = false
  last_is_terminal = false

  for c in str
    if c == '['
      insert_child!(node, "")
      node = node.children[end]
      start_reading_terminal = false
      last_is_terminal = false
    elseif c == ']'
      if last_is_terminal
        node = node.parent.parent
      else
        node = node.parent
      end
      start_reading_terminal = false
      last_is_terminal = false
    elseif c == ' '
      start_reading_terminal = prev != '[' && prev != ']'
    else
      if start_reading_terminal
        if last_is_terminal
          node = node.parent
        end
        insert_child!(node, "")
        node = node.children[end]
        start_reading_terminal = false
        last_is_terminal = true
      end
      node.data = string(node.data, c)
    end
    prev = c
  end

  if prev == ']'
    node
  else
    node.parent
  end
end

parenthesis_to_brackets(str::String) =
  replace(replace(str, "(", "["), ")", "]")

start(tree::TreeNode{T}) where T = [tree]
next(::TreeNode{T}, list::Vector{TreeNode{T}}) where T =
  list[1], prepend!(list[2:end], list[1].children)
done(::TreeNode{T}, list::Vector{TreeNode{T}}) where T = isempty(list)
eltype(::Type{TreeNode{T}}) where T = TreeNode{T}

length(tree::TreeNode) =
    isterminal(tree) ? 1 : 1 + sum(length(c) for c in tree.children)

leafs(tree) = [n for n in tree if isterminal(n)]
leaf_data(tree) = map(x->x.data, leafs(tree))

function Base.map(f, t::Tree)
    n = TreeNode(f(t.data))
    for c in t.children
        insert_child!(n, map(f, c))
    end
    n
end

eltype(::Type{Tree}) = Pointer
IteratorSize(::Type{Tree}) = Base.SizeUnknown()
function iterate(tree::Tree, state = [tree])
    if isempty(state)
        nothing
    else
        state[1], prepend!(state[2:end], collect(values(state[1].children)))
    end
end

# t = tree("[1[2[3][4][6]][5]]")
# t2 = map(x->2parse(x), t)

# print(tree(str))
#
# str = "[S [NP phpSyntaxTree][VP [V creates][NP nice syntax trees]]]"
# print(lisp_tree_structure(str))

# str = "((SBAR (IN Because) (S (NP (NN realEDASHestate) (NNS purchases) (CC and) (NNS leases)) (VP (VBP are) (NP (NP (JJ such) (JJ major) (JJ longEDASHterm) (NNS commitments)) (SBAR (IN that) (S (NP (JJS most) (NNS companies) (CC and) (NNS individuals)) (VP (VBP make) (NP (DT these) (NNS decisions)) (SBAR (WHADVP (RB only) (WRB when)) (S (ADJP (JJ confident) (PP (IN of) (NP (JJ future) (NX (NX (JJ economic) (NN stability)) (CC and) (NX (NN growth)))))))))))))) (EDOT EDOT)))"
# println(parenthesis_to_brackets(str))
# println()
# println(lisp_tree_structure(parenthesis_to_brackets(str)))

end # module


text:
#=
# Cachesim example
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("cachesim");
@useLog("cachesimlog")

#################################################
# For cache sim
cachesim(true); # Select cachesim. Normal solve won't happen
l1 = build_cache_level(1, 32, 8, 64, "LRU");    #
l2 = build_cache_level(2, 256, 4, 64, "LRU");   # Set up the cache specs (this is my laptop)
l3 = build_cache_level(3, 12288, 16, 64, "LRU");# (level, sets, ways, line_size, policy)
build_cache([l1, l2, l3]);
#################################################

n = 10;
ord = 2;

# Set up the configuration
@domain(3)
@functionSpace(LEGENDRE, ord)

# Specify the problem
@mesh(HEXMESH, n)

@variable(u)
@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)

# Write the weak form
@coefficient(f, "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

cachesimSolve(u);

# Dump config to log
log_dump_config();
log_dump_prob();

@finalize()


text:
module GyroKinetics4d

    include("gyro_kinetics_4d/gc_solovev_iter_xpoint.jl")

end


text:
using LightGraphs

"""
    DAGHasLoops

Exception, thrown when loops are detected on DAG construction.

# Examples
```jldoctest
julia> using Dagitty

julia> g = DAG(:A => :B, :B => :A)
ERROR: DAGHasLoops()
```
"""
struct DAGHasLoops <: Exception end

"""
    LabelNotFound

Exception thrown on attempt to resolve label which is not present in the DAG.

# Examples
```jldoctest
julia> using Dagitty

julia> g = DAG(:A => :B)
DAG: {2, 1} directed simple Int64 graph with labels [:A, :B])

julia> node(g, :C)
"""
struct LabelNotFound <: Exception
    label::Symbol
end

"""
    DAGEdge

Pair of symbols representing DAG edge.
"""
DAGEdge = Pair{Symbol, Symbol}

"""
    DAG

Structure with fields:
* `graph::SimpleDiGraph` - underlying directed graph
* `labels::Vector{Symbol}` - vector of each nodes' label

Could be constructed by passing list of `Pair{Symbol, Symbol}`

# Examples
```jldoctest
julia> using Dagitty

julia> g = DAG(:A => :B)
DAG: {2, 1} directed simple Int64 graph with labels [:A, :B])

julia> g.graph
{2, 1} directed simple Int64 graph

julia> g.labels
2-element Vector{Symbol}:
 :A
 :B

julia> g = DAG(:A => :A)
ERROR: DAGHasLoops()
```
"""
struct DAG
    graph::SimpleDiGraph
    labels::Vector{Symbol}
end

DAG(edges::DAGEdge...) = DAG(edges)

function DAG(edges)
    labels = collect(unique(Iterators.flatten(edges)))
    sort!(labels)
    vert_count = length(labels)
    lab_indices = Dict(zip(labels, Base.OneTo(vert_count)))
    graph = SimpleDiGraph(vert_count)
    for (a, b) ∈ edges
        has_edge(graph, a, b) && continue
        add_edge!(graph, lab_indices[a], lab_indices[b])
    end
    simplecyclescount(graph) > 0 && throw(DAGHasLoops())
    DAG(graph, labels)
end

function Base.show(io::IO, g::DAG)
    println(io, "DAG: $(g.graph) with labels $(g.labels))")
end

export
    DAG,
    DAGEdge,
    DAGHasLoops,
    LabelNotFound


text:
module WhereTheWaterFlows

using StaticArrays, Requires

export waterflows

const I11 = CartesianIndex(1,1)
const I22 = CartesianIndex(2,2)
const NOFLOW = 5     # direction number indicating no flow.  Use a constant to better keep track.
const NOFLOWer = 10  # direction number indicating no flow as well as no flow into this cell.

"""
Direction numbers.  E.g. dirnums[1,1] will return the number
corresponding to the direction top-left.

Note, I use the conversion that the x-axis corresponds to the row of
the matrix and the y-axis the columns.  To print them in this "normal"
coordinate system use `showme`
"""
const dirnums = SMatrix{3,3}(reverse([ 7      8 9
                                       4 NOFLOW 6
                                       1      2 3]',
                                     dims=2))

"Translation from dirnums to CartesianIndex"
const cartesian = SMatrix{3,3}(reverse(permutedims([CartesianIndex(-1,1)  CartesianIndex(0,1)  CartesianIndex(1,1)
                                                    CartesianIndex(-1,0)  CartesianIndex(0,0)  CartesianIndex(1,0)
                                                    CartesianIndex(-1,-1) CartesianIndex(0,-1) CartesianIndex(1,-1)]),
                                       dims=2))

"Adjustment factor for gradient in diagonal vs non-diagonal"
diagonal_fac(J::CartesianIndex) = (J.I[1]==0 || J.I[2]==0) ?  1.0 : 1/sqrt(2)

"""
Show an array with its axes oriented such that they correspond
to x and y direction.
"""
showme(ar) = (display(reverse(ar',dims=1)); println(" "))

"Translate a CartesianIndex, giving the offset, into a direction number"
ind2dir(ind::CartesianIndex) = dirnums[ind + I22]

"Translate a D8 direction number into a CartesianIndex (i.e. a flow vector). Maps dir==10 to no-flow also."
dir2ind(dir) = dir==10 ? CartesianIndex(0,0) : cartesian[dir]

"Translate a D8 direction number into a 2D vector."
dir2vec(dir) = [dir2ind(dir).I...]

"""
Tests whether a cell `J` with flowdir `dirJ` flows into cell `I`.
"""
flowsinto(J::CartesianIndex, dirJ::Integer, I::CartesianIndex) = ind2dir(I-J) == dirJ

"Return CartesianIndices corresponding to the 8 neighbors and the point itself"
iterate_D9(I, Iend, I1=I11) = max(I1, I-I1):min(Iend, I+I1)
function iterate_D9(I, ar::AbstractMatrix)
    R = CartesianIndices(size(ar))
    I1, Iend = first(R), last(R)
    return max(I1, I-I1):min(Iend, I+I1)
end

"""
    on_outer_boundary(ar, I::CartesianIndex)

Check whether on outer boundary of ar
"""
function on_outer_boundary(ar, I::CartesianIndex)
    i,j = I.I
    iend, jend = size(ar)
    if i==1 || i==iend || j==1 || j==jend
        return true
    else
        return false
    end
end

"""
    d8dir_feature(dem, bnd_as_pits)

D8 directions of a DEM and drainage features.

Elevations with NaN map to dir==NOFLOW, i.e. just like pits.
However, they are not treated in the pits-filling function `drainpits`.

The argument `bnd_as_pits` determines whether neighboring cells
flow out of the domain or into NaN-cells (`bnd_as_pits==true`)
or not.

Return
- dir  - direction, encoded as `dirnums`
- nout - number of outflow cells of a cell (0 or 1)
- nin  - number of inflow cells of a cell (0-8)
- pits - location of pits as a `Vector{CartesianIndex{2}}` (sorted)
- flowdir_extra_output -- nothing (not used by this function)
"""
function d8dir_feature(dem, bnd_as_pits)
    # outputs
    dir = zeros(Int8, size(dem))
    nout = falses(size(dem))
    nin = zeros(Int8, size(dem))
    pits = CartesianIndex{2}[]

    R = CartesianIndices(size(dem))
    Iend = last(R)

    # get dir for all points
    for I in R
        # make pits on boundary if bnd_as_pits is set
        if bnd_as_pits && on_outer_boundary(dem,I)
            # make it a pit
            dir[I] = NOFLOW
            continue
        end

        ele = dem[I]
        delta_ele = 0.0 # keeps track of biggest elevation change
        dir_ = NOFLOW
        if isnan(ele)
            # just mark as NOFLOW
        else
            for J in iterate_D9(I, Iend)
                I==J && continue
                ele2 = dem[J]
                if isnan(ele2)
                    if bnd_as_pits
                        # flow into first found NaN-cell
                        dir_ = ind2dir(J-I)
                        break
                    else
                        # ignore NaN-Cell
                        continue
                    end
                end
                delta_ele2 = (ele2 - ele) * diagonal_fac(J-I)
                if delta_ele2 < delta_ele
                    # lower elevation found, adjust dir
                    delta_ele = delta_ele2
                    dir_ = ind2dir(J-I)
                end
            end
        end
        dir[I] = dir_
    end
    # flow features
    # (not thread-safe because of push!)
    for I in R
        for J in iterate_D9(I, Iend)
            J==I && continue
            nin[I] += flowsinto(J, dir[J], I)
        end
        if dir[I]==NOFLOW
            nout[I] = false
            if !isnan(dem[I])
                push!(pits, I)
            elseif nin[I]>0
                # mark NaN points only as pits if something is flowing into them
                push!(pits, I)
            end
        else
            nout[I] = true
        end
    end

    return dir, nout, nin, pits, dem, nothing
end

"""
    waterflows(dem, cellarea=ones(size(dem)), flowdir_fn=d8dir_feature;
               calc_streamlength=true, drain_pits=true, bnd_as_pits=false)

Does the water flow routing according the D8 algorithm.

args:
- dem -- the DEM (or hydro-potential); array
- cellarea=ones(size(dem)) -- the source per cell, defaults to 1.  If using physical units
                              then use a volumetric flux per cell, e.g. m3/s.
- flowdir_fn=d8dir_feature -- the routing function.  Defaults to the built-in `d8dir_feature`
                              function but could be customised

kwargs:
- drain_pits -- whether to route through pits (true)
- bnd_as_pits (true) -- whether the domain boundary and NaNs should be pits,
                 i.e. adjacent cells can drain into them,
                 or whether to ignore them.

TODO: if bnd_as_pits routes water along boundary edges first. This would probably
substantially reduce the number of catchments, as currently every boundary point
is a pit and thus a catchment (if bnd_as_pits==true).

Returns
- area -- upslope area
- stream-length -- length of stream to the farthest source
- dir -- flow direction at each location
- nout -- whether the point has outlflow.  I.e. nout[I]==0 --> I is a pit
- nin -- number of inflow cells
- pits -- location of pits as Vector{CartesianIndex{2}}
- c -- catchment map
- bnds -- boundaries between catchments.  The boundary to the exterior/NaNs is not in here.
- flowdir_extra_output -- extra output of the flowdir_fn, which is `nothing` for the default
"""
function waterflows(dem, cellarea=ones(size(dem)), flowdir_fn=d8dir_feature;
                    drain_pits=true, bnd_as_pits=true)
    dir, nout, nin, pits, dem4drainpits, flowdir_extra_output = flowdir_fn(dem, bnd_as_pits)
    area, slen, c = flowrouting_catchments(dir, pits, cellarea)
    bnds = make_boundaries(c, 1:length(pits))
    if drain_pits
        drainpits!(dir, nin, nout, pits, c, bnds, dem4drainpits)
        area, slen, c = flowrouting_catchments(dir, pits, cellarea)
    end
    #area[isnan.(dem)] .= NaN
    return area, slen, dir, nout, nin, pits, c, bnds, flowdir_extra_output
end

# TODO: implement this instead of waterflows
# """
#     flowrouting(dir, nin, cellarea; maxiter=max(size(dir)...)*2, calc_streamlength=true)

# Do the actual routing

# Input:
# - `dir` direction field
# - `nin` number of inputs into cell
# - `cellarea` water input per cell

# KW:
# - `maxiter=max(size(dir)...)*2` -- may need to be increased for very sinous drainage paths
# - `calc_streamlength=true` -- calculate stream length.  Will slow-down calculation somewhat

# Return
# - upstream area
# - stream-length (or something irrelevant if calc_streamlength==false)
# """

"""
    flowrouting_catchments(dir, pits, cellarea)

Recursively calculate flow-routing and catchments from
- dir - direction field
- pits - pit coordinates
- cellarea - water input

Return:
- upstream area
- stream length
- catchments Matrix{Int}.  Value==0 corresponds to NaNs in the DEM
  which are not pits (i.e. where no water flows into).
"""
function flowrouting_catchments(dir, pits, cellarea)
    c = zeros(Int, size(dir))
    area = zeros(size(cellarea))
    slen = zeros(Int, size(dir))
    np = length(pits)
    # recursively traverse the drainage tree in up-flow direction,
    # starting at all pits
    Threads.@threads for color = 1:length(pits)
        pit = pits[color]
        _flowrouting_catchments!(area, slen, c, dir, cellarea, color, pit)
    end
    return area, slen, c
end
# modifies c and area
function _flowrouting_catchments!(area, len, c, dir, cellarea, color, ij)
    c[ij] = color

    # proc upstream points
    uparea = 0.0
    slen = 0
    for IJ in iterate_D9(ij, c)
        if ij==IJ
            uparea += cellarea[ij]
            slen = max(slen, 1)
        elseif flowsinto(IJ, dir[IJ], ij)
            uparea_, slen_ = _flowrouting_catchments!(area, len, c, dir, cellarea, color, IJ)
            uparea += uparea_
            slen = max(slen, slen_+1)
        end
    end
    area[ij] = uparea
    len[ij] = slen
    return uparea, slen
end

"""
    make_boundaries(catchments, colors)

Make vectors of boundary points.  Note that points along the
edge of the domain as well as points bordering NaN-cells with no
inflow do not count as boundaries.
"""
function make_boundaries(catchments, colors)
    bnds = [CartesianIndex{2}[] for c in colors]
    for R in CartesianIndices(size(catchments))
        c = catchments[R]
        c==0 && continue # don't find boundaries for c==0 (NaNs with no inflow)
        bnd = bnds[c]
        for I in iterate_D9(R, catchments)
            co = catchments[I]
            if co!=c && co!=0
                push!(bnd, R)
                break
            end
        end
    end
    bnds
end

"""
Checks that all points in `bnds[color]` are indeed on the
boundary; removes them otherwise.

TODO: this is one of the performance bottlenecks.
"""
function _prune_boundary!(del, bnds, catchments::Matrix, color, colormap)
    empty!(del)
    # loop over all boundary cells
    @inbounds for (i,P) in enumerate(bnds[color])
        keep = false
        # check cells around it
        for PP in iterate_D9(P, catchments)
            # if one is of different color, keep it.
            co = catchments[PP]
            if co!=0
                co = colormap[co]
                if co!=color
                    keep = true
                    break
                end
            end
        end
        if !keep
            push!(del, i)
        end
    end
    deleteat!(bnds[color], del)
    return nothing
end

"""
     drainpits!(dir, nin, nout, pits, c, bnds, dem;
               maxiter=100)

Return an updated direction field which drains (interior) pits.
This is done by reversing the flow connecting the lowest point on the catchment boundary
to the pit for each catchment.

There needs to be a decision on how to treat the outer boundary and boundaries
to NaN-cells.  Possibilities:
- do not treat boundary cells as pits but do treat pits at the boundary as terminal.
- treat boundary cells as pits, i.e. flow reaching such a cell will vanish.  Again
  such cells would be terminal.  This is currently done

TODO: What to do if there are no terminal pits in the DEM?  Fill to the uppermost pit?  Or take
the lowermost as terminal? -> currently it picks a random one

Returns new dir, nin, nout, pits (sorted), c, bnds

TODO: this is the performance bottleneck.
"""
function drainpits!(dir, nin, nout, pits, c, bnds, dem)
    maxiter = 100
    Iend = CartesianIndex(size(dem))

    pits_to_keep = trues(length(pits))

    # Table which translates the old color to the new color.
    # Note that only the currently processed color might change.
    # Initialize to the id-map (0 is used for off-points)
    colormap = collect(1:length(pits))

    colormap_inv = [[i] for i=1:length(pits)] # note that the total length

    no_drainage_across_boundary = false

    # work array as otherwise much time is spent allocating
    # (if below loop is ever multi-threaded, this will need one per thread)
    del_workarray  = Int[]
    sizehint!(del_workarray, maximum(length.(bnds)))

    # iterate until all interior pits are removed (not really needed, I think)
    @inbounds for i=1:maxiter
        n_removed = 0
        for (color, P) in enumerate(pits)
            # Already removed pit, skip
            P==CartesianIndex(-1,-1) && continue
            # Don't process pits on the DEM boundary, because there water disappears.
            (P.I[1]==1 || P.I[2]==1 || P.I[1]==size(dir,1) || P.I[2]==size(dir,2)) && continue
            # Don't process pits which are NaNs, because there water disappears.
            # See also bnd_as_pits.
            isnan(dem[P]) && continue
            # If there are no more boundaries left, stop.  This should only occur when
            # bnd_as_pits==false and when there are only interior pits.
            if all((isempty(b) for b in bnds))
                no_drainage_across_boundary = true # set flag to warn later
                break
            end
            # If there are no boundaries for this point, go to next point
            isempty(bnds[color]) && continue

            # Find point on catchment boundary with minimum elevation
            # Avoid, if possible, picking a NOFLOWer point.
            minn = typemax(eltype(dem))
            Imin = CartesianIndex(-1,-1)
            for I in bnds[color]
                if dem[I] < minn && dir[I]!=NOFLOWer
                    minn = dem[I]
                    Imin = I
                end
            end
            if Imin == CartesianIndex(-1,-1) # no NOFLOWer cells found
                for I in bnds[color]
                    if dem[I] < minn
                        minn = dem[I]
                        Imin = I
                    end
                end
            end

            @assert colormap[c[Imin]]==c[P]==color # Something is amiss if the found minimum is the pit!
            # make the outflow and find the next catchment
            min_ = Inf
            target = Imin
            for J in iterate_D9(Imin, dem)
                Imin==J && continue # do not look at the point itself
                cc = c[J]
                cc = cc==0 ? 0 : colormap[cc]
                cc == color && continue # J is in Imin's catchment
                if dem[J] < min_
                    min_ = dem[J]
                    target = J
                end
            end
            if target==Imin
                # this means that point is on a NaN boundary -> don't process
                continue
            end

            # Reverse directions on path going from Imin to P
            P1 = Imin
            P2 = dir2ind(dir[P1]) + P1
            # first, do the flow across the boundary
            _flow_from_to!(Imin, target, dir, nin, nout)
            while P1!=P
                # get down-downstream point (do ahead lookup as a undisturbed dir is needed)
                P3 = dir2ind(dir[P2]) + P2
                # reverse
                _flow_from_to!(P2, P1, dir, nin, nout)

                P1 = P2
                P2 = P3
            end

            # remove from list of pits
            pits_to_keep[color] = false
            pits[color] = CartesianIndex(-1,-1)

            # update colormap and boundaries
            othercolor = colormap[c[target]]
            ## this is slow
            # Threads.@threads for i in eachindex(colormap)
            #     if color==colormap[i]
            #         colormap[i] = othercolor
            #     end
            # end
            ## thus do this instead:
            append!(colormap_inv[othercolor], colormap_inv[color])
            empty!(colormap_inv[color])
            for color in colormap_inv[othercolor]
                colormap[color] = othercolor
            end

            append!(bnds[othercolor], bnds[color])
            empty!(bnds[color])
            _prune_boundary!(del_workarray, bnds, c, othercolor, colormap)
            n_removed +=1
        end
        n_removed==0 && break
        if i==maxiter
            error("Maximum number of iterations reached in `drainpits`: $i")
        end
    end

    if no_drainage_across_boundary
        @warn """Water cannot leave this DEM!  Instead it drains into a random interior pit.
                 Consider setting `bnd_as_pits=true`."""
    end

    # make new colors consecutive
    newcolors = unique(colormap)
    d = Dict((c=>i for (i,c) in enumerate(newcolors)))
    colormap2 = [d[c] for c in colormap]

    # remove removed pits and bnds
    deleteat!(pits, .!pits_to_keep)
    deleteat!(bnds, .!pits_to_keep)

    # # update catchments -> done with `flowrouting_catchments`
    return nothing
end

"""
Update dir, nin, and nout such that flow at P1 is now
from P1 to P2.

It can potentially modify `dir, nin, nout` at three locations
- P1: dir, nout, nin
- P2: nin
  - if allow_reversion==true, then P2's dir, nout can also be modified.
- P3 (previous receiver cell of P1): nin

Note that if P1 and P2 lie in the same catchment, then P3 is also in that catchment.
"""
function _flow_from_to!(P1, P2, dir, nin, nout, allow_reversion=false)
    # already right
    ind2dir(P2-P1)==dir[P1] && return nothing

    # otherwise update
    P3 = P1 + dir2ind(dir[P1]) # previous receiver cell
    if P3!=P1
        nin[P3] -= 1
    end
    dir[P1] = ind2dir(P2-P1)
    nin[P2] += 1
    nout[P1] = true # definitely not a pit anymore

    # if flow from P2 was into P1, then make P2 a pit
    # (otherwise dir will be inconsistent)
    if (dir2ind(dir[P2]) + P2 == P1)
        if allow_reversion
            dir[P2] = ind2dir(CartesianIndex(0,0))
            nout[P2] = false
            nin[P1] -= 1
        else
            error("Flow direction in P2 would need to be reversed but not allowed (set option `allow_reversion`).")
        end
    end
    return nothing
end

## Plotting
function __init__()
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")
end

## Post processing
include("postproc.jl")

end # module


text:
"""
$(SIGNATURES)

This function is the analog of [`continuation`](@ref) when the two first points on the branch are passed (instead of a single one). Hence `x0` is the first point on the branch (with palc `s=0`) with parameter `par0` and `x1` is the second point with parameter `set(par0, lens, p1)`.
"""
function continuation(F, J, x0::Tv, par0, x1::Tv, p1::Real, lens::Lens, contParams::ContinuationPar; linearAlgo = BorderingBLS(), kwargs...) where Tv
	# Create a bordered linear solver using the newton linear solver provided by the user
	# dont modify it if the user passed its own version
	if isnothing(linearAlgo.solver)
	_linearAlgo = @set linearAlgo.solver = contParams.newtonOptions.linsolver
	else
		_linearAlgo = linearAlgo
	end
	# check the sign of ds
	dsfactor = sign(p1 - get(par0, lens))
	# create an iterable
	_contParams = @set contParams.ds = abs(contParams.ds) * dsfactor
	it = ContIterable(F, J, x0, par0, lens, _contParams, _linearAlgo; kwargs...)
	return continuation(it, x0, get(par0, lens), x1, p1)
end

function continuation(it::ContIterable, x0, p0::Real, x1, p1::Real)
	## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	# The result type of this method
	# is not known at compile time so we
	# need a function barrier to resolve it
	#############################################

	# we compute the cache for the continuation, i.e. state::ContState
	# In this call, we also compute the initial point on the branch (and its stability) and the initial tangent
	state, _ = iterateFromTwoPoints(it, x0, p0, x1, p1)

	# variable to hold the result from continuation, i.e. a branch
	contRes = ContResult(it, state)

	# perform the continuation
	return continuation!(it, state, contRes)
end

"""
$(SIGNATURES)

Automatic branch switching at branch points based on a computation of the normal form. More information is provided in [Branch switching](@ref Branch-switching-page). An example of use is provided in [2d generalized Bratu–Gelfand problem](@ref).

# Arguments
- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`...
- `br` branch result from a call to [`continuation`](@ref)
- `ind_bif` index of the bifurcation point in `br` from which you want to branch from
- `optionsCont` options for the call to [`continuation`](@ref)

# Optional arguments
- `Jᵗ` associated jacobian transpose, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`.
- `δ` used internally to compute derivatives (with finite differences) w.r.t the parameter `p`.
- `δp` used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by `optionsCont.ds`. This allows to use a step larger than `optionsCont.dsmax`.
- `ampfactor = 1` factor which alters the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
- `nev` number of eigenvalues to be computed to get the right eigenvector
- `issymmetric` whether the Jacobian is Symmetric, avoid computing the left eigenvectors in the computation of the reduced equation.
- `usedeflation = true` whether to use nonlinear deflation (see [Deflated problems](@ref Deflated-problems)) to help finding the guess on the bifurcated branch
- `kwargs` optional arguments to be passed to [`continuation`](@ref), the regular `continuation` one.

!!! tip "Advanced use"
    In the case of a very large model and use of special hardware (GPU, cluster), we suggest to discouple the computation of the reduced equation, the predictor and the bifurcated branches. Have a look at `methods(BifurcationKit.multicontinuation)` to see how to call these versions. These methods has been tested on GPU with very high memory pressure.
"""
function continuation(F, dF, d2F, d3F, br::ContResult, ind_bif::Int, optionsCont::ContinuationPar = br.contparams ;
		Jᵗ = nothing,
		δ::Real = 1e-8, δp = nothing, ampfactor::Real = 1,
		nev = optionsCont.nev, issymmetric = false,
		usedeflation::Bool = false,
		Teigvec = getvectortype(br),
		scaleζ = norm,
		verbosedeflation::Bool = false,
		maxIterDeflation::Int = min(50, 15optionsCont.newtonOptions.maxIter),
		perturb = identity,
		kwargs...)
	# The usual branch switching algorithm is described in Keller. Numerical solution of bifurcation and nonlinear eigenvalue problems. We do not use this one but compute the Lyapunov-Schmidt decomposition instead and solve the polynomial equation instead.

	verbose = get(kwargs, :verbosity, 0) > 0 ? true : false
	verbose && println("--> Considering bifurcation point:"); _show(stdout, br.specialpoint[ind_bif], ind_bif)

	if kernelDim(br, ind_bif) > 1
		return multicontinuation(F, dF, d2F, d3F, br, ind_bif, optionsCont; Jᵗ = Jᵗ, δ = δ, δp = δp, ampfactor = ampfactor, nev = nev, issymmetric = issymmetric, scaleζ = scaleζ, verbosedeflation = verbosedeflation, maxIterDeflation = maxIterDeflation, perturb = perturb, Teigvec = Teigvec, kwargs...)
	end

	@assert br.type == :Equilibrium "Error! This bifurcation type is not handled.\n Branch point from $(br.type)"
	@assert br.specialpoint[ind_bif].type == :bp "Error! This bifurcation type is not handled.\n Branch point from $(br.specialpoint[ind_bif].type)"

	# compute predictor for point on new branch
	ds = isnothing(δp) ? optionsCont.ds : δp
	Ty = typeof(ds)

	# compute the normal form of the bifurcation point
	specialpoint = computeNormalForm1d(F, dF, d2F, d3F, br, ind_bif; Jᵗ = Jᵗ, δ = δ, nev = nev, verbose = verbose, issymmetric = issymmetric, Teigvec = Teigvec, scaleζ = scaleζ)

	# compute predictor for a point on new branch
	pred = predictor(specialpoint, ds; verbose = verbose, ampfactor = Ty(ampfactor))
	if isnothing(pred); return nothing, nothing; end

	verbose && printstyled(color = :green, "\n--> Start branch switching. \n--> Bifurcation type = ", type(specialpoint), "\n----> newp = ", pred.p, ", δp = ", br.specialpoint[ind_bif].param - pred.p, "\n")

	if usedeflation
		verbose && println("\n----> Compute point on the current branch with nonlinear deflation...")
		optn = optionsCont.newtonOptions
		bifpt = br.specialpoint[ind_bif]
		# find the bifurcated branch using nonlinear deflation
		solbif, _, flag, _ = newton(F, dF, convert(Teigvec, bifpt.x), pred.x, setParam(br, pred.p), setproperties(optn; verbose = verbose = verbosedeflation); kwargs...)[1]
		copyto!(pred.x, solbif)
	end

	# perform continuation
	branch, u, τ = continuation(F, dF,
			specialpoint.x0, specialpoint.params,	# first point on the branch
			pred.x, pred.p,					# second point on the branch
			br.lens, optionsCont; kwargs...)
	return Branch(branch, specialpoint), u, τ
end

# same but for a Branch
continuation(F, dF, d2F, d3F, br::AbstractBranchResult, ind_bif::Int, optionsCont::ContinuationPar = br.contparams ; kwargs...) = continuation(F, dF, d2F, d3F, getContResult(br), ind_bif, optionsCont ; kwargs...)

"""
$(SIGNATURES)

Automatic branch switching at branch points based on a computation of the normal form. More information is provided in [Branch switching](@ref). An example of use is provided in [2d generalized Bratu–Gelfand problem](@ref).

# Arguments
- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`...
- `br` branch result from a call to [`continuation`](@ref)
- `ind_bif` index of the bifurcation point in `br` from which you want to branch from
- `optionsCont` options for the call to [`continuation`](@ref)

# Optional arguments
- `Jᵗ` associated jacobian transpose, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`.
- `δ` used internally to compute derivatives w.r.t the parameter `p`.
- `δp` used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by `optionsCont.ds`. This allows to use a step larger than `optionsCont.dsmax`.
- `ampfactor = 1` factor which alters the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
- `nev` number of eigenvalues to be computed to get the right eigenvector
- `issymmetric` whether the Jacobian is Symmetric, avoid computing the left eigenvectors in the computation of the reduced equation.
- `verbosedeflation = true` whether to display the nonlinear deflation iterations (see [Deflated problems](@ref Deflated-problems)) to help finding the guess on the bifurcated branch
- `scaleζ` norm used to normalize eigenbasis when computing the reduced equation
- `Teigvec` type of the eigenvector. Useful when `br` was loaded from a file and this information was lost
- `ζs` basis of the kernel
- `perturbGuess = identity` perturb the guess from the predictor just before the deflated-newton correction
- `kwargs` optional arguments to be passed to [`continuation`](@ref), the regular `continuation` one.

!!! tip "Advanced use"
    In the case of a very large model and use of special hardware (GPU, cluster), we suggest to discouple the computation of the reduced equation, the predictor and the bifurcated branches. Have a look at `methods(BifurcationKit.multicontinuation)` to see how to call these versions. These methods has been tested on GPU with very high memory pressure.
"""
function multicontinuation(F, dF, d2F, d3F, br::AbstractBranchResult, ind_bif::Int, optionsCont::ContinuationPar = br.contparams;
		Jᵗ = nothing,
		δ::Real = 1e-8,
		δp = nothing,
		ampfactor::Real = getvectoreltype(br)(1),
		nev::Int = optionsCont.nev,
		issymmetric::Bool = false,
		Teigvec = getvectortype(br),
		ζs = nothing,
		verbosedeflation::Bool = false,
		scaleζ = norm,
		perturbGuess = identity,
		kwargs...)

	verbose = get(kwargs, :verbosity, 0) > 0 ? true : false

	bpnf = computeNormalForm(F, dF, d2F, d3F, br, ind_bif; Jᵗ = Jᵗ, δ = δ, nev = nev, verbose = verbose, issymmetric = issymmetric, Teigvec = Teigvec, ζs = ζs, scaleζ = scaleζ)

	return multicontinuation(F, dF, br, bpnf, optionsCont; Teigvec = Teigvec, δp = δp, ampfactor = ampfactor, verbosedeflation = verbosedeflation, kwargs...)
end

# for AbstractBifurcationPoint (like Hopf, BT, ...), it must return nothing
multicontinuation(F, dF, br::AbstractBranchResult, bpnf::AbstractBifurcationPoint, optionsCont::ContinuationPar; kwargs...) = nothing

# general function for branching from Nd bifurcation points
function multicontinuation(F, dF, br::AbstractBranchResult, bpnf::NdBranchPoint, optionsCont::ContinuationPar = br.contparams; δp = nothing, ampfactor = getvectoreltype(br)(1), perturb = identity, kwargs...)

	verbose = get(kwargs, :verbosity, 0) > 0 ? true & get(kwargs, :verbosedeflation, true) : false

	# compute predictor for point on new branch
	ds = abs(isnothing(δp) ? optionsCont.ds : δp)

	# get prediction from solving the reduced equation
	rootsNFm, rootsNFp = predictor(bpnf, ds;  verbose = verbose, perturb = perturb, ampfactor = ampfactor)

	return multicontinuation(F, dF, br, bpnf, (before = rootsNFm, after = rootsNFp), optionsCont; δp = δp, kwargs...)
end

"""
$(SIGNATURES)

Function to transform predictors `solfromRE` in the normal form coordinates of `bpnf` into solutions. Note that `solfromRE = (before = Vector{vectype}, after = Vector{vectype})`.
"""
function getFirstPointsOnBranch(F, dF, br::AbstractBranchResult,
		bpnf::NdBranchPoint, solfromRE,
		optionsCont::ContinuationPar = br.contparams ;
		δp = nothing,
		Teigvec = getvectortype(br),
		verbosedeflation = false,
		maxIterDeflation = min(50, 15optionsCont.newtonOptions.maxIter),
		lsdefop = DeflatedLinearSolver(),
		perturbGuess = identity,
		kwargs...)
	# compute predictor for point on new branch
	ds = isnothing(δp) ? optionsCont.ds : δp |> abs
	dscont = abs(optionsCont.ds)

	rootsNFm = solfromRE.before
	rootsNFp = solfromRE.after

	# attempting now to convert the guesses from the normal form into true zeros of F
	optn = optionsCont.newtonOptions

	# options for newton
	cbnewton = get(kwargs, :callbackN, cbDefault)
	normn = get(kwargs, :normN, norm)

	printstyled(color = :magenta, "--> Looking for solutions after the bifurcation point...\n")
	defOpp = DeflationOperator(2, 1.0, Vector{typeof(bpnf.x0)}(), _copy(bpnf.x0); autodiff = true)
	for (ind, xsol) in pairs(rootsNFp)
		# print("\n--> attempt to converge zero #$ind")
		solbif, _, flag, _ = newton(F, dF, perturbGuess(bpnf(xsol, ds)), setParam(br, bpnf.p + ds), setproperties(optn; maxIter = maxIterDeflation, verbose = verbosedeflation), defOpp, lsdefop; callback = cbnewton, normN = normn)
		flag && push!(defOpp, solbif)
	end

	printstyled(color = :magenta, "--> Looking for solutions before the bifurcation point...\n")
	defOpm = DeflationOperator(2, 1.0, Vector{typeof(bpnf.x0)}(), _copy(bpnf.x0); autodiff = true)
	for (ind, xsol) in pairs(rootsNFm)
		# print("\n--> attempt to converge zero #$ind")
		solbif, _, flag, _ = newton(F, dF, perturbGuess(bpnf(xsol, ds)), setParam(br, bpnf.p - ds), setproperties(optn; maxIter = maxIterDeflation, verbose = verbosedeflation), defOpm, lsdefop; callback = cbnewton, normN = normn)
		flag && push!(defOpm, solbif)
	end
	printstyled(color=:magenta, "--> we find $(length(defOpm)) (resp. $(length(defOpp))) roots after (resp. before) the bifurcation point.\n")
	return (before = defOpm, after = defOpp, bpm = bpnf.p - ds, bpp = bpnf.p + ds)
end

# In this function, I keep usedeflation although it is not used to simplify the calls
function multicontinuation(F, dF, br::AbstractBranchResult,
		bpnf::NdBranchPoint, solfromRE,
		optionsCont::ContinuationPar = br.contparams ;
		δp = nothing,
		Teigvec = getvectortype(br),
		verbosedeflation = false,
		maxIterDeflation = min(50, 15optionsCont.newtonOptions.maxIter),
		lsdefop = DeflatedLinearSolver(),
		perturbGuess = identity,
		kwargs...)

	defOpm, defOpp, _, _ = getFirstPointsOnBranch(F, dF, br, bpnf, solfromRE, optionsCont; δp = δp, verbosedeflation = verbosedeflation, maxIterDeflation = maxIterDeflation, lsdefop = lsdefop, perturbGuess = perturbGuess, kwargs...)

	multicontinuation(F, dF, br,
			bpnf, defOpm, defOpp, optionsCont;
			δp = δp,
			Teigvec = Teigvec,
			verbosedeflation = verbosedeflation,
			maxIterDeflation = maxIterDeflation,
			lsdefop = lsdefop,
			kwargs...)
end

"""
$(SIGNATURES)

Automatic branch switching at branch points based on a computation of the normal form. More information is provided in [Branch switching](@ref). An example of use is provided in [2d generalized Bratu–Gelfand problem](@ref).

# Arguments
- `F, dF`: function `(x, p) -> F(x, p)` and its differential `(x, p, dx) -> d1F(x, p, dx)`
- `br` branch result from a call to [`continuation`](@ref)
- `bpnf` normal form
- `defOpm::DeflationOperator, defOpp::DeflationOperator` to specify converged points on nonn-trivial branches before/after the bifurcation points.

The rest is as the regular `multicontinuation` function.
"""
function multicontinuation(F, dF, br::AbstractBranchResult,
		bpnf::NdBranchPoint, defOpm::DeflationOperator, defOpp::DeflationOperator,
		optionsCont::ContinuationPar = br.contparams ;
		δp = nothing,
		Teigvec = getvectortype(br),
		verbosedeflation = false,
		maxIterDeflation = min(50, 15optionsCont.newtonOptions.maxIter),
		lsdefop = DeflatedLinearSolver(),
		kwargs...)

	ds = isnothing(δp) ? optionsCont.ds : δp |> abs
	dscont = abs(optionsCont.ds)
	par = bpnf.params

	# compute the different branches
	function _continue(_sol, _dp, _ds)
		# needed to reset the tangent algorithm in case fields are used
		empty!(get(kwargs, :tangentAlgo, nothing))
		println("#"^50)
		continuation(F, dF,
			bpnf.x0, par,		# first point on the branch
			_sol, bpnf.p + _dp, # second point on the branch
			br.lens, (@set optionsCont.ds = _ds); kwargs...)
		# continuation(F, dF, _sol, setParam(br, bpnf.p + _dp), br.lens, (@set optionsCont.ds = _ds); kwargs...)
	end

	branches = Branch[]
	for id in 2:length(defOpm)
		br, = _continue(defOpm[id], -ds, -dscont); push!(branches, Branch(br, bpnf))
		# br, = _continue(defOpm[id], -ds, dscont); push!(branches, Branch(br, bpnf))
	end

	for id in 2:length(defOpp)
		br, = _continue(defOpp[id], ds, dscont); push!(branches, Branch(br, bpnf))
		# br, = _continue(defOpp[id], ds, -dscont); push!(branches, Branch(br, bpnf))
	end

	return branches, (before = defOpm, after = defOpp)
end

# same but for a Branch
multicontinuation(F, dF, d2F, d3F, br::Branch, ind_bif::Int, optionsCont::ContinuationPar = br.contparams; kwargs...) = multicontinuation(F, dF, d2F, d3F, getContResult(br), ind_bif, optionsCont ; kwargs...)


text:
"""
    BacktrackingLineSearch{C, T}

`BacktrackingLineSearch` starts from an initial step length guess and gradually backtracks
until the Armijo sufficient decrease condition is met.

# Fields
- `c::C`: the constant for the Armijo condition in the open interval (0, 1)
- `scale::T`: the scale factor to shrink the step length
"""
struct BacktrackingLineSearch{C, T}
    c::C
    scale::T

    function BacktrackingLineSearch(c, scale)
        0 < c < 1 || error("`c` must be in the open interval (0, 1)")
        0 < scale < 1 || error("Scale factor must be between 0 and 1")
        return new{typeof(c), typeof(scale)}(c, scale)
    end
end

"""
    BacktrackingLineSearch(; kwargs...)

Initiate `BacktrackingLineSearch`.

# Keywords
- `c=1e-4`: the constant for the Armijo condition in the open interval (0, 1)
- `scale=0.5`: the scale factor to shrink the step length
"""
BacktrackingLineSearch(; c=1e-4, scale=0.5) = BacktrackingLineSearch(c, scale)

function (bls::BacktrackingLineSearch)(f, state, p, α₀::T) where {T}
    x = state.x
    y = state.f
    ∇y = state.∇f
    dϕ₀ = ∇y⋅p
    c = bls.c
    scale = T(bls.scale)
    while f(x + α₀*p) > y + c*α₀*dϕ₀
        α₀ *= scale
    end
    return α₀
end


text:
@callbackFunction(
    function adiabatic_bdry(intensity, vg, sx, sy, dir, band, normal)
        # debug = (fid == 1 && t > 0.29);
        result = intensity[dir + (band-1)*ndirs];

        specularity = 0.5;
        #specular_tol = 1e-5; # tolerance for selecting specular reflection direction
        reflect_tol = 1e-6; # tolerance for ignoring reflection (parallel to bdry -> no reflection)
        s = [sx[dir], sy[dir]];
        s_dot_n = s[1] * normal[1] + s[2] * normal[2];
        # if debug println("dir="*string(dir)*", normal="*string(normal)*", s_dot_n="*string(s_dot_n)*", val="*string(result)) end

        # Does the direction point in or out of the domain at this boundary?
        if s_dot_n < -1e-6 # into the interior
            sspec = s .- 2*s_dot_n .* normal; # specular reflection vector
            close_to_specular = 0.5;
            angle = 2*pi/ndirs;

            # For every direction that points out of the boundary, add some reflected intensity to this.
            specular_part = 0;
            diffuse_part = 0;
            ndiffuse = 0;
            for di=1:ndirs
                si = [sx[di], sy[di]];
                si_dot_n = si[1]*normal[1] + si[2]*normal[2];
                if si_dot_n > 1e-6 # out of bdry and not too close to parallel
                    # Specular (takes dos portion of reflection directly)
                    #how_close = abs(si[1] - sspec[1])+abs(si[2] - sspec[2]);
                    #if how_close < close_to_specular
                    if reflect[dir] == di
                        specular_part = specularity * intensity[di + (band-1)*ndirs] * si_dot_n;
                        close_to_specular = how_close;
                        # if debug println("specular got "*string(specular_part)*" from "*string(di)) end
                    end

                    # Diffuse (takes its share of the (1-dos) portion of reflection)
                    diffuse_part += (1-specularity) * intensity[di + (band-1)*ndirs] * si_dot_n;
                elseif si_dot_n < -1e-6
                    ndiffuse += 1; # count directions taking a share
                end
            end
            # Adjust diffuse part share depending on number of directions taking a share
            diffuse_part /= ndiffuse;
            # if debug println("diffuse got "*string(diffuse_part)*" from "*string(ndiffuse)*" dirs") end
            result = vg[band] * (specular_part + diffuse_part);
        elseif s_dot_n > 1e-6 # out
            result = -vg[band] * intensity[dir + (band-1)*ndirs] * s_dot_n;
            # if debug println("lost "*string(intensity[dir + (band-1)*ndirs] * s_dot_n)) end
        else
            result = 0;
        end

        return result;
    end
)
@callbackFunction(
    function isothermal_bdry(intensity, vg, sx, sy, band, dir, normal, temp)
        interior_intensity = intensity[dir + (band-1)*ndirs];
        iso_intensity = equilibrium_intensity(center_freq[band], delta_freq, temp);
        sdotn = sx[dir]*normal[1] + sy[dir]*normal[2];
        if sdotn > 0
            result = -vg[band] * interior_intensity * sdotn;
        else
            result = -vg[band] * iso_intensity * sdotn;
        end
        return result;
    end
)
@callbackFunction(
    function symmetric_bdry(intensity, vg, sx, sy, band, dir, normal)
        # This is essentially the same as pure specular reflection?
        interior_intensity = intensity[dir + (band-1)*ndirs];
        sym_intensity = intensity[reflect[dir] + (band-1)*ndirs];
        sdotn = sx[dir]*normal[1] + sy[dir]*normal[2];
        if sdotn > 0
            result = -vg[band] * interior_intensity * sdotn;
        else
            result = -vg[band] * sym_intensity * sdotn;
        end
        return result;
    end
)

@callbackFunction(
    function isothermal_bdry_3d(intensity, vg, sx, sy, sz, band, dir, normal, temp)
        interior_intensity = intensity[dir + (band-1)*ndirs];
        iso_intensity = equilibrium_intensity(center_freq[band], delta_freq, temp);
        sdotn = sx[dir]*normal[1] + sy[dir]*normal[2] + sz[dir]*normal[3];
        if sdotn > 0
            result = -vg[band] * interior_intensity * sdotn;
        else
            result = -vg[band] * iso_intensity * sdotn;
        end
        return result;
    end
)
@callbackFunction(
    function symmetric_bdry_3d(intensity, vg, sx, sy, sz, band, dir, normal)
        # This is essentially the same as pure specular reflection?
        interior_intensity = intensity[dir + (band-1)*ndirs];
        sym_intensity = intensity[reflect[dir] + (band-1)*ndirs];
        sdotn = sx[dir]*normal[1] + sy[dir]*normal[2] + sz[dir]*normal[3];
        if sdotn > 0
            result = -vg[band] * interior_intensity * sdotn;
        else
            result = -vg[band] * sym_intensity * sdotn;
        end
        return result;
    end
)

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using Genie.Configuration

const config =  Settings(
                  output_length       = 0,
                  suppress_output     = false,
                  log_db              = false,
                  log_queries         = false,
                  log_requests        = false,
                  log_responses       = false,
                  log_router          = false,
                  log_formatted       = false,
                  log_cache           = false,
                  log_level           = "error",
                  log_verbosity       = LOG_LEVEL_VERBOSITY_MINIMAL,
                  assets_path         = "/",
                  cache_duration      = 1_000,
                  session_auto_start  = false
                )

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
## Transforms take values at Chebyshev points of the first and second kinds and produce Chebyshev coefficients

struct ChebyshevTransformPlan{T,kind,inplace,P} <: Plan{T}
    plan::P
end

ChebyshevTransformPlan{k,inp}(plan) where {k,inp} =
    ChebyshevTransformPlan{eltype(plan),k,inp,typeof(plan)}(plan)

function plan_chebyshevtransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
    if kind == 1
        plan = isempty(x) ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.REDFT10)
        ChebyshevTransformPlan{1,true}(plan)
    elseif kind == 2
        plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.REDFT00)
        ChebyshevTransformPlan{2,true}(plan)
    end
end

function plan_chebyshevtransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
    plan = plan_chebyshevtransform!(x;kind=kind)
    ChebyshevTransformPlan{kind,false}(plan)
end

function *(P::ChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T
    n = length(x)
    n ≤ 1 && return x

    x = P.plan*x
    x[1] /= 2
    lmul!(inv(convert(T,n)), x)
end

function *(P::ChebyshevTransformPlan{T,2,true}, x::AbstractVector{T}) where T
    n = length(x)
    n ≤ 1 && return x

    x = P.plan*x
    x[1] /= 2; x[end] /= 2
    lmul!(inv(convert(T,n-1)),x)
end

chebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
    plan_chebyshevtransform!(x;kind=kind)*x

chebyshevtransform(x;kind::Integer=1) = chebyshevtransform!(copy(x);kind=kind)

*(P::ChebyshevTransformPlan{T,k,false}, x::AbstractVector{T}) where {T,k} = P.plan*copy(x)

## Inverse transforms take Chebyshev coefficients and produce values at Chebyshev points of the first and second kinds

struct IChebyshevTransformPlan{T,kind,inplace,P}
    plan::P
end

# second kind Chebyshev transforms share a plan with their inverse
# so we support this via inv
inv(P::ChebyshevTransformPlan{T,2,true}) where T = IChebyshevTransformPlan{T,2,true,typeof(P)}(P)
inv(P::IChebyshevTransformPlan{T,2,true}) where T = P.plan

\(P::ChebyshevTransformPlan, x::AbstractArray) = inv(P) * x
\(P::IChebyshevTransformPlan, x::AbstractArray) = inv(P) * x

function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
    if kind == 1
        plan = isempty(x) ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.REDFT01)
        IChebyshevTransformPlan{T,1,true,typeof(plan)}(plan)
    elseif kind == 2
        inv(plan_chebyshevtransform!(x;kind=2))
    end
end

function plan_ichebyshevtransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
    plan = plan_ichebyshevtransform!(similar(Vector{T},axes(x));kind=kind)
    IChebyshevTransformPlan{T,kind,false,typeof(plan)}(plan)
end

function *(P::IChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T<:fftwNumber
    isempty(x) && return x
    x[1] *=2
    x = lmul!(convert(T,0.5), P.plan*x)
    x
end

function *(P::IChebyshevTransformPlan{T,2,true},x::AbstractVector{T}) where T<:fftwNumber
    n = length(x)
    n ≤ 1 && return x
    x[1] *= 2; x[end] *= 2
    x = P.plan*x
    x[1] *= 2; x[end] *= 2
    lmul!(convert(T,0.5(n-1)),x)
end

ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
    plan_ichebyshevtransform!(x;kind=kind)*x

ichebyshevtransform(x;kind::Integer=1) = ichebyshevtransform!(copy(x); kind=kind)

*(P::IChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)

## Code generation for integer inputs

for func in (:chebyshevtransform,:ichebyshevtransform)
    @eval $func(x::AbstractVector{T};kind::Integer=1) where {T<:Integer} = $func(convert(Float64,x);kind=kind)
end

# Matrix inputs
#
#
function chebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
    if kind == 1
        if size(X) == (1,1)
            X
        else
            X=FFTW.r2r!(X,FFTW.REDFT10)
            X[:,1]/=2;X[1,:]/=2;
            lmul!(1/(size(X,1)*size(X,2)),X)
        end
    elseif kind == 2
        if size(X) == (1,1)
            X
        else
            X=FFTW.r2r!(X,FFTW.REDFT00)
            lmul!(1/((size(X,1)-1)*(size(X,2)-1)),X)
            X[:,1]/=2;X[:,end]/=2
            X[1,:]/=2;X[end,:]/=2
            X
        end
    end
end
#
function ichebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
    if kind == 1
        if size(X) == (1,1)
            X
        else
            X[1,:]*=2;X[:,1]*=2
            X = FFTW.r2r(X,FFTW.REDFT01)
            lmul!(0.25, X)
        end
    elseif kind == 2
        if size(X) == (1,1)
            X
        else
            X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
            X=chebyshevtransform!(X;kind=kind)
            X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
            lmul!((size(X,1)-1)*(size(X,2)-1)/4,X)
        end
    end
end
#

## Chebyshev U

struct ChebyshevUTransformPlan{T,kind,inplace,P} <: Plan{T}
    plan::P
end

ChebyshevUTransformPlan{k,inp}(plan) where {k,inp} =
    ChebyshevUTransformPlan{eltype(plan),k,inp,typeof(plan)}(plan)

function plan_chebyshevutransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
    if kind == 1
        plan = isempty(x) ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.RODFT10)
        ChebyshevUTransformPlan{1,true}(plan)
    elseif kind == 2
        plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.RODFT00)
        ChebyshevUTransformPlan{2,true}(plan)
    end
end

function plan_chebyshevutransform(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
    plan = plan_chebyshevutransform!(x;kind=kind)
    ChebyshevUTransformPlan{kind,false}(plan)
end

function *(P::ChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T
    n = length(x)
    n ≤ 1 && return x

    for k=1:n # sqrt(1-x_j^2) weight
        x[k] *= sinpi(one(T)/(2n) + (k-one(T))/n)/n
    end
    P.plan * x
end

function *(P::ChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T
    n = length(x)
    n ≤ 1 && return x

    c = one(T)/ (n+1)
    for k=1:n # sqrt(1-x_j^2) weight
        x[k] *= sinpi(k*c)
    end
    lmul!(c, P.plan * x)
end

chebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
    plan_chebyshevutransform!(x;kind=kind)*x

chebyshevutransform(x;kind::Integer=1) = chebyshevutransform!(copy(x);kind=kind)

*(P::ChebyshevUTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)

## Inverse transforms take ChebyshevU coefficients and produce values at ChebyshevU points of the first and second kinds

struct IChebyshevUTransformPlan{T,kind,inplace,P}
    plan::P
end

function plan_ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
    if kind == 1
        plan = isempty(x) ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.RODFT01)
        IChebyshevUTransformPlan{T,1,true,typeof(plan)}(plan)
    elseif kind == 2
        plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : FFTW.plan_r2r!(x, FFTW.RODFT00)
        IChebyshevUTransformPlan{T,2,true,typeof(plan)}(plan)
    end
end

function plan_ichebyshevutransform(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
    plan = plan_ichebyshevutransform!(similar(Vector{T},axes(x)); kind=kind)
    IChebyshevUTransformPlan{T,kind,false,typeof(plan)}(plan)
end

function *(P::IChebyshevUTransformPlan{T,1,true}, x::AbstractVector{T}) where T<:fftwNumber
    n = length(x)
    n ≤ 1 && return x

    x = P.plan * x
    for k=1:n # sqrt(1-x_j^2) weight
        x[k] /= 2sinpi(one(T)/(2n) + (k-one(T))/n)
    end
    x
end

function *(P::IChebyshevUTransformPlan{T,2,true}, x::AbstractVector{T}) where T<:fftwNumber
    n = length(x)
    n ≤ 1 && return x

    c = one(T)/ (n+1)
    lmul!((n+1)/(2n+2*one(T)), x)
    x = P.plan * x
    for k=1:n # sqrt(1-x_j^2) weight
        x[k] /= sinpi(k*c)
    end
    x
end

ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
    plan_ichebyshevutransform!(x;kind=kind)*x

ichebyshevutransform(x;kind::Integer=1) = ichebyshevutransform!(copy(x);kind=kind)

*(P::IChebyshevUTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)

## Code generation for integer inputs

for func in (:chebyshevutransform,:ichebyshevutransform)
    @eval $func(x::AbstractVector{T};kind::Integer=1) where {T<:Integer} = $func(convert(Float64,x);kind=kind)
end

## points

function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
    if kind == 1
        T[sinpi((n-2k-one(T))/2n) for k=0:n-1]
    elseif kind == 2
        if n == 1
            zeros(T,1)
        else
	    T[sinpi((n-2k-one(T))/(2n-2)) for k=0:n-1]
        end
    else
        throw(ArgumentError("kind $kind not a valid kind of Chebyshev points"))
    end
end
chebyshevpoints(n::Integer; kind::Int=1) = chebyshevpoints(Float64, n; kind=kind)

# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
#
# struct DSTPlan{T,kind,inplace,P} <: Plan{T}
#     plan::P
# end
#
# DSTPlan{k,inp}(plan) where {k,inp} =
#     DSTPlan{eltype(plan),k,inp,typeof(plan)}(plan)
#
#
# plan_DSTI!(x) = length(x) > 0 ? DSTPlan{1,true}(FFTW.FFTW.plan_r2r!(x, FFTW.FFTW.RODFT00)) :
#                                 fill(one(T),1,length(x))
#
# function *(P::DSTPlan{T,1}, x::AbstractArray) where {T}
#     x = P.plan*x
#     rmul!(x,half(T))
# end