text:
"""
    coordinates!(mtg; angle = 45; force = false)

Compute dummy 3d coordinates for the mtg nodes using an alterning phyllotaxy. Used when
coordinates are missing.
Coordinates are just node attributes with reserved names: :XX, :YY and :ZZ.

# Returns

Nothing, mutates the mtg in-place (adds :XX, :YY and :ZZ to nodes).

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
coordinates!(mtg)
DataFrame(mtg, [:XX, :YY, :ZZ])
```
"""
function coordinates!(mtg; angle = 45, force = false)

    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))

    if !force && any(coord_in_attributes)
        error(
            "Coordinates for $([:XX, :YY, :ZZ][coord_in_attributes]) are already present",
            " in the mtg attributes (at least for some nodes). Use `force = true` if you ",
            "want to overwrite their values."
        )
    end

    phyllotaxy = [-1]
    # This function adds a XX, YY and ZZ coordinates of each node
    traverse!(mtg, new_pos, angle, phyllotaxy)
end

function coordinates_parent!(mtg)
    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))
    if !all(coord_in_attributes)
        error("No coordinates found in MTG. Use `coordinates!(mtg)` first.")
    end

    transform!(
        mtg,
        :XX => (x -> 1.0) => :ZZ,
        (x -> get(ancestors(x, :XX, recursivity_level = 1), 1, nothing)) => :XX_from,
        (x -> get(ancestors(x, :YY, recursivity_level = 1), 1, nothing)) => :YY_from,
        (x -> get(ancestors(x, :ZZ, recursivity_level = 1), 1, nothing)) => :ZZ_from
    )
end

function new_pos(node, angle, phyllotaxy)

    if isroot(node)
        append!(node, (XX = 0.0, YY = 0.0, ZZ = 0.0))
        return nothing
    end

    parent_node = parent(node)
    great_parent_node = parent(parent_node)

    if great_parent_node === nothing
        great_parent_node_XX = 0.0
        great_parent_node_YY = -1.0
    else
        great_parent_node_XX = great_parent_node[:XX]
        great_parent_node_YY = great_parent_node[:YY]
    end

    if node.MTG.link == "/"
        extend_length = 0.2
    else
        extend_length = 1
    end

    point = extend_pos(great_parent_node_XX, great_parent_node_YY, parent_node[:XX], parent_node[:YY], extend_length)

    if node.MTG.link == "+"
        point =
            rotate_point(
                parent_node[:XX],
                parent_node[:YY],
                point[1],
                point[2],
                phyllotaxy[1] * angle
            )

        # Change phyllotaxy for next node:
        if phyllotaxy[1] == 1
            phyllotaxy[1] = -1
        else
            phyllotaxy[1] = 1
        end
    end
    node[:XX] = point[1]
    node[:YY] = point[2]
    node[:ZZ] = 0.0

    return
end

"""
Add a new point after (x1,y1) using same direction and length relative to it
"""
function extend_pos(x0, y0, x1, y1, extend_length)
    lengthAB = sqrt((x0 - x1)^2 + (y0 - y1)^2)
    x = x1 + (x1 - x0) / lengthAB * extend_length
    y = y1 + (y1 - y0) / lengthAB * extend_length
    return x, y
end

"""
Rotate a point (x1,y1) around (x0, y0) with `angle`.
"""
function rotate_point(x0, y0, x1, y1, angle)
    angle = -angle * pi / 180
    x1 = x1 - x0
    y1 = y1 - y0
    cos_a = cos(angle)
    sin_a = sin(angle)

    x = x1 * cos_a - y1 * sin_a + x0
    y = x1 * sin_a + y1 * cos_a + y0

    return x, y
end

"""
    mtg_coordinates_df(mtg, attr; force = false)
    mtg_coordinates_df!(mtg, attr; force = false)

Extract the coordinates of the nodes of the mtg and the
coordinates of their parents (:XX_from, :YY_from, :ZZ_from) and output
a DataFrame. Optionally you can also provide an attribute to add to the
output DataFrame too by passing its name as a symbol to `attr`.

The coordinates are computed using [`coordinates!`](@ref) if missing, or if
`force = true`.
"""
function mtg_coordinates_df(mtg, attr = :YY; force = false)
    mtg_coordinates_df!(deepcopy(mtg), attr; force = force)
end

function mtg_coordinates_df!(mtg, attr = :YY; force = false)
    coord_in_attributes = [:XX, :YY, :ZZ] .∈ Ref(get_attributes(mtg))
    if !all(coord_in_attributes) || force
        if coord_in_attributes == [true, true, false] && !force
            # Only :ZZ is missing, using a dummy value:
            transform!(
                mtg,
                :XX => (x -> 0.0) => :ZZ,
            )
        else
            # :XX and/or :YY missing, recompute them
            coordinates!(mtg; force = force)
        end
    end

    # Get the coordinates of the parent of the node to draw edges:
    coordinates_parent!(mtg)

    DataFrame(mtg, unique([:XX, :YY, :ZZ, :XX_from, :YY_from, :ZZ_from, attr]))
end

function mtg_XYZ_color(mtg, color, edge_color, colormap; color_missing = RGBA(0, 0, 0, 0.3))
    if Symbol(color) in get_attributes(mtg)
        # Color is an attribute from the mtg:
        df_coordinates = mtg_coordinates_df(mtg, color, force = true)
        max_color = maximum(skipmissing(df_coordinates[:, color]))
        color_var = [i === missing ? color_missing : RGBA(get(colormap, i / max_color)) for i in df_coordinates[:, color]]
        text_color = color_var
    elseif typeof(color) <: Colorant || typeof(color) <: String || typeof(color) <: Symbol
        df_coordinates = mtg_coordinates_df(mtg, force = true)
        color_var = fill(color, size(df_coordinates, 1))
        text_color = color_var
    else
        error(
            "color argument should be of type Colorant ",
            "(see [Colors.jl](https://juliagraphics.github.io/Colors.jl/stable/)), or ",
            "an MTG attribute."
        )
    end

    if Symbol(edge_color) in get_attributes(mtg)
        if edge_color == color
            # If edge_color use the same attribute than color, use color_var (avoid recomputing)
            edge_color_var = [[color_var[i-1], color_var[i]] for i in 2:length(color_var)]
            pushfirst!(edge_color_var, edge_color_var[1])
        else
            # If edge_color is different than color, and is an attribute from the mtg
            df_coordinates = mtg_coordinates_df(mtg, edge_color, force = true)

            ecol = df_coordinates[:, edge_color]
            max_edge_color = maximum(skipmissing(ecol))
            edge_color_var = []
            for i in 2:size(df_coordinates, 1)
                if ecol[i-1] === missing || ecol[i] === missing
                    push!(edge_color_var, color_missing)
                else
                    push!(edge_color_var, RGBA.(get(colormap, [ecol[i-1] / max_edge_color, ecol[i] / max_edge_color])))
                end
            end
            pushfirst!(edge_color_var, edge_color_var[1])
        end
    elseif typeof(edge_color) <: Colorant || typeof(edge_color) <: String || typeof(edge_color) <: Symbol
        # edge_color is just a single color
        edge_color_var = fill([edge_color, edge_color], size(df_coordinates, 1))
    else
        error(
            "edge_color argument should be of type Colorant ",
            "(see [Colors.jl](https://juliagraphics.github.io/Colors.jl/stable/)), or ",
            "an MTG attribute."
        )
    end

    return df_coordinates, color_var, edge_color_var, text_color
end

function attribute_range(mtg, attr)
    vals =
        descendants(
            mtg,
            attr,
            ignore_nothing = true
        )

    if length(vals[1]) == 1
        range_val = extrema(vals)
    else
        range_val = (minimum(minimum.(vals)), maximum(maximum.(vals)))
    end

    return range_val
end


text:
#EasyPlot: Tools to manipulate text
#-------------------------------------------------------------------------------
#=TODO
 - Find a more practical font system than InspectDR.
=#

#==Useful constants
===============================================================================#

#==Main types
===============================================================================#
abstract type AbstractFont end

#=
mutable struct Font <: AbstractFont
	name::String
	_size::Float64
	bold::Bool
	color::Colorant
end
Font(name::String, _size::Real; bold::Bool=false, color=COLOR_BLACK) =
	Font(name, _size, bold, color)
Font(::PreDefaultsType) = Font(DEFAULT_FONTNAME, 10) #Construct some object
=#

ext: jl


text:
# Intended to ensure the README stays working (this is a copy)

using AxisArrays, Unitful
import Unitful: s, ms, µs

fs = 40000
y = randn(60*fs+1)*3
for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
           sin.(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
    i = rand(round(Int,.001fs):1fs)
    while i+length(spk)-1 < length(y)
        y[i:i+length(spk)-1] += spk
        i += rand(round(Int,.001fs):1fs)
    end
end

A = AxisArray([y 2y], Axis{:time}(0s:1s/fs:60s), Axis{:chan}([:c1, :c2]))
A = AxisArray(hcat(y, 2 .* y); time = (0s:1s/fs:60s), chan = ([:c1, :c2]))

A[Axis{:time}(4)]
A[time=4]

A[Axis{:chan}(:c2), Axis{:time}(1:5)]
A[chan = :c2, time = 1:5]

ax = A[40µs .. 220µs, :c1]
AxisArrays.axes(ax, 1)
A[atindex(-90µs .. 90µs, 5), :c2]
idxs = findall(diff(A[:,:c1] .< -15) .> 0)
spks = A[atindex(-200µs .. 800µs, idxs), :c1]
A[atvalue(2.5e-5s), :c1]
A[2.5e-5s..2.5e-5s, :c1]
A[atvalue(25.0µs)]

# # A possible "dynamic verification" strategy
# const readmefile = joinpath(dirname(dirname(@__FILE__)), "README.md")

# function extract_julialines(iowr, filein)
#     open(filein) do iord
#         while !eof(iord)
#             line = readline(iord)
#             if startswith(line, "julia>")
#                 print(iowr, line[8:end])
#                 while !eof(iord)
#                     line = readline(iord)
#                     if !startswith(line, "    ")
#                         break
#                     end
#                     print(iowr, line[8:end])
#                 end
#             end
#         end
#     end
# end

# tmpfile, iowr = mktemp()
# extract_julialines(iowr, readmefile)
# close(iowr)

# include(tmpfile)
# rm(tmpfile)


text:
"Tableau for variational Gauss-Lobatto IIIA-IIIB method with two stages"
function getTableauVPRKpLegendreQLobIIIA2()
    ω = [0.5, 0.5]
    ω = reshape(ω_λ, 1, length(ω_λ))
    d = [+1.0, -1.0]
    R∞ = -1
    TableauVPRKpLegendreQ(:LobIIIA2, 2, getCoefficientsLobIIIA2(), R∞, ω, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with three stages"
function getTableauVPRKpLegendreQLobIIIA3()
    d = [+1.0, -2.0, +1.0]
    R∞ = +1
    TableauVPRKpLegendreQ(:LobIIIA3, 4, getCoefficientsLobIIIA3(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with four stages"
function getTableauVPRKpLegendreQLobIIIA4()
    d = [+1.0, -√5, +√5, -1.0]
    R∞ = -1
    TableauVPRKpLegendreQ(:LobIIIA4, 6, getCoefficientsLobIIIA4(), R∞, d)
end

"Tableau for variational Gauss-Legendre method with s stages"
function getTableauVPRKpLegendreQGLRK(s; T=Float64)
    glrk = getCoefficientsGLRK(s, T=T)
    R∞ = (-1)^s
    TableauVPRKpLegendreQ(Symbol("GLRK", s), glrk.o, glrk, R∞)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with two stages"
function getTableauVPRKpLegendrePLobIIIA2()
    d = [+1.0, -1.0]
    R∞ = -1
    TableauVPRKpLegendreP(:LobIIIA2, 2, getCoefficientsLobIIIA2(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with three stages"
function getTableauVPRKpLegendrePLobIIIA3()
    d = [+1.0, -2.0, +1.0]
    R∞ = +1
    TableauVPRKpLegendreP(:LobIIIA3, 4, getCoefficientsLobIIIA3(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with four stages"
function getTableauVPRKpLegendrePLobIIIA4()
    d = [+1.0, -√5, +√5, -1.0]
    R∞ = -1
    TableauVPRKpLegendreP(:LobIIIA4, 6, getCoefficientsLobIIIA4(), R∞, d)
end

"Tableau for variational Gauss-Legendre method with s stages"
function getTableauVPRKpLegendrePGLRK(s; T=Float64)
    glrk = getCoefficientsGLRK(s, T=T)
    R∞ = (-1)^s
    TableauVPRKpLegendreP(Symbol("GLRK", s), glrk.o, glrk, R∞)
end


text:
"""
    initial_conditions(model)

Get initial conditions at all degrees of freedom.
"""
initial_conditions(model) =
    mapreduce((ρ, p) -> fill(ρ, size(p, 2)), vcat, model.ρ, model.mesh.points)


text:
using RadeonProRender, GeometryBasics, Colors
const RPR = RadeonProRender

function translationmatrix(t::Vec{3,T}) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1, T0, T0, T0, T0, T1, T0, T0, T0, T0, T1, T0, t[1], t[2], t[3], T1)
end

context = RPR.Context()
scene = RPR.Scene(context)
matsys = RPR.MaterialSystem(context, 0)

camera = RPR.Camera(context)
lookat!(camera, Vec3f(8, 0, 5), Vec3f(2, 0, 2), Vec3f(0, 0, 1))
RPR.rprCameraSetFocalLength(camera, 45.0)
set!(scene, camera)
set!(context, scene)

env_light = RadeonProRender.EnvironmentLight(context)
image_path = RPR.assetpath("studio026.exr")
img = RPR.Image(context, image_path)
set!(env_light, img)
setintensityscale!(env_light, 1.1)
push!(scene, env_light)

light = RPR.PointLight(context)
transform!(light, translationmatrix(Vec3f(2, 0, 5)))
f = 16
RPR.setradiantpower!(light, 700 / f, 641 / f, 630 / f)
push!(scene, light)

function add_shape!(scene, context, matsys, mesh; material=RPR.RPR_MATERIAL_NODE_DIFFUSE,
                    color=colorant"green", roughness=0.01)
    rpr_mesh = RPR.Shape(context, mesh)
    push!(scene, rpr_mesh)
    m = RPR.MaterialNode(matsys, material)
    set!(m, RPR.RPR_MATERIAL_INPUT_COLOR, color)
    set!(m, RPR.RPR_MATERIAL_INPUT_ROUGHNESS, roughness, roughness, roughness, roughness)
    set!(rpr_mesh, m)
    return rpr_mesh, m
end

mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(-10, -10, -1), Vec3f(20, 20, 1));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, -10, 0), Vec3f(0.1, 20, 5));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, -2, 0), Vec3f(5, 0.1, 5));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, 2, 0), Vec3f(5, 0.1, 5));
                       color=colorant"white")

mesh, mat = add_shape!(scene, context, matsys, Tesselation(Sphere(Point3f(2, 0, 2), 0.95f0), 100);
                       color=colorant"white")
rpr_mesh = RPR.Shape(context, Tesselation(Sphere(Point3f(2, 0, 2), 1.0f0), 100))
push!(scene, rpr_mesh)
material = RPR.MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_UBERV2)

# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_COLOR, 0.358878, 0.497024 , 1.07771, 1.0)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_WEIGHT, Vec4f(0.0)...)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_ROUGHNESS, Vec4f(0.1)...)

# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_COLOR, 0.01f0, 0.01f0, 0.01f0, 0.01f0)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_WEIGHT, Vec4f(0.01)...)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_TRANSPARENCY, Vec4f(0.99)...)

set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_COLOR, Vec4f(0.358878, 0.497024 , 1.07771, 1.0)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_WEIGHT, Vec4f(0.8)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ROUGHNESS, Vec4f(0.01)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_IOR, Vec4f(1.0)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_COLOR, Vec4f(0.7, 0.6, 0.09, 1.0)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_DISTANCE, Vec4f(3.0)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_CAUSTICS, Vec4f(0.0)...)

# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ROUGHNESS, 0.1f0, 0f0, 0f0, 0f0)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_MODE, Int(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR))
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_METALNESS, 0.f0, 0f0, 0f0, 0f0)

set!(material, RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_COLOR, 0.358878, 0.497024 , 1.07771, 1.0)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DISTANCE, Vec4f(3)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DIRECTION, Vec4f(0)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_SSS_WEIGHT, Vec4f(1)...)
set!(material, RPR.RPR_MATERIAL_INPUT_UBER_SSS_MULTISCATTER, 1)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_WEIGHT, Vec4f(1.0)...)
# set!(material, RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_COLOR, 0.4, 0.4, 0.5, 0.1)

set!(rpr_mesh, material)

fb_size = (800, 600)
frame_buffer = RPR.FrameBuffer(context, RGBA, fb_size)
frame_bufferSolved = RPR.FrameBuffer(context, RGBA, fb_size)
set!(context, RPR.RPR_AOV_COLOR, frame_buffer)

set_standard_tonemapping!(context)

begin
    clear!(frame_buffer)
    RPR.rprContextSetParameterByKey1u(context, RPR.RPR_CONTEXT_ITERATIONS, 200)
    RPR.render(context)
    RPR.rprContextResolveFrameBuffer(context, frame_buffer, frame_bufferSolved, false)
    RPR.save(frame_bufferSolved, "test.png")
end


text:
"""
    ReLU(θ, γ)

ReLU layer, with location parameters `θ` and scale parameters `γ`.
"""
struct ReLU{N,A1,A2} <: AbstractLayer{N}
    θ::A1
    γ::A2
    function ReLU(θ::AbstractArray, γ::AbstractArray)
        @assert size(θ) == size(γ)
        return new{ndims(θ), typeof(θ), typeof(γ)}(θ, γ)
    end
end

ReLU(::Type{T}, sz::Int...) where {T} = ReLU(zeros(T, sz...), ones(T, sz...))
ReLU(sz::Int...) = ReLU(Float64, sz...)

Base.repeat(l::ReLU, n::Int...) = ReLU(repeat(l.θ, n...), repeat(l.γ, n...))

function energies(layer::ReLU, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    return relu_energy.(layer.θ, layer.γ, x)
end

free_energies(layer::ReLU, inputs::Union{Real,AbstractArray} = 0) = relu_free.(layer.θ .+ inputs, layer.γ)
sample_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0) = relu_rand.(layer.θ .+ inputs, layer.γ)
mode_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0) = max.((layer.θ .+ inputs) ./ abs.(layer.γ), 0)

function mean_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0)
    g = Gaussian(layer.θ, layer.γ)
    μ = mean_from_inputs(g, inputs)
    σ = sqrt.(var_from_inputs(g, inputs))
    return @. μ + σ * tnmean(-μ / σ)
end

mean_abs_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0) = mean_from_inputs(layer, inputs)

function var_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0)
    g = Gaussian(layer.θ, layer.γ)
    μ = mean_from_inputs(g, inputs)
    ν = var_from_inputs(g, inputs)
    return @. ν * tnvar(-μ / √ν)
end

function meanvar_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0)
    g = Gaussian(layer.θ, layer.γ)
    μ = mean_from_inputs(g, inputs)
    ν = var_from_inputs(g, inputs)
    σ = sqrt.(ν)
    tμ, tν = tnmeanvar(-μ ./ σ)
    return μ + σ .* tμ, ν .* tν
end

std_from_inputs(layer::ReLU, inputs::Union{Real,AbstractArray} = 0) = sqrt.(var_from_inputs(layer, inputs))

function ∂free_energies(layer::ReLU, inputs::Union{Real,AbstractArray} = 0)
    μ, ν = meanvar_from_inputs(layer, inputs)
    return (θ = -μ, γ = sign.(layer.γ) .* (ν .+ μ.^2) / 2)
end

const ReLUStats = GaussStats
suffstats(layer::ReLU, data::AbstractArray; wts = nothing) = ReLUStats(layer, data; wts)

function ∂energy(layer::ReLU, stats::ReLUStats)
    @assert size(layer) == size(stats)
    return (; θ = -stats.x1, γ = sign.(layer.γ) .* stats.x2 / 2)
end

function relu_energy(θ::Real, γ::Real, x::Real)
    E = gauss_energy(θ, γ, x)
    return x < 0 ? inf(E) : E
end

function relu_free(θ::Real, γ::Real)
    abs_γ = abs(γ)
    return -logerfcx(-θ / √(2abs_γ)) + log(2abs_γ/π) / 2
end

function relu_rand(θ::Real, γ::Real)
    abs_γ = abs(γ)
    μ = θ / abs_γ
    σ = √inv(abs_γ)
    return randnt_half(μ, σ)
end


text:
module TestOIFITS

using OIFITS
using Compat
using Compat.Test
using Compat: @debug, @error, @info, @warn

dir = @__DIR__

files = ("contest-2004-obj1.oifits" ,"contest-2004-obj2.oifits",
         "contest-2008-binary.oifits", "contest-2008-obj1-H.oifits",
         "contest-2008-obj1-J.oifits", "contest-2008-obj1-K.oifits",
         "contest-2008-obj2-H.oifits", "contest-2008-obj2-J.oifits",
         "contest-2008-obj2-K.oifits", "contest-2008-obj1-J.oifits")

counter = 0

quiet = true

function tryload(dir, file)
    global counter
    try
        db = OIFITS.load(joinpath(dir, file))
        counter += 1
        quiet || @info "file \"", file, "\" successfully loaded"
        return true
    catch
        @warn "failed to load \"", file, "\""
        return false
    end
end

@testset "OIFITS.load" begin
    for file in files
        @test tryload(dir, file) == true
    end
end

end

nothing


text:
"""
测试温度流的refine
"""

using Plots
using MARY_fRG.Drawers
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

"""
获取极低温的时候的bubble
"""
function get_ult_bubb(Γ4; maxnum=1000000)
    Γ4_ult = fliter_away_surface(Γ4, 2.0)
    Γ4_ult = refine_to_surface(Γ4_ult)
    while length(Γ4_ult.ltris) < maxnum
        Γ4_ult = refine_to_surface(Γ4_ult)
    end
    plt = draw_points([tri.center for tri in Γ4_ult.ltris])
    @savepng plt "ult_surface"
    maxi, mini = engpeak_to_surface(Γ4_ult)
    bubb_pp_u, bubb_fs_u, bubb_ex_u = all_bubble_tf_ult_mt(Γ4_ult, maxi)
    return bubb_pp_u, bubb_fs_u, bubb_ex_u
end

function run_rftf()
    model = common_square_lattice(0.20)
    Γ4 = TFGamma4(
        model, 8.0, 16, 200
    )
    Γ4.V .+= 1.0
    lval = 0.
    lstep = 0.01
    #plt = draw_points([tri.center for tri in Γ4.ltris])
    #@savepng plt "brlu"*string(lval)
    #
    bpu, bfu, beu = get_ult_bubb(Γ4)
    @info "ult结束"
    #
    for idx in 1:1:3351
        #if idx % 10 == 1
        #    blval = lval - 5
        #    blval = max(blval, 0.)
        #    blval = min(blval, 5.)
        #    Γ4b = TFGamma4_addition_ltris_mt(Γ4, blval; maxnum=1000000)
        #    println(length(Γ4.ltris), "->", length(Γ4b.ltris))
        #    #plt = draw_points([tri.center for tri in Γ4b.ltris])
        #    #@savepng plt "brlu"*string(lval)
        #end
        #blval = min(lval, 10.0)
        bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, lval)
        dl = dl_tf_mix_ult_mt(Γ4, bubb_pp, bubb_fs, bubb_ex, bpu, bfu, beu)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1 || idx > 3301
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma4"*string(idx))
        end
        lval += lstep
    end
end

run_rftf()


text:
module VersionBenchmarks

using Dates
using DataFrames
using Statistics: mean, std, median
import AlgebraOfGraphics
const AoG = AlgebraOfGraphics
import CairoMakie
import Pkg

export Config
export benchmark

const hide_output = Ref(true)

function __init__()
    CairoMakie.activate!(px_per_unit = 2)
end

function _run(cmd)
    if hide_output[]
        run(pipeline(cmd, stderr = devnull, stdout = devnull))
    else
        run(cmd)
    end
end

struct Config
    name::String
    pkgspecs::Vector{NamedTuple}
    julia::Cmd
end

Config(name, pkgspec::NamedTuple, julia) = Config(name, [pkgspec], julia)
Config(name, pkgspecs) = Config(name, pkgspecs, `julia`)

benchmark(configs, file::String; kwargs...) = benchmark(configs, [file]; kwargs...)

function benchmark(configs::Vector{Config}, files::AbstractVector{String};
        repetitions = 1
    )

    tmpdir = mktempdir()

    df = DataFrame()

    date = now()

    tmpdir_dict = Dict()

    try
        # repetition should be the outer loop so that neither julia version nor code
        # version is blocked in time, which should give better robustness against
        # performance fluctuations of the system over time
        n_runs = repetitions * length(configs)
        start = time()
        i = 0
        for repetition in 1:repetitions
            for (iconfig, config) in enumerate(configs)
                i += 1
                julia_cmd = config.julia
                julia_version = get_julia_version(julia_cmd)

                println("""
                    Repetition $repetition of $repetitions, config $iconfig of $(length(configs)).
                     ├ Name: $(config.name)
                     └ Julia: $julia_version""")

                tmpenvdir = prepare_julia_environment(config, tmpdir_dict)

                for file in files
                    time_of_run = now()

                    print("   Executing \"$file\"...")
                    t = time()
                    resultdf = execute_file(file, julia_cmd, tmpenvdir, repetition)
                    println(" Done. ($(round(Int, time() - t))s)")

                    duration_of_run = now() - time_of_run
                    resultdf.config_name .= config.name
                    resultdf.pkgspecs .= Ref(config.pkgspecs)
                    resultdf.file .= file
                    resultdf.date .= date
                    resultdf.time_of_run .= time_of_run
                    resultdf.duration_of_run .= duration_of_run
                    resultdf.repetition .= repetition
                    resultdf.julia_version .= julia_version
                    append!(df, resultdf, cols = :union)
                end

                elapsed = time() - start
                estimated = elapsed / i * (n_runs - i)
                println("   Time elapsed: $(round(Int, elapsed))s.")
                println("   Estimated remaining: $(round(Int, estimated))s.")
                println("   Estimated finish time: $(Dates.format(now() + Second(round(Int, estimated)), "HH:MM")).")
            end
        end
    catch e
        rethrow(e)
    finally
        rm(tmpdir, recursive = true)
    end

    return df
end

get_julia_version(julia_cmd) = read(`$julia_cmd --startup-file=no -e 'print(VERSION)'`, String)

function prepare_julia_environment(config, tmpdir_dict)
    if haskey(tmpdir_dict, config)
        tmpenvdir = tmpdir_dict[config]

        code = """
        "@stdlib" ∉ LOAD_PATH && push!(LOAD_PATH, "@stdlib")
        using Pkg
        Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true # avoid updating every time
        Pkg.activate("$tmpenvdir")
        Pkg.precompile()
        """
    else
        # create a directory for the environment in which to install the version
        tmpenvdir = mktempdir()
        tmpdir_dict[config] = tmpenvdir

        specs = [config.pkgspecs; (;name = "BenchmarkTools")]
        specstring = string(specs) # TODO: other way to transfer package specs to the new process?

        code = """
        "@stdlib" ∉ LOAD_PATH && push!(LOAD_PATH, "@stdlib")
        using Pkg
        Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true # avoid updating every time
        Pkg.activate("$tmpenvdir")
        Pkg.add($specstring; $(hide_output[] ? "io = devnull, " : ""))
        Pkg.precompile()
        """
    end
    print("   Preparing Julia environment...")
    t = time()
    _run(`$(config.julia) --startup-file=no -e $code`)
    println(" Done. ($(round(Int, time() - t))s)")
    return tmpenvdir
end

function execute_file(file, julia_cmd, tmpenvdir, repetition)
    resultpath, resultio = mktemp()

    basecode = get_basecode(tmpenvdir, resultpath, repetition)

    testcode = read(file, String)

    fullcode = join([basecode, testcode], "\n")

    path, io = mktemp()
    open(path, "w") do file
        println(file, fullcode)
    end

    # execute the modified code, this should write results to the temp file at `resultpath`
    _run(`$julia_cmd --startup-file=no $path`)

    results = foldl(((a, b) -> push!(a, b, cols = :union)),
        map(parse_result_line, readlines(resultpath)),
        init = DataFrame())
    return results
end

function parse_result_line(line)
    # the file should only have lines with serialized NamedTuples
    lineresult = Dict(pairs(eval(Meta.parse(line))))
end

function get_basecode(tmpenvdir, resultpath, repetition)
    """
    "@stdlib" ∉ LOAD_PATH && push!(LOAD_PATH, "@stdlib")
    using Pkg
    Pkg.activate("$tmpenvdir")
    import BenchmarkTools
    import Statistics

    macro vbtime(name, expr)
        @assert name isa String
        quote
            io = open("$resultpath", "a")
            @timed 1 + 1
            tstart = time_ns()
            timed = @timed(\$(esc(expr)))
            dt = (time_ns() - tstart) / 1_000_000_000
            println(io, (type = "vbtime", name = \$name, time_s = dt, alloc_bytes = timed.bytes, gctime_s = timed.gctime))
            close(io)
        end
    end

    macro vbbenchmark(name, exprs...)
        @assert name isa String
        quote
            # only run vbbenchmark on first repetition
            if $repetition == 1
                io = open("$resultpath", "a")
                bm = BenchmarkTools.@benchmark \$(exprs...)
                println(io, (
                    type = "vbbenchmark",
                    name = \$name,
                    min_time_ns = minimum(bm.times),
                    max_time_ns = maximum(bm.times),
                    median_time_ns = Statistics.median(bm.times),
                    mean_time_ns = Statistics.mean(bm.times),
                ))
                close(io)
            end
        end
    end
"""
end

function summarize_repetitions(df)
    gdf = groupby(df, [:config_name, :name, :julia_version])
    combine(gdf, :repetition => (x -> length(x)) => :n, [:time_s :alloc_bytes :gctime_s] .=> [mean, std, minimum, maximum])
end

function plot_summary(df, variable = :time_s)
    df = dropmissing(df, variable)
    plt = AoG.data(df) *
        AoG.mapping(
            :config_name, variable;
            row = :name,
            col = :julia_version => x -> "Julia $x",
        ) *
        (AoG.expectation() + AoG.visual(strokewidth = 1, strokecolor = :black, color = :transparent))

    AoG.draw(plt,
        facet = (; linkyaxes = :minimal),
        axis = (;
            limits = (nothing, nothing, 0, nothing),
            xticklabelrotation = pi/6,
        ),
    )
end

end


text:
# From Smallwood and Sondik, 1973
# 3 states       - 0, 1 or 2 faulty internal components
# 4 actions
#    - manufacture without examination (M)
#    - manufacture with examination (E)
#    - interrupt the line and replace failed components (I)
#    - interrupt the line and replace both components without inspection (R)
# 2 observations - defective product or not
#
# NOTE: The original paper includes "salvage values", or terminal rewards
#       equal to the number of working parts. We do not include that here.

# Component breakdowns are independent
# 10% probability that either breaks down during the manufacture of a product
# Failed components cause 50% of part being defective

@with_kw struct MachineReplacement
    γ::Float64=1.0
end

function DiscretePOMDP(pomdp::MachineReplacement; γ::Float64=pomdp.γ)
    T = Array{Float64}(undef, 3, 4, 3)
    R = Array{Float64}(undef, 3, 4)
    O = Array{Float64}(undef, 4, 3, 2)

    s_0 = 1 # none broken
    s_1 = 2 # one broken
    s_2 = 3 # two broken

    a_m = 1 # manufacture
    a_e = 2 # manufacture + examine
    a_i = 3 # interrupt the line, inspect components, replace failed components
    a_r = 4 # interrupt the line, replace both components

    o_n = 1 # nondefective
    o_d = 2 # defective

    T[s_0, a_m, :] = [0.81, 0.18, 0.01] # 10% independent chance of part breaking
    T[s_0, a_e, :] = [0.81, 0.18, 0.01]
    T[s_0, a_i, :] = [1.00, 0.00, 0.00]
    T[s_0, a_r, :] = [1.00, 0.00, 0.00]
    T[s_1, a_m, :] = [0.00, 0.90, 0.10] # 10% chance of remaining part breaking
    T[s_1, a_e, :] = [0.00, 0.90, 0.10]
    T[s_1, a_i, :] = [1.00, 0.00, 0.00]
    T[s_1, a_r, :] = [1.00, 0.00, 0.00]
    T[s_2, a_m, :] = [0.00, 0.00, 1.00] # stay broken
    T[s_2, a_e, :] = [0.00, 0.00, 1.00]
    T[s_2, a_i, :] = [1.00, 0.00, 0.00]
    T[s_2, a_r, :] = [1.00, 0.00, 0.00]

    # There is a profit of 1 for producing a nondefective product.
    # Thus, the expected profit for beginning with 0, 1, or 2 defective parts is
    # 0.9025, 0.475, and 0.25, respectively.
    # Examining the finished product costs 0.25.
    # The inspect action incurs a 0.5 penalty plus replacement cost for each unit of 1.
    # The straight-up replacement action has no inspection cost but does incur a 2 unit cost.
    r_examine = -0.25
    r_inspect = -0.5
    r_replace = -2.0

    R[s_0, a_m] = 0.9025
    R[s_1, a_m] = 0.475
    R[s_2, a_m] = 0.25
    R[s_0, a_e] = 0.9025 + r_examine
    R[s_1, a_e] = 0.475 + r_examine
    R[s_2, a_e] = 0.25 + r_examine
    R[s_0, a_i] = r_inspect
    R[s_1, a_i] = r_inspect - 1.0 # replace 1
    R[s_2, a_i] = r_inspect - 2.0 # replace 2
    R[s_0, a_r] = r_replace
    R[s_1, a_r] = r_replace
    R[s_2, a_r] = r_replace

    # Probabilities of observing a nondefective product are 1.0, 0.5, and 0.25 if
    # there are 0, 1, or 2 faulty internal components.
    # If we don't examine, we always observe nondefective.
    O[a_m, s_0, :] = [1.00, 0.00]
    O[a_m, s_1, :] = [1.00, 0.00]
    O[a_m, s_2, :] = [1.00, 0.00]
    O[a_e, s_0, :] = [1.00, 0.00]
    O[a_e, s_1, :] = [0.50, 0.50]
    O[a_e, s_2, :] = [0.25, 0.75]
    O[a_i, s_0, :] = [1.00, 0.00]
    O[a_i, s_1, :] = [1.00, 0.00]
    O[a_i, s_2, :] = [1.00, 0.00]
    O[a_r, s_0, :] = [1.00, 0.00]
    O[a_r, s_1, :] = [1.00, 0.00]
    O[a_r, s_2, :] = [1.00, 0.00]

    return DiscretePOMDP(T, R, O, γ)
end

function POMDP(pomdp::MachineReplacement; γ::Float64=pomdp.γ)
    T = Array{Float64}(undef, 3, 4, 3)
    R = Array{Float64}(undef, 3, 4)
    O = Array{Float64}(undef, 4, 3, 2)

    s_0 = 1 # none broken
    s_1 = 2 # one broken
    s_2 = 3 # two broken

    a_m = 1 # manufacture
    a_e = 2 # manufacture + examine
    a_i = 3 # interrupt the line, inspect components, replace failed components
    a_r = 4 # interrupt the line, replace both components

    o_n = 1 # nondefective
    o_d = 2 # defective

    T[s_0, a_m, :] = [0.81, 0.18, 0.01] # 10% independent chance of part breaking
    T[s_0, a_e, :] = [0.81, 0.18, 0.01]
    T[s_0, a_i, :] = [1.00, 0.00, 0.00]
    T[s_0, a_r, :] = [1.00, 0.00, 0.00]
    T[s_1, a_m, :] = [0.00, 0.90, 0.10] # 10% chance of remaining part breaking
    T[s_1, a_e, :] = [0.00, 0.90, 0.10]
    T[s_1, a_i, :] = [1.00, 0.00, 0.00]
    T[s_1, a_r, :] = [1.00, 0.00, 0.00]
    T[s_2, a_m, :] = [0.00, 0.00, 1.00] # stay broken
    T[s_2, a_e, :] = [0.00, 0.00, 1.00]
    T[s_2, a_i, :] = [1.00, 0.00, 0.00]
    T[s_2, a_r, :] = [1.00, 0.00, 0.00]

    # There is a profit of 1 for producing a nondefective product.
    # Thus, the expected profit for beginning with 0, 1, or 2 defective parts is
    # 0.9025, 0.475, and 0.25, respectively.
    # Examining the finished product costs 0.25.
    # The inspect action incurs a 0.5 penalty plus replacement cost for each unit of 1.
    # The straight-up replacement action has no inspection cost but does incur a 2 unit cost.
    r_examine = -0.25
    r_inspect = -0.5
    r_replace = -2.0

    R[s_0, a_m] = 0.9025
    R[s_1, a_m] = 0.475
    R[s_2, a_m] = 0.25
    R[s_0, a_e] = 0.9025 + r_examine
    R[s_1, a_e] = 0.475 + r_examine
    R[s_2, a_e] = 0.25 + r_examine
    R[s_0, a_i] = r_inspect
    R[s_1, a_i] = r_inspect - 1.0 # replace 1
    R[s_2, a_i] = r_inspect - 2.0 # replace 2
    R[s_0, a_r] = r_replace
    R[s_1, a_r] = r_replace
    R[s_2, a_r] = r_replace

    # Probabilities of observing a nondefective product are 1.0, 0.5, and 0.25 if
    # there are 0, 1, or 2 faulty internal components.
    # If we don't examine, we always observe nondefective.
    O[a_m, s_0, :] = [1.00, 0.00]
    O[a_m, s_1, :] = [1.00, 0.00]
    O[a_m, s_2, :] = [1.00, 0.00]
    O[a_e, s_0, :] = [1.00, 0.00]
    O[a_e, s_1, :] = [0.50, 0.50]
    O[a_e, s_2, :] = [0.25, 0.75]
    O[a_i, s_0, :] = [1.00, 0.00]
    O[a_i, s_1, :] = [1.00, 0.00]
    O[a_i, s_2, :] = [1.00, 0.00]
    O[a_r, s_0, :] = [1.00, 0.00]
    O[a_r, s_1, :] = [1.00, 0.00]
    O[a_r, s_2, :] = [1.00, 0.00]

    return POMDP(DiscretePOMDP(T, R, O, γ))
end

# MachineReplacement = generate_machine_replacement_pomdp(1.0)

MACHINE_REPLACEMENT_ACTION_COLORS = Dict(
    1 => "pastelRed",
    2 => "pastelGreen",
    3 => "pastelBlue",
    4 => "pastelPurple",
)
MACHINE_REPLACEMENT_ACTION_NAMES = Dict(
    1 => "manufacture",
    2 => "examine",
    3 => "interrupt",
    4 => "replace",
)


text:
# Fast mode
fast_mode(X::Vector{T}) where T <: Integer = fast_mode(X, minimum(X):maximum(X))
function fast_mode(X, range)
    m, cm = (0, 0)
    for i in range
       c = count(isequal(i), X)
       if c > cm
        cm = c
        m = i
        end
    end
    return m
end

# Summary stats
function average(V::Vector{A}) where A <: AbstractArray
    N = length(V)
    M = V[1]
    for i in 2:N
       M .+= V[i]
    end
    M ./= N
    return M
end
average(V::Vector{T}) where T <: AbstractFloat = mean(V)
average(V::Vector{T}) where T <: Integer = fast_mode(V)

function average(SV::StructVector{TV}) where TV
    return TV(average.(values(StructArrays.components(SV))))
end

# Posterior expectation
function expectation(P::Posterior)
    return (hyperparameters = average(P.hyperparameters),
            latent = average(P.latent))
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
module Lifting
using DelimitedFiles
using DataFrames
using Dates
using Dierckx
using Plots

include("LiftingTypes.jl")
export AbstractProgression,
    AbstractProgramme,
    LinearProgression,
    DoubleProgression,
    PeriodProgression,
    BlockProgression
export calcIntensity,
    calcRPE,
    calcReps,
    calcIntensityRatio,
    calcRepRatio,
    calcRPERatio,
    calcRepMax,
    adjustMaxes,
    updateMaxes,
    adjustMaxes!,
    updateMaxes!,
    calcTrainingMaxLogs,
    intensityArb,
    SetScheme,
    Progression,
    Exercise,
    calcWeights,
    Programme,
    makeDays,
    numDays
include("LiftingIO.jl")
export numDays, loadLogFile, loadLogFile, plotData, plotData!, scatterData, scatterData!

include("./assets/LiftingAssets.jl")
# Default dictionaries, uncomment to auto export.
# export Lifting_Aux
# export Lifting_Exercise_Names
# export Lifting_Progressions
# export Lifting_Programmes

end # module

max_stars_repo_name: dcelisgarza/Lifting


text:
@testset "cxsc.intervaladdsub" begin
    @test +(interval(10.0, 20.0), interval(13.0, 17.0)) == interval(23.0, 37.0)
    @test +(interval(13.0, 17.0), interval(10.0, 20.0)) == interval(23.0, 37.0)
    @test -(interval(10.0, 20.0), interval(13.0, 16.0)) == interval(-6.0, 7.0)
    @test -(interval(13.0, 16.0), interval(10.0, 20.0)) == interval(-7.0, 6.0)
    @test -(interval(10.0, 20.0)) == interval(-20.0, -10.0)
    @test +(interval(10.0, 20.0)) == interval(10.0, 20.0)
end

@testset "cxsc.intervalmuldiv" begin
    @test *(interval(1.0, 2.0), interval(3.0, 4.0)) == interval(3.0, 8.0)
    @test *(interval(-1.0, 2.0), interval(3.0, 4.0)) == interval(-4.0, 8.0)
    @test *(interval(-2.0, 1.0), interval(3.0, 4.0)) == interval(-8.0, 4.0)
    @test *(interval(-2.0, -1.0), interval(3.0, 4.0)) == interval(-8.0, -3.0)
    @test *(interval(1.0, 2.0), interval(-3.0, 4.0)) == interval(-6.0, 8.0)
    @test *(interval(-1.0, 2.0), interval(-3.0, 4.0)) == interval(-6.0, 8.0)
    @test *(interval(-2.0, 1.0), interval(-3.0, 4.0)) == interval(-8.0, 6.0)
    @test *(interval(-2.0, -1.0), interval(-3.0, 4.0)) == interval(-8.0, 6.0)
    @test *(interval(1.0, 2.0), interval(-4.0, 3.0)) == interval(-8.0, 6.0)
    @test *(interval(-1.0, 2.0), interval(-4.0, 3.0)) == interval(-8.0, 6.0)
    @test *(interval(-2.0, 1.0), interval(-4.0, 3.0)) == interval(-6.0, 8.0)
    @test *(interval(-2.0, -1.0), interval(-4.0, 3.0)) == interval(-6.0, 8.0)
    @test *(interval(1.0, 2.0), interval(-4.0, -3.0)) == interval(-8.0, -3.0)
    @test *(interval(-1.0, 2.0), interval(-4.0, -3.0)) == interval(-8.0, 4.0)
    @test *(interval(-2.0, -1.0), interval(-4.0, -3.0)) == interval(3.0, 8.0)
    @test /(interval(1.0, 2.0), interval(4.0, 8.0)) == interval(0.125, 0.5)
    @test /(interval(-1.0, 2.0), interval(4.0, 8.0)) == interval(-0.25, 0.5)
    @test /(interval(-2.0, 1.0), interval(4.0, 8.0)) == interval(-0.5, 0.25)
    @test /(interval(-2.0, -1.0), interval(4.0, 8.0)) == interval(-0.5, -0.125)
    @test /(interval(1.0, 2.0), interval(-4.0, 8.0)) == entireinterval()
    @test /(interval(-1.0, 2.0), interval(-4.0, 8.0)) == entireinterval()
    @test /(interval(-2.0, 1.0), interval(-4.0, 8.0)) == entireinterval()
    @test /(interval(-2.0, -1.0), interval(-4.0, 8.0)) == entireinterval()
    @test /(interval(1.0, 2.0), interval(-8.0, 4.0)) == entireinterval()
    @test /(interval(-1.0, 2.0), interval(-8.0, 4.0)) == entireinterval()
    @test /(interval(-2.0, 1.0), interval(-8.0, 4.0)) == entireinterval()
    @test /(interval(-2.0, -1.0), interval(-8.0, 4.0)) == entireinterval()
    @test /(interval(1.0, 2.0), interval(-8.0, -4.0)) == interval(-0.5, -0.125)
    @test /(interval(-1.0, 2.0), interval(-8.0, -4.0)) == interval(-0.5, 0.25)
    @test /(interval(-2.0, 1.0), interval(-8.0, -4.0)) == interval(-0.25, 0.5)
    @test /(interval(-2.0, -1.0), interval(-8.0, -4.0)) == interval(0.125, 0.5)
end

@testset "cxsc.intervalsetops" begin
    @test hull(interval(-2.0, 2.0), interval(-4.0, -3.0)) == interval(-4.0, 2.0)
    @test hull(interval(-2.0, 2.0), interval(-4.0, -1.0)) == interval(-4.0, 2.0)
    @test hull(interval(-2.0, 2.0), interval(-4.0, 4.0)) == interval(-4.0, 4.0)
    @test hull(interval(-2.0, 2.0), interval(-1.0, 1.0)) == interval(-2.0, 2.0)
    @test hull(interval(-2.0, 2.0), interval(1.0, 4.0)) == interval(-2.0, 4.0)
    @test hull(interval(-2.0, 2.0), interval(3.0, 4.0)) == interval(-2.0, 4.0)
    @test hull(interval(-4.0, -3.0), interval(-2.0, 2.0)) == interval(-4.0, 2.0)
    @test hull(interval(-4.0, -1.0), interval(-2.0, 2.0)) == interval(-4.0, 2.0)
    @test hull(interval(-4.0, 4.0), interval(-2.0, 2.0)) == interval(-4.0, 4.0)
    @test hull(interval(-1.0, 1.0), interval(-2.0, 2.0)) == interval(-2.0, 2.0)
    @test hull(interval(1.0, 4.0), interval(-2.0, 2.0)) == interval(-2.0, 4.0)
    @test hull(interval(3.0, 4.0), interval(-2.0, 2.0)) == interval(-2.0, 4.0)
    @test intersect(interval(-2.0, 2.0), interval(-4.0, -3.0)) == emptyinterval()
    @test intersect(interval(-2.0, 2.0), interval(-4.0, -1.0)) == interval(-2.0, -1.0)
    @test intersect(interval(-2.0, 2.0), interval(-4.0, 4.0)) == interval(-2.0, 2.0)
    @test intersect(interval(-2.0, 2.0), interval(-1.0, 1.0)) == interval(-1.0, 1.0)
    @test intersect(interval(-2.0, 2.0), interval(1.0, 4.0)) == interval(1.0, 2.0)
    @test intersect(interval(-2.0, 2.0), interval(3.0, 4.0)) == emptyinterval()
    @test intersect(interval(-4.0, -3.0), interval(-2.0, 2.0)) == emptyinterval()
    @test intersect(interval(-4.0, -1.0), interval(-2.0, 2.0)) == interval(-2.0, -1.0)
    @test intersect(interval(-4.0, 4.0), interval(-2.0, 2.0)) == interval(-2.0, 2.0)
    @test intersect(interval(-1.0, 1.0), interval(-2.0, 2.0)) == interval(-1.0, 1.0)
    @test intersect(interval(1.0, 4.0), interval(-2.0, 2.0)) == interval(1.0, 2.0)
    @test intersect(interval(3.0, 4.0), interval(-2.0, 2.0)) == emptyinterval()
end

@testset "cxsc.intervalmixsetops" begin
    @test hull(interval(-2.0, 2.0), interval(-4.0, -4.0)) == interval(-4.0, 2.0)
    @test hull(interval(-2.0, 2.0), interval(1.0, 1.0)) == interval(-2.0, 2.0)
    @test hull(interval(-2.0, 2.0), interval(4.0, 4.0)) == interval(-2.0, 4.0)
    @test hull(interval(-4.0, -4.0), interval(-2.0, 2.0)) == interval(-4.0, 2.0)
    @test hull(interval(1.0, 1.0), interval(-2.0, 2.0)) == interval(-2.0, 2.0)
    @test hull(interval(4.0, 4.0), interval(-2.0, 2.0)) == interval(-2.0, 4.0)
    @test intersect(interval(-2.0, 2.0), interval(-4.0, -4.0)) == emptyinterval()
    @test intersect(interval(-2.0, 2.0), interval(1.0, 1.0)) == interval(1.0, 1.0)
    @test intersect(interval(-2.0, 2.0), interval(4.0, 4.0)) == emptyinterval()
    @test intersect(interval(-4.0, -4.0), interval(-2.0, 2.0)) == emptyinterval()
    @test intersect(interval(1.0, 1.0), interval(-2.0, 2.0)) == interval(1.0, 1.0)
    @test intersect(interval(4.0, 4.0), interval(-2.0, 2.0)) == emptyinterval()
end

@testset "cxsc.scalarmixsetops" begin
    @test hull(interval(-2.0, -2.0), interval(-4.0, -4.0)) == interval(-4.0, -2.0)
    @test hull(interval(-2.0, -2.0), interval(-2.0, -2.0)) == interval(-2.0, -2.0)
    @test hull(interval(-2.0, -2.0), interval(2.0, 2.0)) == interval(-2.0, 2.0)
    @test hull(interval(-4.0, -4.0), interval(-2.0, -2.0)) == interval(-4.0, -2.0)
    @test hull(interval(-2.0, -2.0), interval(-2.0, -2.0)) == interval(-2.0, -2.0)
    @test hull(interval(2.0, 2.0), interval(-2.0, -2.0)) == interval(-2.0, 2.0)
end

@testset "cxsc.intervalsetcompops" begin
    @test !isinterior(interval(-1.0, 2.0), interval(-1.0, 2.0))
    @test isinterior(interval(-2.0, 1.0), interval(-3.0, 2.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-1.0, 1.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-1.0, 2.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-2.0, 1.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-2.0, 3.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-3.0, 2.0))
    @test !isinterior(interval(-1.0, 2.0), interval(-1.0, 2.0))
    @test !isinterior(interval(-3.0, 2.0), interval(-2.0, 1.0))
    @test isinterior(interval(-1.0, 1.0), interval(-2.0, 2.0))
    @test !isinterior(interval(-1.0, 2.0), interval(-2.0, 2.0))
    @test !isinterior(interval(-2.0, 1.0), interval(-2.0, 2.0))
    @test !isinterior(interval(-2.0, 3.0), interval(-2.0, 2.0))
    @test !isinterior(interval(-3.0, 2.0), interval(-2.0, 2.0))
    @test issubset(interval(-1.0, 2.0), interval(-1.0, 2.0))
    @test issubset(interval(-2.0, 1.0), interval(-3.0, 2.0))
    @test !issubset(interval(-2.0, 2.0), interval(-1.0, 1.0))
    @test !issubset(interval(-2.0, 2.0), interval(-1.0, 2.0))
    @test !issubset(interval(-2.0, 2.0), interval(-2.0, 1.0))
    @test issubset(interval(-2.0, 2.0), interval(-2.0, 3.0))
    @test issubset(interval(-2.0, 2.0), interval(-3.0, 2.0))
    @test !issubset(interval(-3.0, 2.0), interval(-2.0, 1.0))
    @test issubset(interval(-1.0, 1.0), interval(-2.0, 2.0))
    @test issubset(interval(-1.0, 2.0), interval(-2.0, 2.0))
    @test issubset(interval(-2.0, 1.0), interval(-2.0, 2.0))
    @test !issubset(interval(-2.0, 3.0), interval(-2.0, 2.0))
    @test !issubset(interval(-3.0, 2.0), interval(-2.0, 2.0))
    @test ==(interval(-1.0, 2.0), interval(-1.0, 2.0))
    @test !==(interval(-2.0, 1.0), interval(-3.0, 2.0))
    @test !==(interval(-2.0, 2.0), interval(-1.0, 1.0))
    @test !==(interval(-2.0, 2.0), interval(-1.0, 2.0))
    @test !==(interval(-2.0, 2.0), interval(-2.0, 1.0))
    @test !==(interval(-2.0, 2.0), interval(-2.0, 3.0))
    @test !==(interval(-2.0, 2.0), interval(-3.0, 2.0))
end

@testset "cxsc.intervalscalarsetcompops" begin
    @test !isinterior(interval(-1.0, 2.0), interval(-2.0, -2.0))
    @test !isinterior(interval(-2.0, 2.0), interval(-2.0, -2.0))
    @test !isinterior(interval(-2.0, 2.0), interval(0.0, 0.0))
    @test !isinterior(interval(-2.0, 2.0), interval(2.0, 2.0))
    @test !isinterior(interval(-2.0, 2.0), interval(3.0, 3.0))
    @test !isinterior(interval(-1.0, -1.0), interval(1.0, 1.0))
    @test !isinterior(interval(-1.0, -1.0), interval(-1.0, -1.0))
    @test !isinterior(interval(-2.0, -2.0), interval(-1.0, 2.0))
    @test !isinterior(interval(-2.0, -2.0), interval(-2.0, 2.0))
    @test isinterior(interval(0.0, 0.0), interval(-2.0, 2.0))
    @test !isinterior(interval(2.0, 2.0), interval(-2.0, 2.0))
    @test !isinterior(interval(3.0, 3.0), interval(-2.0, 2.0))
    @test !isinterior(interval(1.0, 1.0), interval(-1.0, -1.0))
    @test !isinterior(interval(-1.0, -1.0), interval(-1.0, -1.0))
    @test !issubset(interval(-1.0, 2.0), interval(-2.0, -2.0))
    @test !issubset(interval(-2.0, 2.0), interval(-2.0, -2.0))
    @test !issubset(interval(-2.0, 2.0), interval(0.0, 0.0))
    @test !issubset(interval(-2.0, 2.0), interval(2.0, 2.0))
    @test !issubset(interval(-2.0, 2.0), interval(3.0, 3.0))
    @test !issubset(interval(-1.0, -1.0), interval(1.0, 1.0))
    @test issubset(interval(-1.0, -1.0), interval(-1.0, -1.0))
    @test !issubset(interval(-2.0, -2.0), interval(-1.0, 2.0))
    @test issubset(interval(-2.0, -2.0), interval(-2.0, 2.0))
    @test issubset(interval(0.0, 0.0), interval(-2.0, 2.0))
    @test issubset(interval(2.0, 2.0), interval(-2.0, 2.0))
    @test !issubset(interval(3.0, 3.0), interval(-2.0, 2.0))
    @test !issubset(interval(1.0, 1.0), interval(-1.0, -1.0))
    @test issubset(interval(-1.0, -1.0), interval(-1.0, -1.0))
    @test !==(interval(-1.0, 2.0), interval(-2.0, -2.0))
    @test !==(interval(-2.0, 2.0), interval(-2.0, -2.0))
    @test !==(interval(-2.0, 2.0), interval(0.0, 0.0))
    @test !==(interval(-2.0, 2.0), interval(2.0, 2.0))
    @test !==(interval(-2.0, 2.0), interval(3.0, 3.0))
    @test !==(interval(-1.0, -1.0), interval(1.0, 1.0))
    @test ==(interval(-1.0, -1.0), interval(-1.0, -1.0))
end

@testset "cxsc.intervalstdfunc" begin
    @test abs2(interval(11.0, 11.0)) == interval(121.0, 121.0)
    @test abs2(interval(0.0, 0.0)) == interval(0.0, 0.0)
    @test abs2(interval(-9.0, -9.0)) == interval(81.0, 81.0)
    @test sqrt(interval(121.0, 121.0)) == interval(11.0, 11.0)
    @test sqrt(interval(0.0, 0.0)) == interval(0.0, 0.0)
    @test sqrt(interval(81.0, 81.0)) == interval(9.0, 9.0)
    # rootn [27.0, 27.0] 3 = [3.0, 3.0]
    # rootn [0.0, 0.0] 4 = [0.0, 0.0]
    # rootn [1024.0, 1024.0] 10 = [2.0, 2.0]
    @test ^(interval(2.0, 2.0), interval(2.0, 2.0)) ⊇ interval(4.0, 4.0)
    @test ^(interval(4.0, 4.0), interval(5.0, 5.0)) ⊇ interval(1024.0, 1024.0)
    @test ^(interval(2.0, 2.0), interval(3.0, 3.0)) ⊇ interval(8.0, 8.0)
end


text:
@doc """
`path = finddir(depname)` find the (versioned) path for a particular
dependency specified by `depname`.  `depname` might be a string, e.g.,
`arpack-ng`, or a tuple representing a subdirectory, e.g.,
`("SuiteSparse","CHOLMOD")`. On output, `path` is the full pathname to
the most recent version of the dependency.
""" ->
finddir(depsdir, depname::Tuple) = joinpath(finddir(depsdir, depname[1]), depname[2:end]...)

function finddir(depsdir, topdir::AbstractString)
    depsnames = readdir(depsdir)
    names = filter(x->isdir(joinpath(depsdir, x)), filter(x->startswith(x, topdir), depsnames))
    versions = map(x->Base.VersionNumber(x[length(topdir)+2:end]), names)
    ind = indmax(versions)
    names[ind]
end


text:
# explore warping functions.

using FFTW
import PyPlot
import BSON
import Optim
import Random
using LinearAlgebra

import Interpolations

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

Random.seed!(25)

a = [0.5; 5.0; 1/5] # higher makes sharper transition.
b = [1.0; 1.0; 1.0] .*1

c = [1.0; 1.0; 1.0] # higher makes sharper transition.
d = [0.0; 0.0; 0.0]

fs = collect( xx->( evalf2(evalf1(xx, a[i], b[i]), c[i], c[i]) ) for i = 1:length(a) )
#fs = collect( xx->( evalf1(0.5*evalf2(xx, a[i], b[i])+0.5, c[i], c[i]) ) for i = 1:length(a) )

function evalf1(x::T, a, b)::T where T
    return a*log(x/(1-x))+b
end

## algebraic sigmoid function. f: ℝ → [-1, 1]
# function evalf2(x::T, a, b)::T where T
#     return (x-b)/sqrt(a+(x-b)^2)
# end

## logistic function. f: ℝ → [-1, 1]
function evalf2(x::T, a, b)::T where T
    return 1/(1+exp(-a*(x-b)))
end

x = LinRange(0 + 1e-2, 1 - 1e-2, 500)
#x = LinRange(-30 , 30, 500)

PyPlot.figure(fig_num)
fig_num += 1

for i = 1:length(fs)
    PyPlot.plot(x, fs[i].(x), label = "f[$(i)]")
end

PyPlot.legend()
PyPlot.xlabel("x")
PyPlot.ylabel("f")
PyPlot.title("target warp func")

##### streamline evalf2 ∘ evalf1.

include("../src/MonotoneMaps.jl")
import .MonotoneMaps

q = tt->MonotoneMaps.evalcompositelogisticprobit(tt, 0.25, 0.83, -1.0, 1.0)

x = LinRange(-1 + 1e-2, 1 - 1e-2, 500)

PyPlot.figure(fig_num)
fig_num += 1

ax = PyPlot.axes()
ax[:set_ylim]([-1,1])
PyPlot.plot(x, q.(x), label = "q")

PyPlot.legend()
PyPlot.xlabel("x")
PyPlot.ylabel("q")
PyPlot.title("evalcompositelogisticprobit()")

## inverse.

q_inv = yy->MonotoneMaps.evalinversecompositelogisticprobit(yy, 0.25, 0.83, -1.0, 1.0)

x = LinRange(-1 + 1e-2, 1 - 1e-2, 500)

PyPlot.figure(fig_num)
fig_num += 1

ax = PyPlot.axes()
ax[:set_ylim]([-1,1])
PyPlot.plot(x, q_inv.(x), label = "q_inv")

PyPlot.legend()
PyPlot.xlabel("x")
PyPlot.ylabel("q")
PyPlot.title("evalinversecompositelogisticprobit()")

q_inv(q(-0.94)) # check.


text:
"""
    VectorSpaceType

Abstract type for tangent spaces, cotangent spaces, their tensor products,
exterior products, etc.

Every vector space `fiber` is supposed to provide:
* a method of constructing vectors,
* basic operations: addition, subtraction, multiplication by a scalar
  and negation (unary minus),
* `zero_vector(fiber, p)` to construct zero vectors at point `p`,
* `allocate(X)` and `allocate(X, T)` for vector `X` and type `T`,
* `copyto!(X, Y)` for vectors `X` and `Y`,
* `number_eltype(v)` for vector `v`,
* [`vector_space_dimension`](@ref).

Optionally:
* inner product via `inner` (used to provide Riemannian metric on vector
  bundles),
* [`flat`](https://juliamanifolds.github.io/Manifolds.jl/stable/features/atlases.html#Manifolds.flat-Tuple{AbstractManifold,%20Any,%20Any}) and [`sharp`](https://juliamanifolds.github.io/Manifolds.jl/stable/features/atlases.html#Manifolds.sharp-Tuple{AbstractManifold,%20Any,%20Any}),
* `norm` (by default uses `inner`),
* [`project`](@ref) (for embedded vector spaces),
* [`representation_size`](@ref),
* broadcasting for basic operations.
"""
abstract type VectorSpaceType end

struct TangentSpaceType <: VectorSpaceType end

struct CotangentSpaceType <: VectorSpaceType end

TCoTSpaceType = Union{TangentSpaceType,CotangentSpaceType}

const TangentSpace = TangentSpaceType()
const CotangentSpace = CotangentSpaceType()

"""
    AbstractBasis{𝔽,VST<:VectorSpaceType}

Abstract type that represents a basis of vector space of type `VST` on a manifold or
a subset of it.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
abstract type AbstractBasis{𝔽,VST<:VectorSpaceType} end

"""
    DefaultBasis{𝔽,VST<:VectorSpaceType}

An arbitrary basis of vector space of type `VST` on a manifold. This will usually
be the fastest basis available for a manifold.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
struct DefaultBasis{𝔽,VST<:VectorSpaceType} <: AbstractBasis{𝔽,VST}
    vector_space::VST
end
function DefaultBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace)
    return DefaultBasis{𝔽,typeof(vs)}(vs)
end
function DefaultBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽}
    return DefaultBasis{𝔽,typeof(vs)}(vs)
end
function DefaultBasis{𝔽,TangentSpaceType}() where {𝔽}
    return DefaultBasis{𝔽,TangentSpaceType}(TangentSpace)
end

"""
    AbstractOrthogonalBasis{𝔽,VST<:VectorSpaceType}

Abstract type that represents an orthonormal basis of vector space of type `VST` on a
manifold or a subset of it.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
abstract type AbstractOrthogonalBasis{𝔽,VST<:VectorSpaceType} <: AbstractBasis{𝔽,VST} end

"""
    DefaultOrthogonalBasis{𝔽,VST<:VectorSpaceType}

An arbitrary orthogonal basis of vector space of type `VST` on a manifold. This will usually
be the fastest orthogonal basis available for a manifold.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
struct DefaultOrthogonalBasis{𝔽,VST<:VectorSpaceType} <: AbstractOrthogonalBasis{𝔽,VST}
    vector_space::VST
end
function DefaultOrthogonalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace)
    return DefaultOrthogonalBasis{𝔽,typeof(vs)}(vs)
end
function DefaultOrthogonalBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽}
    return DefaultOrthogonalBasis{𝔽,typeof(vs)}(vs)
end
function DefaultOrthogonalBasis{𝔽,TangentSpaceType}() where {𝔽}
    return DefaultOrthogonalBasis{𝔽,TangentSpaceType}(TangentSpace)
end

struct VeeOrthogonalBasis{𝔽} <: AbstractOrthogonalBasis{𝔽,TangentSpaceType} end
VeeOrthogonalBasis(𝔽::AbstractNumbers = ℝ) = VeeOrthogonalBasis{𝔽}()

"""
    AbstractOrthonormalBasis{𝔽,VST<:VectorSpaceType}

Abstract type that represents an orthonormal basis of vector space of type `VST` on a
manifold or a subset of it.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
abstract type AbstractOrthonormalBasis{𝔽,VST<:VectorSpaceType} <:
              AbstractOrthogonalBasis{𝔽,VST} end

"""
    DefaultOrthonormalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace)

An arbitrary orthonormal basis of vector space of type `VST` on a manifold. This will usually
be the fastest orthonormal basis available for a manifold.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# See also

[`VectorSpaceType`](@ref)
"""
struct DefaultOrthonormalBasis{𝔽,VST<:VectorSpaceType} <: AbstractOrthonormalBasis{𝔽,VST}
    vector_space::VST
end

function DefaultOrthonormalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace)
    return DefaultOrthonormalBasis{𝔽,typeof(vs)}(vs)
end
function DefaultOrthonormalBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽}
    return DefaultOrthonormalBasis{𝔽,typeof(vs)}(vs)
end
function DefaultOrthonormalBasis{𝔽,TangentSpaceType}() where {𝔽}
    return DefaultOrthonormalBasis{𝔽,TangentSpaceType}(TangentSpace)
end

"""
    ProjectedOrthonormalBasis(method::Symbol, 𝔽::AbstractNumbers = ℝ)

An orthonormal basis that comes from orthonormalization of basis vectors
of the ambient space projected onto the subspace representing the tangent space
at a given point.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

Available methods:
  - `:gram_schmidt` uses a modified Gram-Schmidt orthonormalization.
  - `:svd` uses SVD decomposition to orthogonalize projected vectors.
    The SVD-based method should be more numerically stable at the cost of
    an additional assumption (local metric tensor at a point where the
    basis is calculated has to be diagonal).
"""
struct ProjectedOrthonormalBasis{Method,𝔽} <: AbstractOrthonormalBasis{𝔽,TangentSpaceType} end

function ProjectedOrthonormalBasis(method::Symbol, 𝔽::AbstractNumbers = ℝ)
    return ProjectedOrthonormalBasis{method,𝔽}()
end

@doc raw"""
    GramSchmidtOrthonormalBasis{𝔽} <: AbstractOrthonormalBasis{𝔽}

An orthonormal basis obtained from a basis.

# Constructor
    GramSchmidtOrthonormalBasis(𝔽::AbstractNumbers = ℝ)
"""
struct GramSchmidtOrthonormalBasis{𝔽} <: AbstractOrthonormalBasis{𝔽,TangentSpaceType} end
GramSchmidtOrthonormalBasis(𝔽::AbstractNumbers = ℝ) = GramSchmidtOrthonormalBasis{𝔽}()

@doc raw"""
    DiagonalizingOrthonormalBasis{𝔽,TV} <: AbstractOrthonormalBasis{𝔽,TangentSpaceType}

An orthonormal basis `Ξ` as a vector of tangent vectors (of length determined by
[`manifold_dimension`](@ref)) in the tangent space that diagonalizes the curvature
tensor $R(u,v)w$ and where the direction `frame_direction` $v$ has curvature `0`.

The type parameter `𝔽` denotes the [`AbstractNumbers`](@ref) that will be used
for the vectors elements.

# Constructor

    DiagonalizingOrthonormalBasis(frame_direction, 𝔽::AbstractNumbers = ℝ)
"""
struct DiagonalizingOrthonormalBasis{𝔽,TV} <: AbstractOrthonormalBasis{𝔽,TangentSpaceType}
    frame_direction::TV
end
function DiagonalizingOrthonormalBasis(X, 𝔽::AbstractNumbers = ℝ)
    return DiagonalizingOrthonormalBasis{𝔽,typeof(X)}(X)
end
struct DiagonalizingBasisData{D,V,ET}
    frame_direction::D
    eigenvalues::ET
    vectors::V
end

const DefaultOrDiagonalizingBasis{𝔽} =
    Union{DefaultOrthonormalBasis{𝔽,TangentSpaceType},DiagonalizingOrthonormalBasis{𝔽}}

"""
    CachedBasis{𝔽,V,<:AbstractBasis{𝔽}} <: AbstractBasis{𝔽}

A cached version of the given `basis` with precomputed basis vectors. The basis vectors
are stored in `data`, either explicitly (like in cached variants of
[`ProjectedOrthonormalBasis`](@ref)) or implicitly.

# Constructor

    CachedBasis(basis::AbstractBasis, data)
"""
struct CachedBasis{𝔽,B,V} <:
       AbstractBasis{𝔽,TangentSpaceType} where {B<:AbstractBasis{𝔽,TangentSpaceType},V}
    data::V
end
function CachedBasis(::B, data::V) where {V,𝔽,B<:AbstractBasis{𝔽,<:TangentSpaceType}}
    return CachedBasis{𝔽,B,V}(data)
end
function CachedBasis(basis::CachedBasis) # avoid double encapsulation
    return basis
end
function CachedBasis(
    basis::DiagonalizingOrthonormalBasis,
    eigenvalues::ET,
    vectors::T,
) where {ET<:AbstractVector,T<:AbstractVector}
    data = DiagonalizingBasisData(basis.frame_direction, eigenvalues, vectors)
    return CachedBasis(basis, data)
end

# forward declarations
function get_coordinates end
function get_vector end

const all_uncached_bases{T} = Union{
    AbstractBasis{<:Any,T},
    DefaultBasis{<:Any,T},
    DefaultOrthogonalBasis{<:Any,T},
    DefaultOrthonormalBasis{<:Any,T},
}

"""
    allocate_coordinates(M::AbstractManifold, p, T, n::Int)

Allocate vector of coordinates of length `n` of type `T` of a vector at point `p`
on manifold `M`.
"""
allocate_coordinates(::AbstractManifold, p, T, n::Int) = allocate(p, T, n)
function allocate_coordinates(M::AbstractManifold, p::Int, T, n::Int)
    return (representation_size(M) == () && n == 0) ? zero(T) : zeros(T, n)
end
function allocate_result(
    M::AbstractManifold,
    f::typeof(get_coordinates),
    p,
    X,
    basis::AbstractBasis,
)
    T = allocate_result_type(M, f, (p, X))
    return allocate_coordinates(M, p, T, number_of_coordinates(M, basis))
end

@inline function allocate_result_type(
    M::AbstractManifold,
    f::Union{typeof(get_coordinates),typeof(get_vector)},
    args::Tuple{Any,Vararg{Any}},
)
    apf = allocation_promotion_function(M, f, args)
    return apf(
        invoke(allocate_result_type, Tuple{AbstractManifold,Any,typeof(args)}, M, f, args),
    )
end

"""
    allocation_promotion_function(M::AbstractManifold, f, args::Tuple)

Determine the function that must be used to ensure that the allocated representation is of
the right type. This is needed for [`get_vector`](@ref) when a point on a complex manifold
is represented by a real-valued vectors with a real-coefficient basis, so that
a complex-valued vector representation is allocated.
"""
allocation_promotion_function(M::AbstractManifold, f, args::Tuple) = identity

function combine_allocation_promotion_functions(f::T, ::T) where {T}
    return f
end
function combine_allocation_promotion_functions(::typeof(complex), ::typeof(identity))
    return complex
end
function combine_allocation_promotion_functions(::typeof(identity), ::typeof(complex))
    return complex
end

@doc raw"""
    dual_basis(M::AbstractManifold, p, B::AbstractBasis)

Get the dual basis to `B`, a basis of a vector space at point `p` from manifold `M`.

The dual to the $i$th vector $v_i$ from basis `B` is a vector $v^i$ from the dual space
such that $v^i(v_j) = δ^i_j$, where $δ^i_j$ is the Kronecker delta symbol:
````math
δ^i_j = \begin{cases}
1 & \text{ if } i=j, \\
0 & \text{ otherwise.}
\end{cases}
````
"""
dual_basis(M::AbstractManifold, p, B::AbstractBasis) = _dual_basis(M, p, B)

function _dual_basis(
    ::AbstractManifold,
    p,
    ::DefaultOrthonormalBasis{𝔽,TangentSpaceType},
) where {𝔽}
    return DefaultOrthonormalBasis{𝔽}(CotangentSpace)
end
function _dual_basis(
    ::AbstractManifold,
    p,
    ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    return DefaultOrthonormalBasis{𝔽}(TangentSpace)
end

function _euclidean_basis_vector(p::StridedArray, i)
    X = zero(p)
    X[i] = 1
    return X
end
function _euclidean_basis_vector(p, i)
    # when p is for example a SArray
    X = similar(p)
    copyto!(X, zero(p))
    X[i] = 1
    return X
end

"""
    get_basis(M::AbstractManifold, p, B::AbstractBasis; kwargs...) -> CachedBasis

Compute the basis vectors of the tangent space at a point on manifold `M`
represented by `p`.

Returned object derives from [`AbstractBasis`](@ref) and may have a field `.vectors`
that stores tangent vectors or it may store them implicitly, in which case
the function [`get_vectors`](@ref) needs to be used to retrieve the basis vectors.

See also: [`get_coordinates`](@ref), [`get_vector`](@ref)
"""
function get_basis(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    return _get_basis(M, p, B; kwargs...)
end

function _get_basis(::AbstractManifold, ::Any, B::CachedBasis)
    return B
end
function _get_basis(M::AbstractManifold, p, B::ProjectedOrthonormalBasis{:svd,ℝ})
    S = representation_size(M)
    PS = prod(S)
    dim = manifold_dimension(M)
    # projection
    # TODO: find a better way to obtain a basis of the ambient space
    Xs = [
        convert(Vector, reshape(project(M, p, _euclidean_basis_vector(p, i)), PS)) for
        i in eachindex(p)
    ]
    O = reduce(hcat, Xs)
    # orthogonalization
    # TODO: try using rank-revealing QR here
    decomp = svd(O)
    rotated = Diagonal(decomp.S) * decomp.Vt
    vecs = [collect(reshape(rotated[i, :], S)) for i in 1:dim]
    # normalization
    for i in 1:dim
        i_norm = norm(M, p, vecs[i])
        vecs[i] /= i_norm
    end
    return CachedBasis(B, vecs)
end
function _get_basis(
    M::AbstractManifold,
    p,
    B::ProjectedOrthonormalBasis{:gram_schmidt,ℝ};
    kwargs...,
)
    E = [project(M, p, _euclidean_basis_vector(p, i)) for i in eachindex(p)]
    V = gram_schmidt(M, p, E; kwargs...)
    return CachedBasis(B, V)
end
function _get_basis(M::AbstractManifold, p, B::VeeOrthogonalBasis)
    return get_basis_vee(M, p, number_system(B))
end
function get_basis_vee(M::AbstractManifold, p, N)
    return get_basis(M, p, DefaultOrthogonalBasis(N))
end

function _get_basis(M::AbstractManifold, p, B::DefaultBasis)
    return get_basis_default(M, p, number_system(B))
end
function get_basis_default(M::AbstractManifold, p, N)
    return get_basis(M, p, DefaultOrthogonalBasis(N))
end

function _get_basis(M::AbstractManifold, p, B::DefaultOrthogonalBasis)
    return get_basis_orthogonal(M, p, number_system(B))
end
function get_basis_orthogonal(M::AbstractManifold, p, N)
    return get_basis(M, p, DefaultOrthonormalBasis(N))
end

function _get_basis(M::AbstractManifold, p, B::DiagonalizingOrthonormalBasis)
    return get_basis_diagonalizing(M, p, B)
end
function get_basis_diagonalizing end

function _get_basis(M::AbstractManifold, p, B::DefaultOrthonormalBasis)
    return get_basis_orthonormal(M, p, number_system(B))
end
function get_basis_orthonormal(M::AbstractManifold, p, N::AbstractNumbers; kwargs...)
    B = DefaultOrthonormalBasis(N)
    dim = number_of_coordinates(M, B)
    return CachedBasis(
        B,
        [get_vector(M, p, [ifelse(i == j, 1, 0) for j in 1:dim], B) for i in 1:dim],
    )
end

@doc raw"""
    get_coordinates(M::AbstractManifold, p, X, B::AbstractBasis)
    get_coordinates(M::AbstractManifold, p, X, B::CachedBasis)

Compute a one-dimensional vector of coefficients of the tangent vector `X`
at point denoted by `p` on manifold `M` in basis `B`.

Depending on the basis, `p` may not directly represent a point on the manifold.
For example if a basis transported along a curve is used, `p` may be the coordinate
along the curve. If a [`CachedBasis`](@ref) is provided, their stored vectors are used,
otherwise the user has to provide a method to compute the coordinates.

For the [`CachedBasis`](@ref) keep in mind that the reconstruction with [`get_vector`](@ref)
requires either a dual basis or the cached basis to be selfdual, for example orthonormal

See also: [`get_vector`](@ref), [`get_basis`](@ref)
"""
function get_coordinates(M::AbstractManifold, p, X, B::AbstractBasis)
    return _get_coordinates(M, p, X, B)
end

function _get_coordinates(M::AbstractManifold, p, X, B::VeeOrthogonalBasis)
    return get_coordinates_vee(M, p, X, number_system(B))
end
function get_coordinates_vee(M::AbstractManifold, p, X, N)
    return get_coordinates(M, p, X, DefaultOrthogonalBasis(N))
end

function _get_coordinates(M::AbstractManifold, p, X, B::DefaultBasis)
    return get_coordinates_default(M, p, X, number_system(B))
end
function get_coordinates_default(M::AbstractManifold, p, X, N::AbstractNumbers)
    return get_coordinates(M, p, X, DefaultOrthogonalBasis(N))
end

function _get_coordinates(M::AbstractManifold, p, X, B::DefaultOrthogonalBasis)
    return get_coordinates_orthogonal(M, p, X, number_system(B))
end
function get_coordinates_orthogonal(M::AbstractManifold, p, X, N)
    return get_coordinates_orthonormal(M, p, X, N)
end

function _get_coordinates(M::AbstractManifold, p, X, B::DefaultOrthonormalBasis)
    return get_coordinates_orthonormal(M, p, X, number_system(B))
end
function get_coordinates_orthonormal(M::AbstractManifold, p, X, N)
    Y = allocate_result(M, get_coordinates, p, X, DefaultOrthonormalBasis(N))
    return get_coordinates_orthonormal!(M, Y, p, X, N)
end

function _get_coordinates(M::AbstractManifold, p, X, B::DiagonalizingOrthonormalBasis)
    return get_coordinates_diagonalizing(M, p, X, B)
end
function get_coordinates_diagonalizing(
    M::AbstractManifold,
    p,
    X,
    B::DiagonalizingOrthonormalBasis,
)
    Y = allocate_result(M, get_coordinates, p, X, B)
    return get_coordinates_diagonalizing!(M, Y, p, X, B)
end

function _get_coordinates(M::AbstractManifold, p, X, B::CachedBasis)
    return get_coordinates_cached(M, number_system(M), p, X, B, number_system(B))
end
function get_coordinates_cached(
    M::AbstractManifold,
    ::ComplexNumbers,
    p,
    X,
    B::CachedBasis,
    ::RealNumbers,
)
    return map(vb -> conj(inner(M, p, X, vb)), get_vectors(M, p, B))
end
function get_coordinates_cached(
    M::AbstractManifold,
    ::𝔽,
    p,
    X,
    C::CachedBasis,
    ::𝔽,
) where {𝔽}
    return map(vb -> real(inner(M, p, X, vb)), get_vectors(M, p, C))
end

function get_coordinates!(M::AbstractManifold, Y, p, X, B::AbstractBasis)
    return _get_coordinates!(M, Y, p, X, B)
end
function _get_coordinates!(M::AbstractManifold, Y, p, X, B::VeeOrthogonalBasis)
    return get_coordinates_vee!(M, Y, p, X, number_system(B))
end
function get_coordinates_vee!(M::AbstractManifold, Y, p, X, N)
    return get_coordinates!(M, Y, p, X, DefaultOrthogonalBasis(N))
end

function _get_coordinates!(M::AbstractManifold, Y, p, X, B::DefaultBasis)
    return get_coordinates_default!(M, Y, p, X, number_system(B))
end
function get_coordinates_default!(M::AbstractManifold, Y, p, X, N)
    return get_coordinates!(M, Y, p, X, DefaultOrthogonalBasis(N))
end

function _get_coordinates!(M::AbstractManifold, Y, p, X, B::DefaultOrthogonalBasis)
    return get_coordinates_orthogonal!(M, Y, p, X, number_system(B))
end
function get_coordinates_orthogonal!(M::AbstractManifold, Y, p, X, N)
    return get_coordinates!(M, Y, p, X, DefaultOrthonormalBasis(N))
end

function _get_coordinates!(M::AbstractManifold, Y, p, X, B::DefaultOrthonormalBasis)
    return get_coordinates_orthonormal!(M, Y, p, X, number_system(B))
end
function get_coordinates_orthonormal! end

function _get_coordinates!(M::AbstractManifold, Y, p, X, B::DiagonalizingOrthonormalBasis)
    return get_coordinates_diagonalizing!(M, Y, p, X, B)
end
function get_coordinates_diagonalizing! end

function _get_coordinates!(M::AbstractManifold, Y, p, X, B::CachedBasis)
    return get_coordinates_cached!(M, number_system(M), Y, p, X, B, number_system(B))
end
function get_coordinates_cached!(
    M::AbstractManifold,
    ::ComplexNumbers,
    Y,
    p,
    X,
    B::CachedBasis,
    ::RealNumbers,
)
    map!(vb -> conj(inner(M, p, X, vb)), Y, get_vectors(M, p, B))
    return Y
end
function get_coordinates_cached!(
    M::AbstractManifold,
    ::𝔽,
    Y,
    p,
    X,
    C::CachedBasis,
    ::𝔽,
) where {𝔽}
    map!(vb -> real(inner(M, p, X, vb)), Y, get_vectors(M, p, C))
    return Y
end

"""
    X = get_vector(M::AbstractManifold, p, c, B::AbstractBasis)

Convert a one-dimensional vector of coefficients in a basis `B` of
the tangent space at `p` on manifold `M` to a tangent vector `X` at `p`.

Depending on the basis, `p` may not directly represent a point on the manifold.
For example if a basis transported along a curve is used, `p` may be the coordinate
along the curve.

For the [`CachedBasis`](@ref) keep in mind that the reconstruction from [`get_coordinates`](@ref)
requires either a dual basis or the cached basis to be selfdual, for example orthonormal

See also: [`get_coordinates`](@ref), [`get_basis`](@ref)
"""
function get_vector(M::AbstractManifold, p, c, B::AbstractBasis)
    return _get_vector(M, p, c, B)
end

function _get_vector(M::AbstractManifold, p, c, B::VeeOrthogonalBasis)
    return get_vector_vee(M, p, c, number_system(B))
end
function get_vector_vee(M::AbstractManifold, p, c, N)
    return get_vector(M, p, c, DefaultOrthogonalBasis(N))
end

function _get_vector(M::AbstractManifold, p, c, B::DefaultBasis)
    return get_vector_default(M, p, c, number_system(B))
end
function get_vector_default(M::AbstractManifold, p, c, N)
    return get_vector(M, p, c, DefaultOrthogonalBasis(N))
end

function _get_vector(M::AbstractManifold, p, c, B::DefaultOrthogonalBasis)
    return get_vector_orthogonal(M, p, c, number_system(B))
end
function get_vector_orthogonal(M::AbstractManifold, p, c, N)
    return get_vector_orthonormal(M, p, c, N)
end

function _get_vector(M::AbstractManifold, p, c, B::DefaultOrthonormalBasis)
    return get_vector_orthonormal(M, p, c, number_system(B))
end
function get_vector_orthonormal(M::AbstractManifold, p, c, N)
    B = DefaultOrthonormalBasis(N)
    Y = allocate_result(M, get_vector, p, c)
    return get_vector!(M, Y, p, c, B)
end

function _get_vector(M::AbstractManifold, p, c, B::DiagonalizingOrthonormalBasis)
    return get_vector_diagonalizing(M, p, c, B)
end
function get_vector_diagonalizing(
    M::AbstractManifold,
    p,
    c,
    B::DiagonalizingOrthonormalBasis,
)
    Y = allocate_result(M, get_vector, p, c)
    return get_vector!(M, Y, p, c, B)
end

function _get_vector(M::AbstractManifold, p, c, B::CachedBasis)
    return get_vector_cached(M, p, c, B)
end
_get_vector_cache_broadcast(::Any) = Val(true)
function get_vector_cached(M::AbstractManifold, p, X, B::CachedBasis)
    # quite convoluted but:
    #  1) preserves the correct `eltype`
    #  2) guarantees a reasonable array type `Y`
    #     (for example scalar * `SizedValidation` is an `SArray`)
    bvectors = get_vectors(M, p, B)
    if _get_vector_cache_broadcast(bvectors[1]) === Val(false)
        Xt = X[1] * bvectors[1]
        for i in 2:length(X)
            copyto!(Xt, Xt + X[i] * bvectors[i])
        end
    else
        Xt = X[1] .* bvectors[1]
        for i in 2:length(X)
            Xt .+= X[i] .* bvectors[i]
        end
    end
    return Xt
end
function get_vector!(M::AbstractManifold, Y, p, c, B::AbstractBasis)
    return _get_vector!(M, Y, p, c, B)
end

function _get_vector!(M::AbstractManifold, Y, p, c, B::VeeOrthogonalBasis)
    return get_vector_vee!(M, Y, p, c, number_system(B))
end
get_vector_vee!(M, Y, p, c, N) = get_vector!(M, Y, p, c, DefaultOrthogonalBasis(N))

function _get_vector!(M::AbstractManifold, Y, p, c, B::DefaultBasis)
    return get_vector_default!(M, Y, p, c, number_system(B))
end
function get_vector_default!(M::AbstractManifold, Y, p, c, N)
    return get_vector!(M, Y, p, c, DefaultOrthogonalBasis(N))
end

function _get_vector!(M::AbstractManifold, Y, p, c, B::DefaultOrthogonalBasis)
    return get_vector_orthogonal!(M, Y, p, c, number_system(B))
end
function get_vector_orthogonal!(M::AbstractManifold, Y, p, c, N)
    return get_vector!(M, Y, p, c, DefaultOrthonormalBasis(N))
end

function _get_vector!(M::AbstractManifold, Y, p, c, B::DefaultOrthonormalBasis)
    return get_vector_orthonormal!(M, Y, p, c, number_system(B))
end
function get_vector_orthonormal! end

function _get_vector!(M::AbstractManifold, Y, p, c, B::DiagonalizingOrthonormalBasis)
    return get_vector_diagonalizing!(M, Y, p, c, B)
end
function get_vector_diagonalizing! end

function _get_vector!(M::AbstractManifold, Y, p, c, B::CachedBasis)
    return get_vector_cached!(M, Y, p, c, B)
end
function get_vector_cached!(M::AbstractManifold, Y, p, X, B::CachedBasis)
    # quite convoluted but:
    #  1) preserves the correct `eltype`
    #  2) guarantees a reasonable array type `Y`
    #     (for example scalar * `SizedValidation` is an `SArray`)
    bvectors = get_vectors(M, p, B)
    return if _get_vector_cache_broadcast(bvectors[1]) === Val(false)
        Xt = X[1] * bvectors[1]
        copyto!(Y, Xt)
        for i in 2:length(X)
            copyto!(Y, Y + X[i] * bvectors[i])
        end
        return Y
    else
        Xt = X[1] .* bvectors[1]
        copyto!(Y, Xt)
        for i in 2:length(X)
            Y .+= X[i] .* bvectors[i]
        end
        return Y
    end
end

"""
    get_vectors(M::AbstractManifold, p, B::AbstractBasis)

Get the basis vectors of basis `B` of the tangent space at point `p`.
"""
function get_vectors(M::AbstractManifold, p, B::AbstractBasis)
    return _get_vectors(M, p, B)
end
function _get_vectors(::AbstractManifold, ::Any, B::CachedBasis)
    return _get_vectors(B)
end
_get_vectors(B::CachedBasis{𝔽,<:AbstractBasis,<:AbstractArray}) where {𝔽} = B.data
function _get_vectors(B::CachedBasis{𝔽,<:AbstractBasis,<:DiagonalizingBasisData}) where {𝔽}
    return B.data.vectors
end

@doc raw"""
    gram_schmidt(M::AbstractManifold{𝔽}, p, B::AbstractBasis{𝔽}) where {𝔽}
    gram_schmidt(M::AbstractManifold, p, V::AbstractVector)

Compute an ONB in the tangent space at `p` on the [`AbstractManifold`](@ref} `M` from either an
[`AbstractBasis`](@ref) basis ´B´ or a set of (at most) [`manifold_dimension`](@ref)`(M)`
many vectors.
Note that this method requires the manifold and basis to work on the same
[`AbstractNumbers`](@ref) `𝔽`, i.e. with real coefficients.

The method always returns a basis, i.e. linearly dependent vectors are removed.

# Keyword arguments

* `warn_linearly_dependent` (`false`) – warn if the basis vectors are not linearly
  independent
* `skip_linearly_dependent` (`false`) – whether to just skip (`true`) a vector that
  is linearly dependent to the previous ones or to stop (`false`, default) at that point
* `return_incomplete_set` (`false`) – throw an error if the resulting set of vectors is not
  a basis but contains less vectors

further keyword arguments can be passed to set the accuracy of the independence test.
Especially `atol` is raised slightly by default to `atol = 5*1e-16`.

# Return value

When a set of vectors is orthonormalized a set of vectors is returned.
When an [`AbstractBasis`](@ref) is orthonormalized, a [`CachedBasis`](@ref) is returned.
"""
function gram_schmidt(
    M::AbstractManifold{𝔽},
    p,
    B::AbstractBasis{𝔽};
    warn_linearly_dependent = false,
    return_incomplete_set = false,
    skip_linearly_dependent = false,
    kwargs...,
) where {𝔽}
    V = gram_schmidt(
        M,
        p,
        get_vectors(M, p, B);
        warn_linearly_dependent = warn_linearly_dependent,
        skip_linearly_dependent = skip_linearly_dependent,
        return_incomplete_set = return_incomplete_set,
        kwargs...,
    )
    return CachedBasis(GramSchmidtOrthonormalBasis(𝔽), V)
end
function gram_schmidt(
    M::AbstractManifold,
    p,
    V::AbstractVector;
    atol = eps(number_eltype(first(V))),
    warn_linearly_dependent = false,
    return_incomplete_set = false,
    skip_linearly_dependent = false,
    kwargs...,
)
    N = length(V)
    Ξ = empty(V)
    dim = manifold_dimension(M)
    N < dim && @warn "Input only has $(N) vectors, but manifold dimension is $(dim)."
    @inbounds for n in 1:N
        Ξₙ = copy(M, p, V[n])
        for k in 1:length(Ξ)
            Ξₙ .-= real(inner(M, p, Ξ[k], Ξₙ)) .* Ξ[k]
        end
        nrmΞₙ = norm(M, p, Ξₙ)
        if isapprox(nrmΞₙ / dim, 0; atol = atol, kwargs...)
            warn_linearly_dependent &&
                @warn "Input vector $(n) lies in the span of the previous ones."
            !skip_linearly_dependent && throw(
                ErrorException("Input vector $(n) lies in the span of the previous ones."),
            )
        else
            push!(Ξ, Ξₙ ./ nrmΞₙ)
        end
        if length(Ξ) == dim
            (n < N) &&
                @warn "More vectors ($(N)) entered than the dimension of the manifold ($dim). All vectors after the $n th ignored."
            return Ξ
        end
    end
    if return_incomplete_set # if we rech this point - length(Ξ) < dim
        return Ξ
    else
        throw(
            ErrorException(
                "gram_schmidt found only $(length(Ξ)) orthonormal basis vectors, but manifold dimension is $(dim).",
            ),
        )
    end
end

@doc raw"""
    hat(M::AbstractManifold, p, Xⁱ)

Given a basis $e_i$ on the tangent space at a point `p` and tangent
component vector $X^i$, compute the equivalent vector representation
$X=X^i e_i$, where Einstein summation notation is used:

````math
∧ : X^i ↦ X^i e_i
````

For array manifolds, this converts a vector representation of the tangent
vector to an array representation. The [`vee`](@ref) map is the `hat` map's
inverse.
"""
@inline hat(M::AbstractManifold, p, X) = get_vector(M, p, X, VeeOrthogonalBasis())
@inline hat!(M::AbstractManifold, Y, p, X) = get_vector!(M, Y, p, X, VeeOrthogonalBasis())

"""
    number_of_coordinates(M::AbstractManifold{𝔽}, B::AbstractBasis)
    number_of_coordinates(M::AbstractManifold{𝔽}, ::𝔾)

Compute the number of coordinates in basis of field type `𝔾` on a manifold `M`.
This also corresponds to the number of vectors represented by `B`,
or stored within `B` in case of a [`CachedBasis`](@ref).
"""
function number_of_coordinates(M::AbstractManifold{𝔽}, ::AbstractBasis{𝔾}) where {𝔽,𝔾}
    return number_of_coordinates(M, 𝔾)
end
function number_of_coordinates(M::AbstractManifold{𝔽}, f::𝔾) where {𝔽,𝔾}
    # for odd manifolds this first case has to match.
    (real_dimension(𝔽) == real_dimension(f)) && return manifold_dimension(M)
    return div(manifold_dimension(M), real_dimension(𝔽)) * real_dimension(f)
end

"""
    number_system(::AbstractBasis)

The number system for the vectors of the given basis.
"""
number_system(::AbstractBasis{𝔽}) where {𝔽} = 𝔽

function _show_basis_vector(io::IO, X; pre = "", head = "")
    sX = sprint(show, "text/plain", X, context = io, sizehint = 0)
    sX = replace(sX, '\n' => "\n$(pre)")
    return print(io, head, pre, sX)
end
function _show_basis_vector_range(io::IO, Ξ, range; pre = "", sym = "E")
    for i in range
        _show_basis_vector(io, Ξ[i]; pre = pre, head = "\n$(sym)$(i) =\n")
    end
    return nothing
end
function _show_basis_vector_range_noheader(io::IO, Ξ; max_vectors = 4, pre = "", sym = "E")
    nv = length(Ξ)
    return if nv ≤ max_vectors
        _show_basis_vector_range(io, Ξ, 1:nv; pre = "  ", sym = " E")
    else
        halfn = div(max_vectors, 2)
        _show_basis_vector_range(io, Ξ, 1:halfn; pre = "  ", sym = " E")
        print(io, "\n ⋮")
        _show_basis_vector_range(io, Ξ, (nv - halfn + 1):nv; pre = "  ", sym = " E")
    end
end

function show(io::IO, ::DefaultBasis{𝔽}) where {𝔽}
    return print(io, "DefaultBasis($(𝔽))")
end
function show(io::IO, ::DefaultOrthogonalBasis{𝔽}) where {𝔽}
    return print(io, "DefaultOrthogonalBasis($(𝔽))")
end
function show(io::IO, ::DefaultOrthonormalBasis{𝔽}) where {𝔽}
    return print(io, "DefaultOrthonormalBasis($(𝔽))")
end
function show(io::IO, ::GramSchmidtOrthonormalBasis{𝔽}) where {𝔽}
    return print(io, "GramSchmidtOrthonormalBasis($(𝔽))")
end
function show(io::IO, ::ProjectedOrthonormalBasis{method,𝔽}) where {method,𝔽}
    return print(io, "ProjectedOrthonormalBasis($(repr(method)), $(𝔽))")
end
function show(io::IO, ::MIME"text/plain", onb::DiagonalizingOrthonormalBasis)
    println(
        io,
        "DiagonalizingOrthonormalBasis($(number_system(onb))) with eigenvalue 0 in direction:",
    )
    sk = sprint(show, "text/plain", onb.frame_direction, context = io, sizehint = 0)
    sk = replace(sk, '\n' => "\n ")
    return print(io, sk)
end
function show(
    io::IO,
    ::MIME"text/plain",
    B::CachedBasis{𝔽,T,D},
) where {𝔽,T<:AbstractBasis,D}
    vectors = _get_vectors(B)
    print(
        io,
        "Cached basis of type $T with $(length(vectors)) basis vector$(length(vectors) == 1 ? "" : "s"):",
    )
    return _show_basis_vector_range_noheader(
        io,
        vectors;
        max_vectors = 4,
        pre = "  ",
        sym = " E",
    )
end
function show(
    io::IO,
    ::MIME"text/plain",
    B::CachedBasis{𝔽,T,D},
) where {𝔽,T<:DiagonalizingOrthonormalBasis,D<:DiagonalizingBasisData}
    vectors = _get_vectors(B)
    nv = length(vectors)
    sk = sprint(show, "text/plain", T(B.data.frame_direction), context = io, sizehint = 0)
    sk = replace(sk, '\n' => "\n ")
    print(io, sk)
    println(io, "\nand $(nv) basis vector$(nv == 1 ? "" : "s").")
    print(io, "Basis vectors:")
    _show_basis_vector_range_noheader(io, vectors; max_vectors = 4, pre = "  ", sym = " E")
    println(io, "\nEigenvalues:")
    sk = sprint(show, "text/plain", B.data.eigenvalues, context = io, sizehint = 0)
    sk = replace(sk, '\n' => "\n ")
    return print(io, ' ', sk)
end

@doc raw"""
    vee(M::AbstractManifold, p, X)

Given a basis $e_i$ on the tangent space at a point `p` and tangent
vector `X`, compute the vector components $X^i$, such that $X = X^i e_i$, where
Einstein summation notation is used:

````math
\vee : X^i e_i ↦ X^i
````

For array manifolds, this converts an array representation of the tangent
vector to a vector representation. The [`hat`](@ref) map is the `vee` map's
inverse.
"""
vee(M::AbstractManifold, p, X) = get_coordinates(M, p, X, VeeOrthogonalBasis())
vee!(M::AbstractManifold, Y, p, X) = get_coordinates!(M, Y, p, X, VeeOrthogonalBasis())


text:
using Smg2s
using Test
using SparseArrays
using Plots

using LinearAlgebra
using Random

n = 10000
diag_l = -9900
diag_u = -800
nbOne = 70
nilpdiag = 3500
sp=0.0002
figname = "fig/pattern12"

nilpvec = zeros(Int64, n-1)

cnt = 0

while cnt < n-1
    l₁ = rand(1:nbOne)
    l₀ = rand(1:nbOne)

    for i = cnt+1:l₁+cnt
        if i <= n-1
            nilpvec[i] = 1
        end
    end

    global cnt += l₀ + l₁

end

#nilpMatrix = sparse(diagm(nilpdiag=>nilpvec[1:n-nilpdiag]))

#nilp = Nilp(nilpMatrix, n)
nilp = Nilp(nilpvec, nilpdiag, n)

function f(idx::Integer, n = n)
    return 10 * cos(((idx-1) / 2) * 2 * π / n) + 5
end

spec = zeros(ComplexF64, n)
Spectrum!(spec, f, n)

Am = spzeros(Float64, n, n)
initMat!(Am, diag_l, diag_u, n; scale=0.1, sparsity=sp)
@time genMat = nonsym(n, diag_l, diag_u, spec, Am, nilp)
sparsity = nnz(genMat) / (n * n)
@info "sparsity=" sparsity
spy(genMat, legend = :none, title="size:$n, diag_l:$diag_l, diag_u:$diag_u, nbOne:$nbOne, nilpdiag: $nilpdiag, sp: $sp",titlefont = font(8))