text:
using MRIgeneralizedBloch
using Test

@testset "Greens_Interpolation" begin
    include("Greens_Interpolation.jl")
end

@testset "apply_hamiltonian_gbloch!" begin
    for i = 1:10 # test different random initializations
        include("gBloch_Hamiltonian.jl")
    end
end

@testset "apply_hamiltonian_freeprecession!" begin
    for i = 1:10 # test different random initializations
        include("FreePrecession_Hamiltonian.jl")
    end
end

@testset "apply_hamiltonian_graham_superlorentzian!" begin
    for i = 1:10 # test different random initializations
        include("Graham_Hamiltonian_superLorentzian.jl")
    end
end

@testset "apply_hamiltonian_sled!" begin
    for i = 1:10 # test different random initializations
        include("Sled_Hamiltonian_superLorentzian.jl")
    end
end

@testset "Linear_Approximation" begin
    include("Linear_Approximation.jl")
end

@testset "apply_hamiltonian_gbloch! gradients" begin
    include("gBloch_Hamiltonian_Gradients.jl")
end

@testset "apply_hamiltonian_freeprecession! gradients" begin
    include("FreePrecession_Hamiltonian_Gradients.jl")
end

@testset "Linear_Approx gradients" begin
    include("Linear_Approx_Gradients.jl")
end

@testset "apply_hamiltonian_graham_superlorentzian! gradients" begin
    include("Graham_Hamiltonian_Gradients.jl")
end

@testset "inversionpulse" begin
    include("inversionpulse.jl")
end

@testset "cf. Solvers" begin
    include("Solvers_cf.jl")
end

@testset "Solvers_Gradients" begin
    include("Solvers_Gradients.jl")
end

@testset "OCT Gradients" begin
    for i = 1:10 # test different random initializations
        include("OCT.jl")
    end
end


text:
using lmp_str

file_out = "test.data"

data_cell = genr_cell([5, 5, 5])
str = structure_si3n4()
data = genr(data_cell, str)
write_data(data, file_out)

ext: jl


text:
"""
    backsub(A,B)

adalah fungsi untuk menyelesaikan SPL `Ax=B` dengan `A` adalah matriks segitiga
atas dan `B` adalah vektor ruas kanan.

# Examples
```jldoctest
julia> A = [ 4 -1 2 3
             0 -2 7 -4
             0 0 6 5
             0 0 0 3];

julia> B = [20;-7; 4; 6];

julia> X = backsub(A,B)
4×1 Array{Float64,2}:
  3.0
 -4.0
 -1.0
  2.0
```
return solusi SPL `X`
"""
function backsub(A,B)
    # Hitung ukuran matriks dan inisiasi solusi X
    n = length(B);
    X = zeros(n,1);
    # Hitung nilai solusi X ke-n
    X[n] = B[n]/A[n,n];
    # Hitung nilai solusi X ke- n-1 sampai X ke-1
    for i = n-1:-1:1
        X[i] = (B[i] - A[i,i+1:n]'X[i+1:n])/A[i,i];
    end
    return X
end


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

julia_version = v"1.6.3"

name = "OpenSpiel"
version = v"1.1.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/deepmind/open_spiel/archive/v1.1.1.tar.gz", "73d46e2d9ce7a86a420aad141b789e98237463b5c98608401faf86b15ec5eaf9"),
    ArchiveSource("https://github.com/findmyway/dds/archive/v0.1.1.tar.gz", "fd13ee77feb5b5c3dfcc3333a0523266beb2a3d27715703cf508313af25306e5"),
    ArchiveSource("https://github.com/abseil/abseil-cpp/archive/20211102.0.tar.gz", "dcf71b9cba8dc0ca9940c4b316a0c796be8fab42b070bb6b7cab62b48f0e66c4"),
    ArchiveSource("https://github.com/findmyway/hanabi-learning-environment/archive/v0.1.0.tar.gz", "6126936fd13a95f8cadeacaa69dfb38a960eaf3bd588aacc8893a6e07e4791a3"),
    ArchiveSource("https://github.com/findmyway/project_acpc_server/archive/v0.1.0.tar.gz", "e29f969dd62ba354b7019cae3f7f1dbfbd9a744687ea4a8f7494c2bb1ee87382"),
]

# Bash recipe for building across all platforms
script = raw"""

mv open_spiel-* open_spiel
mv abseil-cpp-* open_spiel/open_spiel/abseil-cpp
mv dds-* open_spiel/open_spiel/games/bridge/double_dummy_solver
mv hanabi-learning-environment-* open_spiel/open_spiel/games/hanabi/hanabi-learning-environment
mv project_acpc_server-* open_spiel/open_spiel/games/universal_poker/acpc

mkdir open_spiel/build
cd open_spiel/build
export OPEN_SPIEL_BUILD_WITH_JULIA=ON OPEN_SPIEL_BUILD_WITH_PYTHON=OFF OPEN_SPIEL_BUILD_WITH_HANABI=ON OPEN_SPIEL_BUILD_WITH_ACPC=OFF
cmake \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    ../open_spiel/
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/open_spiel/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libspieljl", :libspieljl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("libjulia_jll"),
    Dependency("libcxxwrap_julia_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"9",
    julia_compat="$(julia_version.major).$(julia_version.minor)"
)


text:
# This file was generated, do not modify it. # hide
ŷ = MLJ.predict(LinearModel, X)
yhatResponse = [ŷ[i,1].μ for i in 1:nrow(y)]
residuals = y .- yhatResponse
r = report(LinearModel)

k = collect(keys(fp.fitted_params_given_machine))[3]
println("\n Coefficients:  ", fp.fitted_params_given_machine[k].coef)
println("\n y \n ", y[1:5,1])
println("\n ŷ \n ", ŷ[1:5])
println("\n yhatResponse \n ", yhatResponse[1:5])
println("\n Residuals \n ", y[1:5,1] .- yhatResponse[1:5])
println("\n Standard Error per Coefficient \n", r.report_given_machine[k].stderror)


text:
# This file contains functions for the Radial DEA model
"""
    AbstractRadialDEAModel
An abstract type representing a radial DEA model.
"""
abstract type AbstractRadialDEAModel <: AbstractTechnicalDEAModel end

"""
    RadialDEAModel
An data structure representing a radial DEA model.
"""
struct RadialDEAModel <: AbstractRadialDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    disposX::Symbol
    disposY::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dea(X, Y)
Compute the radial model using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `disposX=:Strong`: chooses strong disposability of inputs. For weak disposability choose `:Weak`.
- `disposY=:Strong`: chooses strong disposability of outputs. For weak disposability choose `:Weak`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> dea(X, Y)
Radial DEA Model
DMUs = 11; Inputs = 2; Outputs = 1
Orientation = Input; Returns to Scale = CRS
──────────────────────────────────────────────────
    efficiency       slackX1      slackX2  slackY1
──────────────────────────────────────────────────
1     1.0        0.0          0.0              0.0
2     0.62229   -4.41868e-15  0.0              0.0
3     0.819856   0.0          8.17926e-15      0.0
4     1.0       -8.03397e-16  0.0              0.0
5     0.310371   1.80764e-15  0.0              0.0
6     0.555556   4.44444      0.0              0.0
7     1.0        0.0          0.0              0.0
8     0.757669   1.60679e-15  0.0              0.0
9     0.820106   1.64021      0.0              0.0
10    0.490566   9.68683e-15  0.0              0.0
11    1.0        0.0          4.0              0.0
──────────────────────────────────────────────────
```
"""
function dea(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    orient::Symbol = :Input, rts::Symbol = :CRS, slack::Bool = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    disposX::Symbol = :Strong, disposY::Symbol = :Strong,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing, progress::Bool = false)::RadialDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    if disposX != :Strong && disposX != :Weak
        throw(ArgumentError("`disposX` must be :Strong or :Weak"));
    end

    if disposY != :Strong && disposY != :Weak
        throw(ArgumentError("`disposY` must be :Strong or :Weak"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    lambdaeff = spzeros(n, nref)

    p = progressbarmeter(n, desc = "Computing Radial DEA Model...", progress = progress)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, eff)
        @variable(deamodel, lambda[1:nref] >= 0)

        if orient == :Input
            # Input orientation
            @objective(deamodel, Min, eff)

            # Inequality or equality restrictions based on disposability
            if disposX == :Strong
                @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= eff * x0[j])
            elseif disposX == :Weak
                @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) == eff * x0[j])
            end

            if disposY == :Strong
                @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j])
            elseif disposY == :Weak
                @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) == y0[j])
            end

        elseif orient == :Output
            # Output orientation
            @objective(deamodel, Max, eff)

            # Inequality or equality restrictions based on disposability
            if disposX == :Strong
                @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= x0[j])
            elseif disposX == :Weak
                @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) == x0[j])
            end

            if disposY == :Strong
                @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= eff * y0[j])
            elseif disposY == :Weak
                @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) == eff * y0[j])
            end

        else
            throw(ArgumentError("`orient` must be :Input or :Output"));
        end

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

        next!(p)
    end

    # Get first-stage X and Y targets
    if orient == :Input
        Xtarget = X .* effi
        Ytarget = Y
    elseif orient == :Output
        Xtarget = X
        Ytarget = Y .* effi
    end

    # Compute slacks
    if slack == true

        # Use additive model with radial efficient X and Y to get slacks
        if disposX == :Strong
            rhoX = ones(size(X))
        elseif disposX == :Weak
            rhoX = zeros(size(X))
        end

        if disposY == :Strong
            rhoY = ones(size(Y))
        elseif disposY == :Weak
            rhoY = zeros(size(Y))
        end

        slacksmodel = deaadd(Xtarget, Ytarget, rhoX = rhoX, rhoY = rhoY, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
        slackX = slacks(slacksmodel, :X)
        slackY = slacks(slacksmodel, :Y)

        # Get second-stage X and Y targets
        Xtarget = Xtarget - slackX
        Ytarget = Ytarget + slackY
    else
        if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
        if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    return RadialDEAModel(n, m, s, orient, rts, disposX, disposY, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::RadialDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    disposX = x.disposX
    disposY = x.disposY
    dmunames = names(x)

    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Radial DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if disposX == :Weak print(io, "Weak disposability of inputs \n") end
        if disposY == :Weak print(io, "Weak disposability of outputs \n") end

        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], dmunames))
        end
    end

end


text:
"""
    AdamsBashforthTimeStepper(float_type, arch, grid, tracers, χ=0.125;
                              Gⁿ = TendencyFields(arch, grid, tracers),
                              G⁻ = TendencyFields(arch, grid, tracers))

Return an AdamsBashforthTimeStepper object with tendency fields on `arch` and
`grid` with AB2 parameter `χ`. The tendency fields can be specified via optional
kwargs.
"""
struct AdamsBashforthTimeStepper{T, TG, P} <: AbstractTimeStepper
     χ :: T
    Gⁿ :: TG
    G⁻ :: TG
    predictor_velocities :: P
end

function AdamsBashforthTimeStepper(float_type, arch, grid, velocities, tracers, χ=0.1;
                                   Gⁿ = TendencyFields(arch, grid, tracers),
                                   G⁻ = TendencyFields(arch, grid, tracers))

    u★ = Field{Face, Cell, Cell}(data(velocities.u), grid, velocities.u.boundary_conditions)
    v★ = Field{Cell, Face, Cell}(data(velocities.v), grid, velocities.v.boundary_conditions)
    w★ = Field{Cell, Cell, Face}(data(velocities.w), grid, velocities.w.boundary_conditions)

    U★ = (u=u★, v=v★, w=w★)

    return AdamsBashforthTimeStepper{float_type, typeof(Gⁿ), typeof(U★)}(χ, Gⁿ, G⁻, U★)
end

#####
##### Time steppping
#####

"""
    time_step!(model::IncompressibleModel{<:AdamsBashforthTimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and
pressure-correction substep. Setting `euler=true` will take a forward Euler time step.
"""
function time_step!(model::IncompressibleModel{<:AdamsBashforthTimeStepper}, Δt; euler=false)
    χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ)

    # Convert NamedTuples of Fields to NamedTuples of OffsetArrays
    velocities, tracers, pressures, diffusivities, Gⁿ, G⁻, predictor_velocities =
        datatuples(model.velocities, model.tracers, model.pressures, model.diffusivities,
                   model.timestepper.Gⁿ, model.timestepper.G⁻, model.timestepper.predictor_velocities)

    calculate_explicit_substep!(Gⁿ, velocities, tracers, pressures, diffusivities, model)

    ab2_time_step_tracers!(tracers, model.architecture, model.grid, Δt, χ, Gⁿ, G⁻)

    # Fractional step. Note that predictor velocities share memory space with velocities.
    ab2_update_predictor_velocities!(predictor_velocities, model.architecture, model.grid, Δt, χ, Gⁿ, G⁻)

    calculate_pressure_correction!(pressures.pNHS, Δt, predictor_velocities, model)

    fractional_step_velocities!(velocities, tracers, model.architecture,
                                model.grid, Δt, pressures.pNHS)

    ab2_store_source_terms!(G⁻, model.architecture, model.grid, χ, Gⁿ)

    # Compute w from recontinuity
    fill_halo_regions!(model.velocities, model.architecture,
                       boundary_condition_function_arguments(model)...)

    compute_w_from_continuity!(model)

    tick!(model.clock, Δt)

    return nothing
end

#####
##### Source term storage
#####

""" Store previous source terms before updating them. """
function ab2_store_source_terms!(G⁻, arch, grid, χ, Gⁿ)

    # Velocity fields
    @launch(device(arch), config=launch_config(grid, :xyz),
            ab2_store_velocity_source_terms!(G⁻, grid, χ, Gⁿ))

    # Tracer fields
    for i in 4:length(G⁻)
        @inbounds Gc⁻ = G⁻[i]
        @inbounds Gcⁿ = Gⁿ[i]
        @launch(device(arch), config=launch_config(grid, :xyz),
                ab2_store_tracer_source_term!(Gc⁻, grid, χ, Gcⁿ))
    end

    return nothing
end

""" Store source terms for `u`, `v`, and `w`. """
function ab2_store_velocity_source_terms!(G⁻, grid::AbstractGrid{FT}, χ, Gⁿ) where FT
    @loop_xyz i j k grid begin
        @inbounds G⁻.u[i, j, k] = Gⁿ.u[i, j, k]
        @inbounds G⁻.v[i, j, k] = Gⁿ.v[i, j, k]
        @inbounds G⁻.w[i, j, k] = Gⁿ.w[i, j, k]
    end
    return nothing
end

""" Store previous source terms for a tracer before updating them. """
function ab2_store_tracer_source_term!(Gc⁻, grid::AbstractGrid{FT}, χ, Gcⁿ) where FT
    @loop_xyz i j k grid begin
        @inbounds Gc⁻[i, j, k] = Gcⁿ[i, j, k]
    end
    return nothing
end

function ab2_time_step_tracers!(C, arch, grid, Δt, χ, Gⁿ, G⁻)
    for i in 1:length(C)
        @inbounds c = C[i]
        @inbounds Gcⁿ = Gⁿ[i+3]
        @inbounds Gc⁻ = G⁻[i+3]
        @launch device(arch) config=launch_config(grid, :xyz) ab2_time_step_tracer!(c, grid, Δt, χ, Gcⁿ, Gc⁻)
    end

    return nothing
end

#####
##### Tracer time stepping and predictor velocity updating
#####

"""
Time step tracers via

    `c^{n+1} = c^n + Δt ( (3/2 + χ) * Gc^{n} - (1/2 + χ) G^{n-1} )`

"""
function ab2_time_step_tracer!(c, grid::AbstractGrid{FT}, Δt, χ, Gcⁿ, Gc⁻) where FT
    @loop_xyz i j k grid begin
        @inbounds c[i, j, k] += Δt * ((FT(1.5) + χ) * Gcⁿ[i, j, k] - (FT(0.5) + χ) * Gc⁻[i, j, k])
    end
    return nothing
end

function ab2_update_predictor_velocities!(U★, arch, grid, Δt, χ, Gⁿ, G⁻)
    @launch(device(arch), config=launch_config(grid, :xyz),
            _ab2_update_predictor_velocities!(U★, grid, Δt, χ, Gⁿ, G⁻))
    return nothing
end

""" Update predictor velocity field. """
function _ab2_update_predictor_velocities!(U★, grid::AbstractGrid{FT}, Δt, χ, Gⁿ, G⁻) where FT
    @loop_xyz i j k grid begin
        @inbounds begin
            U★.u[i, j, k] += Δt * (   (FT(1.5) + χ) * Gⁿ.u[i, j, k]
                                    - (FT(0.5) + χ) * G⁻.u[i, j, k] )

            U★.v[i, j, k] += Δt * (   (FT(1.5) + χ) * Gⁿ.v[i, j, k]
                                    - (FT(0.5) + χ) * G⁻.v[i, j, k] )

            U★.w[i, j, k] += Δt * (   (FT(1.5) + χ) * Gⁿ.w[i, j, k]
                                    - (FT(0.5) + χ) * G⁻.w[i, j, k] )
        end
    end
    return nothing
end


text:
using Test
using PotentialCalculation.IdenticalTools
using PotentialCalculation.Atoms
using PotentialCalculation.Clusters

@testset "Clusters" begin
    formic_acid=Cluster(
     [-6.7041359778      1.3501192944      0.0102209137
    -5.3688853815      1.2229556023      0.0440598937
    -7.2470157373      2.4374213225      0.0651311769
     -5.0398812618      2.1435406993      0.1155201154
     -7.2001330967      0.3718768293     -0.0703451879]',
      AtomOnlySymbol.(["C", "O", "O", "H", "H"]) )

     Ar = Cluster(rand(3), AtomOnlySymbol("Ar"))

     nf = ClusterNoSymbols(rand(3,5))
     nar = ClusterNoSymbols(rand(3))

     @test length(nf) == length(formic_acid)

     t1 = formic_acid[5]
     t2 = formic_acid[end]

     @test t1.atoms == t2.atoms && t1.xyz == t2.xyz

     tn = nf[3]

     @test all(tn.xyz .== nf.xyz[:,3])

     @test length(formic_acid[1:2]) == 2
     @test length(nf[2:4]) == 3

     @test length(nf+nar) == 6
     @test length(formic_acid + Ar) == 6

     distances(formic_acid, Ar)
     distances(formic_acid, nf)
     distances(formic_acid, 1:2, 3:4)

     distances(formic_acid, 1, nf, 2)
     distances(formic_acid, 1:2, nf, 2:4)

     show(devnull, nf)
     print(devnull, formic_acid)

     move!(nf, [10,0,0])

     for x in (formic_acid, nf)
         cc = center_coordinates(x)
         center_cluster!(x)
         @test all(abs.(center_coordinates(x)) .< 1E-12 .* ones(size(cc)))
     end

     l = distances(formic_acid, 2,4)
     a = cluster_angle(formic_acid, 1,2,4)
     d = dihedral_angle(formic_acid, 3,1,2,4)

     cluster_angle(formic_acid, 1,2, nf, 3)

     @test l != distances(formic_acid, 1,2)
     @test a != cluster_angle(formic_acid, 1,2,3)
     @test d != dihedral_angle(formic_acid, 1,2,3,5)

     for f in [rotate_x!, rotate_y!, rotate_z!]
         f(formic_acid, rand())
         tl = distances(formic_acid, 2,4)
         ta = cluster_angle(formic_acid, 1,2,4)
         td = dihedral_angle(formic_acid, 3,1,2,4)
         @test l ≈ tl && a ≈ ta && d ≈ td
     end

     distances(nf, 2,4)
     cluster_angle(nf, 1,2,4)
     dihedral_angle(nf, 3,1,2,4)

     @test_throws DimensionMismatch ClusterNoSymbols(rand(4))
     @test_throws DimensionMismatch ClusterNoSymbols(rand(2))
     @test_throws DimensionMismatch ClusterNoSymbols(rand(2,3))
     @test_throws DimensionMismatch ClusterNoSymbols(rand(4,2))

     @test_throws DimensionMismatch Cluster(rand(2), [AtomOnlySymbol("H")])
     @test_throws DimensionMismatch Cluster(rand(2), AtomOnlySymbol("H"))
     @test_throws DimensionMismatch Cluster(rand(3),AtomOnlySymbol.(["H", "O"]))
     @test_throws DimensionMismatch Cluster(rand(3,3),AtomOnlySymbol.(["H", "O"]))
     @test_throws DimensionMismatch Cluster(rand(4,2),AtomOnlySymbol.(["H", "O"]))
end


text:
""" Working with the LowerBound (prefix is LB for Lower Bound) """

# Public functions :
# - LB_dichoSearch : computes X_SE1
# - LB_update! : add potential new points to the given LB
# - LB_getNadirs : computes the local nadir points of a given LB

#-------------------------------------------------------------------#
#------------------------ PUBLIC FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Returns the set LowerBound object (`X_SE1`) associated with the given problem `prob` and `assignment`, computed by the dichotomic method.

Two solvers are available : `Combo` (EXACT_COMBO) and `GLPK` (EXACT_JUMP).

To speed up the search, a parameter `parentToChild` containing a list of already known solutions can be provided.
"""
function LB_dichoSearch(prob::BiOKP, assignment::Assignment; parentToChild = ParentToChild())

    CONFIG.debug && DEBUG_parentToChild(prob, parentToChild)

    LB = LowerBound() # init empty lower bound

    LB.sols = copy(parentToChild.knownSols)

    # compute left and right solutions
    leftSol, rightSol = getLexicoSolutions(prob, assignment, LB.sols, getSolver1OKP(COMPONENTS.firstDicho), parentToChild = parentToChild)

    # In the event that the two solutions found are identical, this solution is optimal and thus added to the lower bound
    if !differentFloatSols(leftSol,rightSol)
        if LB.sols == nil(Sol) # if the LB contains already one solution, we don't need to add this one because it is the same
            LB.sols = cons(leftSol, LB.sols)
        end
    else # the two lexicographic solutions are different

        # get the list of pairs of solutions to study
        toStudy, LB.sols = computesPairsToStudy(prob, LB.sols, parentToChild, leftSol, rightSol)

		# while some pairs of solutions are found by the dichotomic search, we keep going
        while toStudy != []

			currPair = pop!(toStudy)

			leftSol = currPair.solL
			rightSol = currPair.solR

	        λ = [leftSol.y[2] - rightSol.y[2], rightSol.y[1] - leftSol.y[1]]

            solver1OKP = getSolver1OKP(COMPONENTS.firstDicho)
            midSol = solver1OKP(prob, λ, assignment)

			# if the solution is beyond the pair studied, it's added to the LB
			if round(sum(λ .* midSol.y), digits = 7) - round(sum(λ .* rightSol.y), digits = 7) > 0.00001
                # add the new solution found
				strictInsertInFront!(prob, LB.sols, midSol)
                # add the two newly found pairs
				push!(toStudy, PairOfSolution(leftSol, midSol))
				push!(toStudy, PairOfSolution(midSol, rightSol))
            end
        end
    end

    CONFIG.debug && DEBUG_LB(prob, LB)

	return LB

end

"""
Update the given `LB` and returns the new vector of nadir points to study `newNadirsToStudy` according to the given `UBSols` and `dominatedNadirs`.

`dominatedNadirs` are the nadir points to study.
"""
function LB_update!(prob::BiOKP, LB::LowerBound, dominatedNadirs::Vector{PairOfSolution}, UB::UpperBound)

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_UB(prob, UB)
    CONFIG.debug && DEBUG_nadirs(prob, dominatedNadirs)

    for sol in UB.sols
        # check if the solution found in the upper bound is binary (thus feasible)
        if !(COMPONENTS.methodUB == RELAX_LIN_CLASSIC) || sol.isBinary

            LB_removeAllDominatedSols!(prob, LB, sol)

            # insert sol at the right place (if the solution is not already dominated by the current LB)
            inserted = LB_insert!(prob, LB, sol)

        end
    end

    newNadirsToStudy = LB_getNadirsWithNadirsToStudy(prob, LB, dominatedNadirs)

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_nadirs(prob, newNadirsToStudy)

    return newNadirsToStudy
end

"""
Returns the nadir points `nadirs` associated with the given `LB`.
"""
function LB_getNadirs(prob::BiOKP, LB::LowerBound)

    CONFIG.debug && DEBUG_LB(prob, LB)

    sols = LB.sols
    nadirs = Vector{PairOfSolution}(undef, length(sols) - 1)

    iterNadirs = 1
    # go through sols and simply create each nadir point
    while sols.tail != nil(Sol)
    	nadirs[iterNadirs] = PairOfSolution(sols.head, sols.tail.head)
    	sols = sols.tail
    	iterNadirs += 1
    end

    CONFIG.debug && DEBUG_nadirs(prob, nadirs)

    return nadirs
end

#-------------------------------------------------------------------#
#----------------------- PRIVATE FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Inserts in place (if possible) the given solution `sol` into the given `LB`.

The solution `sol` is assumed to `NOT DOMINATE` solutions inside `LB.sols`.
"""
function LB_insert!(prob::BiOKP, LB::LowerBound, sol::Sol)

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_feasibleBinarySolution(prob, sol)
    CONFIG.debug && DEBUG_isntDominatingAny(sol, LB.sols)

    ll = LB.sols
    while ll != nil(Sol)
        # sol is dominated by the head of LB.sols or is the same, we shouldn't insert it
        if dominate(ll.head,sol) || (prob.isInteger && sol.y == ll.head.y) || (!prob.isInteger && !differentFloatSols(sol,ll.head))
            CONFIG.debug && DEBUG_LB(prob, LB)
            return false

        # we have to insert before the head
        elseif sol.y[1] < ll.head.y[1]
            oldSol = ll.head
            ll.head = sol
            CONFIG.debug && DEBUG_LB(prob, LB)
            return LB_insert!(prob, ll, oldSol)

        else

            if ll.tail != nil(Sol) # we can keep going

                # sol is dominated by the following solution or is the same, we shouldn't insert it
                if dominate(ll.tail.head,sol) || (prob.isInteger && sol.y == ll.tail.head.y) || (!prob.isInteger && !differentFloatSols(sol,ll.tail.head))
                    CONFIG.debug && DEBUG_LB(prob, LB)
                    return false

                # we have to insert sol between ll and ll.tail
                elseif sol.y[1] < ll.tail.head.y[1] # we have to insert sol between ll and ll.tail
                    oldTail = ll.tail
                    ll.tail = cons(sol,oldTail)
                    CONFIG.debug && DEBUG_LB(prob, LB)
                    return true
                end

            else # we are at the end of the list, we have to insert sol after ll
                ll.tail = cons(sol,nil(Sol))
                CONFIG.debug && DEBUG_LB(prob, LB)
                return true
            end

        end

        ll = ll.tail
    end

    CONFIG.debug && DEBUG_LB(prob, LB)

    return false
end

"""
Remove in place from the given `LB`, solutions that are dominated by the given solution `sol`, and returns wether or not a solution from `LB` has been removed.
"""
function LB_removeAllDominatedSols!(prob::BiOKP, LB::LowerBound, sol::Sol)

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_feasibleBinarySolution(prob, sol)

    function LB_removeOneDominatedSol!(prob::BiOKP, front::LinkedList{Sol}, sol::Sol)

        ll = front

        removed = false
        if dominate(sol,ll.head)
            if ll.tail != nil(Sol)
                if ll.tail.tail != nil(Sol)
                    ll.head = ll.tail.head
                    ll.tail = ll.tail.tail
                else
                    ll.head = ll.tail.head
                    ll.tail = nil(Sol)
                end
                removed = true
            end
        end

        if ll != nil(Sol) && ll.tail != nil(Sol) && ll.tail.tail == nil(Sol)
            if dominate(sol,ll.tail.head)
                removed = true
                ll.tail = nil(Sol)
            end
        end

        CONFIG.debug && DEBUG_front(prob, front)

        return removed

    end

    # delete in LB the solutions dominated by sol
    front = LB.sols
    while front != nil(Sol)
        removed = LB_removeOneDominatedSol!(prob, front, sol)
        if !removed
            front = front.tail
        end
    end

    CONFIG.debug && DEBUG_LB(prob, LB)

end

"""
Returns the nadir points `nadirs` associated with the given `LB`, but this time we care only about nadirs that are within the `nadirsToStudy` list.
"""
function LB_getNadirsWithNadirsToStudy(prob::BiOKP, LB::LowerBound, nadirsToStudy::Vector{PairOfSolution})

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_nadirs(prob, nadirsToStudy)

    # if the nadirsToStudy list is empty, just call the default function LB_getNadirs
    if length(nadirsToStudy) == 0
        return LB_getNadirs(prob, LB)
    end

    nadirs = Vector{PairOfSolution}()

    sols = LB.sols

    iterNadirs = 1
    # go through the sols and push certain nadirs
    while sols.tail != nil(Sol) && iterNadirs <= length(nadirsToStudy)

        # store the two current solutions from the sols
        solLBLeft = sols.head
        solLBRight = sols.tail.head
        # store the coordinates of the nadir associated with the two solutions
        z1,z2 = solLBLeft.y[1],solLBRight.y[2]

        # check if the nadir is to be stored
        if z1 >= nadirsToStudy[iterNadirs].solL.y[1] && z2 >= nadirsToStudy[iterNadirs].solR.y[2]
            push!(nadirs,PairOfSolution(solLBLeft,solLBRight))
            sols = sols.tail

        # check if the nadir won't even be stored
        elseif z1 < nadirsToStudy[iterNadirs].solL.y[1] && z2 > nadirsToStudy[iterNadirs].solR.y[2]
            sols = sols.tail

        # check if the nadir will eventually be stored
        elseif z1 > nadirsToStudy[iterNadirs].solL.y[1] && z2 < nadirsToStudy[iterNadirs].solR.y[2]
            iterNadirs += 1
        end
    end

    CONFIG.debug && DEBUG_nadirs(prob, nadirs)

    return nadirs
end


text:
import LinearAlgebra

n = 2
nProc = numProc(comm)
Data = Float32

######Build matrix######
map = BlockMap(nProc*n, comm)

numMyElts = numMyElements(map)
numGlobalElts = numGlobalElements(map)
myGlobalElts = myGlobalElements(map)

numNz = Array{GID, 1}(undef, numMyElts)
for i = 1:numMyElts
    if myGlobalElts[i] == 1 || myGlobalElts[i] == numGlobalElts
        numNz[i] = 2
    else
        numNz[i] = 3
    end
end

const A = CSRMatrix{Data}(map, numNz, STATIC_PROFILE)

const values = Data[-1, -1]
two = Data[2]

for i = 1:numMyElts
    if myGlobalElts[i] == 1
        indices = GID[2]
    elseif myGlobalElts[i] == numGlobalElts
        indices = GID[numGlobalElts-2]
    else
        indices = GID[myGlobalElts[i]-1, myGlobalElts[i]+1]
    end

    insertGlobalValues(A, myGlobalElts[i], indices, values)
    insertGlobalValues(A, myGlobalElts[i], GID[myGlobalElts[i]], two)
end

fillComplete(A, map, map)

Y = DenseMultiVector(map, LinearAlgebra.diagm(0 => Data(1):n))
X = DenseMultiVector(map, fill(Data(2), n, n))

@test Y === apply!(Y, A, X, NO_TRANS, Float32(3), Float32(.5))

@test fill(2, n, n) == X.data #ensure X isn't mutated

exp = LinearAlgebra.diagm(0 => Data(1):n)*.5
for i in 1:n
    if i == 1 && pid == 1
        exp[1, :] .+= 6
    elseif i == n && pid == nProc
        exp[i, :] .+= 6
    #else
        #exp[i, :] .+= -6 +12-6
    end
end

@test exp == Y.data

X = DenseMultiVector(map, fill(Data(2), n, n))
Y = A*X

@test Y isa DenseMultiVector
@test numVectors(Y) == n
@test fill(2, n, n) == X.data #ensure X isn't mutated
exp = fill(Data(0), n, n)
for i in 1:n
    if i == 1 && pid == 1
        exp[1, :] .= 2
    elseif i == n && pid == nProc
        exp[i, :] .= 2
    #else
        #exp[i, :] .= -2 +4 -2
    end
end

@test exp == Y.data

Y = DenseMultiVector(map, fill(Data(2), n, n))
X = DenseMultiVector(map, fill(Data(2), n, n))

@test Y === mul!(Y, A, X)

@test fill(2, n, n) == X.data #ensure X isn't mutated

exp = fill(Data(0), n, n)
for i in 1:n
    if i == 1 && pid == 1
        exp[1, :] .+= 2
    elseif i == n && pid == nProc
        exp[i, :] .+= 2
    #else
        #exp[i, :] .+= -2 +4 -2
    end
end

@test exp == Y.data


text:
module TestData

# using Revise
using Test
using MLJBase
using DataFrames
using TypedTables
using StatsBase
# using JuliaDB
using SparseArrays
using CategoricalArrays

using Random
import Random.seed!
seed!(1234)

import MLJBase: decoder, int, raw, classes

# ## RECONSTRUCT

# C = categorical(rand(UInt8,20,20))[1:2,1:3]
# @test length(levels(C)) > 6
# @test reconstruct(broadcast(identity, C)) == C

# C = categorical(string.(rand(UInt8,20,20)))[1:2,1:3]
# @test length(levels(C)) > 6
# @test reconstruct(broadcast(identity, C)) == C

## DECODER

N = 10
mix = shuffle(0:N - 1)

Xraw = broadcast(x->mod(x,N), rand(Int, 20, 100))
Yraw = string.(Xraw)

X = categorical(Xraw)
x = X[1]
Y = categorical(Yraw)
y = Y[1]
V = broadcast(identity, X)
W = broadcast(identity, Y)
Xo = levels!(deepcopy(X), mix)
xo = Xo[1]
Yo = levels!(deepcopy(Y), string.(mix))
yo = Yo[1]
Vo = broadcast(identity, Xo)
Wo = broadcast(identity, Yo)

# classes:
@test raw.(classes(xo)) == xo.pool.levels
@test raw.(classes(yo)) == yo.pool.levels

# getting all possible elements from one:
@test raw.(X) == Xraw
@test raw.(Y) == Yraw
@test raw.(classes(xo)) == levels(Xo)
@test raw.(classes(yo)) == levels(Yo)

# broadcasted encoding:
@test int(X) == int(V)
@test int(Y) == int(W)

# encoding and decoding are inverses:
d = decoder(xo)
@test d(int(Vo)) == Vo
e = decoder(yo)
@test e(int(Wo)) == Wo
A = sample(xo.pool.order, 100)
@test int(d(A)) == A
@test int(e(A)) == A

## MATRIX

B = rand(UInt8, (4, 5))
@test MLJBase.matrix(DataFrame(B)) == B

## ISTABLE

A = broadcast(x->Char(65+mod(x,5)), rand(Int, 10, 5))
X = CategoricalArrays.categorical(A)

df = DataFrame(A)
tt = Table(df)
# uncomment 4 lines to restore JuliaDB testing:
# db = JuliaDB.table(tt)
# nd = ndsparse((document=[6, 1, 1, 2, 3, 4],
#                word=[:house, :house, :sofa, :sofa, :chair, :house]), (values=["big", "small", 17, 34, 4, "small"],))

@test !MLJBase.istable(A)
@test !MLJBase.istable([1,2,3])
@test MLJBase.istable(df)
@test MLJBase.istable(tt)
# uncomment 2 lines to restore JuliaDB testing:
# @test MLJBase.istable(db)
# @test MLJBase.isndsparse(nd)

## TABLE INDEXING

df.z  = 1:10
tt = Table(df)
# uncomment 1 line to restore JuliaDB testing:
# db = JuliaDB.table(tt)

@test selectcols(df, 4:6) == selectcols(df[4:6], :)
@test selectcols(df, [:x1, :z]) == selectcols(df[[:x1, :z]], :)
@test selectcols(df, :x2) == df.x2
@test selectcols(df, 2) == df.x2
@test selectrows(df, 4:6) == selectrows(df[4:6, :], :)
@test selectrows(df, 1) == selectrows(df[1:1, :], :)
@test MLJBase.select(df, 2, :x2) == df[2,:x2]
s = schema(df)
@test s.names == tuple(names(df)...)
@test s.types == (Char, Char, Char, Char, Char, Int64)
@test nrows(df) == size(df, 1)

@test selectcols(tt, 4:6) == selectcols(Table(x4=tt.x4, x5=tt.x5, z=tt.z), :)
@test selectcols(tt, [:x1, :z]) == selectcols(Table(x1=tt.x1, z=tt.z), :)
@test selectcols(tt, :x2) == tt.x2
@test selectcols(tt, 2) == tt.x2
@test selectrows(tt, 4:6) == selectrows(tt[4:6], :)
s = schema(tt)
@test s.names == tuple(names(df)...)
@test s.types == (Char, Char, Char, Char, Char, Int64)
@test nrows(tt) == length(tt.x1)
@test MLJBase.select(tt, 2, :x2) == tt.x2[2]

v = rand(Int, 4)
@test selectrows(v, 2:3) == v[2:3]
@test nrows(v) == 4

v = categorical(collect("asdfasdf"))
@test selectrows(v, 2:3) == v[2:3]
@test nrows(v) == 8

df = DataFrame(v=v, w=v)
@test selectcols(df, :w) == v
tt = TypedTables.Table(df)
# uncomment 1 line to restore JuliaDB testing:
# db = JuliaDB.table(tt)
@test selectcols(tt, :w) == v
# uncomment 1 line to restore JuliaDB testing:
# @test selectcols(db, :w) == v

# uncomment 9 lines to restore JuliaDB testing:
# @test MLJBase.select(nd, 2, :house) isa Missing
# @test MLJBase.select(nd, 1, :house) == "small"
# @test all(MLJBase.select(nd, :, :house) .=== ["small", missing, missing, "small", missing, "big"])
# @test all(MLJBase.select(nd, [2,4], :house) .=== [missing, "small"])
# @test all(selectcols(nd, :house) .=== MLJBase.select(nd, :, :house))
# @test nrows(nd) == 6
# s = schema(nd)
# @test s.names == (:chair, :house, :sofa)
# @test s.types == (Union{Missing,Int64}, Union{Missing,String}, Union{Missing,Int64})

## MANIFESTING ARRAYS AS TABLES

A = hcat(v, v)
tab = MLJBase.table(A)
tab[1] == v
MLJBase.matrix(tab) == A

# uncomment 3 lines to restore JuliaDB testing:
# sparsearray = sparse([6, 1, 1, 2, 3, 4], [2, 2, 3, 3, 1, 2],
#                      ["big", "small", 17, 34, 4, "small"])
# @test MLJBase.matrix(nd) == sparsearray

end # module

true


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Symmetric and Hermitian matrices
struct Symmetric{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
    data::S
    uplo::Char

    function Symmetric{T,S}(data, uplo) where {T,S<:AbstractMatrix{<:T}}
        require_one_based_indexing(data)
        new{T,S}(data, uplo)
    end
end
"""
    Symmetric(A, uplo=:U)

Construct a `Symmetric` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`)
triangle of the matrix `A`.

# Examples
```jldoctest
julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4]
5×5 Array{Int64,2}:
 1  0  2  0  3
 0  4  0  5  0
 6  0  7  0  8
 0  9  0  1  0
 2  0  3  0  4

julia> Supper = Symmetric(A)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  2  0  3
 0  4  0  5  0
 2  0  7  0  8
 0  5  0  1  0
 3  0  8  0  4

julia> Slower = Symmetric(A, :L)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  6  0  2
 0  4  0  9  0
 6  0  7  0  3
 0  9  0  1  0
 2  0  3  0  4
```

Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if `A == transpose(A)`).
"""
function Symmetric(A::AbstractMatrix, uplo::Symbol=:U)
    checksquare(A)
    return symmetric_type(typeof(A))(A, char_uplo(uplo))
end

"""
    symmetric(A, uplo=:U)

Construct a symmetric view of `A`. If `A` is a matrix, `uplo` controls whether the upper
(if `uplo = :U`) or lower (if `uplo = :L`) triangle of `A` is used to implicitly fill the
other one. If `A` is a `Number`, it is returned as is.

If a symmetric view of a matrix is to be constructed of which the elements are neither
matrices nor numbers, an appropriate method of `symmetric` has to be implemented. In that
case, `symmetric_type` has to be implemented, too.
"""
symmetric(A::AbstractMatrix, uplo::Symbol) = Symmetric(A, uplo)
symmetric(A::Number, ::Symbol) = A

"""
    symmetric_type(T::Type)

The type of the object returned by `symmetric(::T, ::Symbol)`. For matrices, this is an
appropriately typed `Symmetric`, for `Number`s, it is the original type. If `symmetric` is
implemented for a custom type, so should be `symmetric_type`, and vice versa.
"""
function symmetric_type(::Type{T}) where {S, T<:AbstractMatrix{S}}
    return Symmetric{Union{S, promote_op(transpose, S), symmetric_type(S)}, T}
end
function symmetric_type(::Type{T}) where {S<:Number, T<:AbstractMatrix{S}}
    return Symmetric{S, T}
end
function symmetric_type(::Type{T}) where {S<:AbstractMatrix, T<:AbstractMatrix{S}}
    return Symmetric{AbstractMatrix, T}
end
symmetric_type(::Type{T}) where {T<:Number} = T

struct Hermitian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
    data::S
    uplo::Char

    function Hermitian{T,S}(data, uplo) where {T,S<:AbstractMatrix{<:T}}
        require_one_based_indexing(data)
        new{T,S}(data, uplo)
    end
end
"""
    Hermitian(A, uplo=:U)

Construct a `Hermitian` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`)
triangle of the matrix `A`.

# Examples
```jldoctest
julia> A = [1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4];

julia> Hupper = Hermitian(A)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  2+2im  0+0im  3-3im
 0+0im  4+0im  0+0im  5+0im  0+0im
 2-2im  0+0im  7+0im  0+0im  8+8im
 0+0im  5+0im  0+0im  1+0im  0+0im
 3+3im  0+0im  8-8im  0+0im  4+0im

julia> Hlower = Hermitian(A, :L)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  6+6im  0+0im  2-2im
 0+0im  4+0im  0+0im  9+0im  0+0im
 6-6im  0+0im  7+0im  0+0im  3+3im
 0+0im  9+0im  0+0im  1+0im  0+0im
 2+2im  0+0im  3-3im  0+0im  4+0im
```

Note that `Hupper` will not be equal to `Hlower` unless `A` is itself Hermitian (e.g. if `A == adjoint(A)`).

All non-real parts of the diagonal will be ignored.

```julia
Hermitian(fill(complex(1,1), 1, 1)) == fill(1, 1, 1)
```
"""
function Hermitian(A::AbstractMatrix, uplo::Symbol=:U)
    n = checksquare(A)
    return hermitian_type(typeof(A))(A, char_uplo(uplo))
end

"""
    hermitian(A, uplo=:U)

Construct a hermitian view of `A`. If `A` is a matrix, `uplo` controls whether the upper
(if `uplo = :U`) or lower (if `uplo = :L`) triangle of `A` is used to implicitly fill the
other one. If `A` is a `Number`, its real part is returned converted back to the input
type.

If a hermitian view of a matrix is to be constructed of which the elements are neither
matrices nor numbers, an appropriate method of `hermitian` has to be implemented. In that
case, `hermitian_type` has to be implemented, too.
"""
hermitian(A::AbstractMatrix, uplo::Symbol) = Hermitian(A, uplo)
hermitian(A::Number, ::Symbol) = convert(typeof(A), real(A))

"""
    hermitian_type(T::Type)

The type of the object returned by `hermitian(::T, ::Symbol)`. For matrices, this is an
appropriately typed `Hermitian`, for `Number`s, it is the original type. If `hermitian` is
implemented for a custom type, so should be `hermitian_type`, and vice versa.
"""
function hermitian_type(::Type{T}) where {S, T<:AbstractMatrix{S}}
    return Hermitian{Union{S, promote_op(adjoint, S), hermitian_type(S)}, T}
end
function hermitian_type(::Type{T}) where {S<:Number, T<:AbstractMatrix{S}}
    return Hermitian{S, T}
end
function hermitian_type(::Type{T}) where {S<:AbstractMatrix, T<:AbstractMatrix{S}}
    return Hermitian{AbstractMatrix, T}
end
hermitian_type(::Type{T}) where {T<:Number} = T

for (S, H) in ((:Symmetric, :Hermitian), (:Hermitian, :Symmetric))
    @eval begin
        $S(A::$S) = A
        function $S(A::$S, uplo::Symbol)
            if A.uplo == char_uplo(uplo)
                return A
            else
                throw(ArgumentError("Cannot construct $($S); uplo doesn't match"))
            end
        end
        $S(A::$H) = $S(A, sym_uplo(A.uplo))
        function $S(A::$H, uplo::Symbol)
            if A.uplo == char_uplo(uplo)
                if $H === Hermitian && !(eltype(A) <: Real) &&
                    any(!isreal, A.data[i] for i in diagind(A.data))

                    throw(ArgumentError("Cannot construct $($S)($($H))); diagonal contains complex values"))
                end
                return $S(A.data, sym_uplo(A.uplo))
            else
                throw(ArgumentError("Cannot construct $($S); uplo doesn't match"))
            end
        end
    end
end

convert(T::Type{<:Symmetric}, m::Union{Symmetric,Hermitian}) = m isa T ? m : T(m)
convert(T::Type{<:Hermitian}, m::Union{Symmetric,Hermitian}) = m isa T ? m : T(m)

const HermOrSym{T,        S} = Union{Hermitian{T,S}, Symmetric{T,S}}
const RealHermSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}}
const RealHermSymComplexHerm{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}
const RealHermSymComplexSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Symmetric{Complex{T},S}}

size(A::HermOrSym, d) = size(A.data, d)
size(A::HermOrSym) = size(A.data)
@inline function getindex(A::Symmetric, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    @inbounds if i == j
        return symmetric(A.data[i, j], sym_uplo(A.uplo))::symmetric_type(eltype(A.data))
    elseif (A.uplo == 'U') == (i < j)
        return A.data[i, j]
    else
        return transpose(A.data[j, i])
    end
end
@inline function getindex(A::Hermitian, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    @inbounds if i == j
        return hermitian(A.data[i, j], sym_uplo(A.uplo))::hermitian_type(eltype(A.data))
    elseif (A.uplo == 'U') == (i < j)
        return A.data[i, j]
    else
        return adjoint(A.data[j, i])
    end
end

function setindex!(A::Symmetric, v, i::Integer, j::Integer)
    i == j || throw(ArgumentError("Cannot set a non-diagonal index in a symmetric matrix"))
    setindex!(A.data, v, i, j)
end

function setindex!(A::Hermitian, v, i::Integer, j::Integer)
    if i != j
        throw(ArgumentError("Cannot set a non-diagonal index in a Hermitian matrix"))
    elseif !isreal(v)
        throw(ArgumentError("Cannot set a diagonal entry in a Hermitian matrix to a nonreal value"))
    else
        setindex!(A.data, v, i, j)
    end
end

# For A<:Union{Symmetric,Hermitian}, similar(A[, neweltype]) should yield a matrix with the same
# symmetry type, uplo flag, and underlying storage type as A. The following methods cover these cases.
similar(A::Symmetric, ::Type{T}) where {T} = Symmetric(similar(parent(A), T), ifelse(A.uplo == 'U', :U, :L))
# If the the Hermitian constructor's check ascertaining that the wrapped matrix's
# diagonal is strictly real is removed, the following method can be simplified.
function similar(A::Hermitian, ::Type{T}) where T
    B = similar(parent(A), T)
    for i in 1:size(B, 1) B[i, i] = 0 end
    return Hermitian(B, ifelse(A.uplo == 'U', :U, :L))
end
# On the other hand, similar(A, [neweltype,] shape...) should yield a matrix of the underlying
# storage type of A (not wrapped in a symmetry type). The following method covers these cases.
similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = similar(parent(A), T, dims)

# Conversion
function Matrix(A::Symmetric)
    B = copytri!(convert(Matrix, copy(A.data)), A.uplo)
    for i = 1:size(A, 1)
        B[i,i] = symmetric(A[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data))
    end
    return B
end
function Matrix(A::Hermitian)
    B = copytri!(convert(Matrix, copy(A.data)), A.uplo, true)
    for i = 1:size(A, 1)
        B[i,i] = hermitian(A[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data))
    end
    return B
end
Array(A::Union{Symmetric,Hermitian}) = convert(Matrix, A)

parent(A::HermOrSym) = A.data
Symmetric{T,S}(A::Symmetric{T,S}) where {T,S<:AbstractMatrix} = A
Symmetric{T,S}(A::Symmetric) where {T,S<:AbstractMatrix} = Symmetric{T,S}(convert(S,A.data),A.uplo)
AbstractMatrix{T}(A::Symmetric) where {T} = Symmetric(convert(AbstractMatrix{T}, A.data), sym_uplo(A.uplo))
Hermitian{T,S}(A::Hermitian{T,S}) where {T,S<:AbstractMatrix} = A
Hermitian{T,S}(A::Hermitian) where {T,S<:AbstractMatrix} = Hermitian{T,S}(convert(S,A.data),A.uplo)
AbstractMatrix{T}(A::Hermitian) where {T} = Hermitian(convert(AbstractMatrix{T}, A.data), sym_uplo(A.uplo))

copy(A::Symmetric{T,S}) where {T,S} = (B = copy(A.data); Symmetric{T,typeof(B)}(B,A.uplo))
copy(A::Hermitian{T,S}) where {T,S} = (B = copy(A.data); Hermitian{T,typeof(B)}(B,A.uplo))

function copyto!(dest::Symmetric, src::Symmetric)
    if src.uplo == dest.uplo
        copyto!(dest.data, src.data)
    else
        transpose!(dest.data, src.data)
    end
    return dest
end

function copyto!(dest::Hermitian, src::Hermitian)
    if src.uplo == dest.uplo
        copyto!(dest.data, src.data)
    else
        adjoint!(dest.data, src.data)
    end
    return dest
end

# fill[stored]!
fill!(A::HermOrSym, x) = fillstored!(A, x)
function fillstored!(A::HermOrSym{T}, x) where T
    xT = convert(T, x)
    if isa(A, Hermitian)
        isreal(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a nonreal value"))
    end
    if A.uplo == 'U'
        fillband!(A.data, xT, 0, size(A,2)-1)
    else # A.uplo == 'L'
        fillband!(A.data, xT, 1-size(A,1), 0)
    end
    return A
end

function Base.isreal(A::HermOrSym)
    n = size(A, 1)
    @inbounds if A.uplo == 'U'
        for j in 1:n
            for i in 1:(j - (A isa Hermitian))
                if !isreal(A.data[i,j])
                    return false
                end
            end
        end
    else
        for j in 1:n
            for i in (j + (A isa Hermitian)):n
                if !isreal(A.data[i,j])
                    return false
                end
            end
        end
    end
    return true
end

ishermitian(A::Hermitian) = true
ishermitian(A::Symmetric{<:Real}) = true
ishermitian(A::Symmetric{<:Complex}) = isreal(A)
issymmetric(A::Hermitian{<:Real}) = true
issymmetric(A::Hermitian{<:Complex}) = isreal(A)
issymmetric(A::Symmetric) = true

adjoint(A::Hermitian) = A
transpose(A::Symmetric) = A
adjoint(A::Symmetric{<:Real}) = A
transpose(A::Hermitian{<:Real}) = A
adjoint(A::Symmetric) = Adjoint(A)
transpose(A::Hermitian) = Transpose(A)

real(A::Symmetric{<:Real}) = A
real(A::Hermitian{<:Real}) = A
real(A::Symmetric) = Symmetric(real(A.data), sym_uplo(A.uplo))
real(A::Hermitian) = Hermitian(real(A.data), sym_uplo(A.uplo))
imag(A::Symmetric) = Symmetric(imag(A.data), sym_uplo(A.uplo))

Base.copy(A::Adjoint{<:Any,<:Hermitian}) = copy(A.parent)
Base.copy(A::Transpose{<:Any,<:Symmetric}) = copy(A.parent)
Base.copy(A::Adjoint{<:Any,<:Symmetric}) =
    Symmetric(copy(adjoint(A.parent.data)), ifelse(A.parent.uplo == 'U', :L, :U))
Base.copy(A::Transpose{<:Any,<:Hermitian}) =
    Hermitian(copy(transpose(A.parent.data)), ifelse(A.parent.uplo == 'U', :L, :U))

tr(A::Hermitian) = real(tr(A.data))

Base.conj(A::HermOrSym) = typeof(A)(conj(A.data), A.uplo)
Base.conj!(A::HermOrSym) = typeof(A)(conj!(A.data), A.uplo)

# tril/triu
function tril(A::Hermitian, k::Integer=0)
    if A.uplo == 'U' && k <= 0
        return tril!(copy(A.data'),k)
    elseif A.uplo == 'U' && k > 0
        return tril!(copy(A.data'),-1) + tril!(triu(A.data),k)
    elseif A.uplo == 'L' && k <= 0
        return tril(A.data,k)
    else
        return tril(A.data,-1) + tril!(triu!(copy(A.data')),k)
    end
end

function tril(A::Symmetric, k::Integer=0)
    if A.uplo == 'U' && k <= 0
        return tril!(copy(transpose(A.data)),k)
    elseif A.uplo == 'U' && k > 0
        return tril!(copy(transpose(A.data)),-1) + tril!(triu(A.data),k)
    elseif A.uplo == 'L' && k <= 0
        return tril(A.data,k)
    else
        return tril(A.data,-1) + tril!(triu!(copy(transpose(A.data))),k)
    end
end

function triu(A::Hermitian, k::Integer=0)
    if A.uplo == 'U' && k >= 0
        return triu(A.data,k)
    elseif A.uplo == 'U' && k < 0
        return triu(A.data,1) + triu!(tril!(copy(A.data')),k)
    elseif A.uplo == 'L' && k >= 0
        return triu!(copy(A.data'),k)
    else
        return triu!(copy(A.data'),1) + triu!(tril(A.data),k)
    end
end

function triu(A::Symmetric, k::Integer=0)
    if A.uplo == 'U' && k >= 0
        return triu(A.data,k)
    elseif A.uplo == 'U' && k < 0
        return triu(A.data,1) + triu!(tril!(copy(transpose(A.data))),k)
    elseif A.uplo == 'L' && k >= 0
        return triu!(copy(transpose(A.data)),k)
    else
        return triu!(copy(transpose(A.data)),1) + triu!(tril(A.data),k)
    end
end

function dot(A::Symmetric, B::Symmetric)
    n = size(A, 2)
    if n != size(B, 2)
        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
    end

    dotprod = zero(dot(first(A), first(B)))
    @inbounds if A.uplo == 'U' && B.uplo == 'U'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    elseif A.uplo == 'L' && B.uplo == 'L'
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
            end
        end
    elseif A.uplo == 'U' && B.uplo == 'L'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i]))
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    else
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i]))
            end
        end
    end
    return dotprod
end

function dot(A::Hermitian, B::Hermitian)
    n = size(A, 2)
    if n != size(B, 2)
        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
    end

    dotprod = zero(dot(first(A), first(B)))
    @inbounds if A.uplo == 'U' && B.uplo == 'U'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * real(dot(A.data[i, j], B.data[i, j]))
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    elseif A.uplo == 'L' && B.uplo == 'L'
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * real(dot(A.data[i, j], B.data[i, j]))
            end
        end
    elseif A.uplo == 'U' && B.uplo == 'L'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i])))
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    else
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i])))
            end
        end
    end
    return dotprod
end

(-)(A::Symmetric) = Symmetric(-A.data, sym_uplo(A.uplo))
(-)(A::Hermitian) = Hermitian(-A.data, sym_uplo(A.uplo))

## Addition/subtraction
for f in (:+, :-)
    @eval $f(A::Symmetric, B::Symmetric) = Symmetric($f(A.data, B), sym_uplo(A.uplo))
    @eval $f(A::Hermitian, B::Hermitian) = Hermitian($f(A.data, B), sym_uplo(A.uplo))
    @eval $f(A::Hermitian, B::Symmetric{<:Real}) = Hermitian($f(A.data, B), sym_uplo(A.uplo))
    @eval $f(A::Symmetric{<:Real}, B::Hermitian) = Hermitian($f(A.data, B), sym_uplo(A.uplo))
end

## Matvec
@inline mul!(y::StridedVector{T}, A::Symmetric{T,<:StridedMatrix}, x::StridedVector{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasFloat} =
    BLAS.symv!(A.uplo, alpha, A.data, x, beta, y)
@inline mul!(y::StridedVector{T}, A::Hermitian{T,<:StridedMatrix}, x::StridedVector{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasReal} =
    BLAS.symv!(A.uplo, alpha, A.data, x, beta, y)
@inline mul!(y::StridedVector{T}, A::Hermitian{T,<:StridedMatrix}, x::StridedVector{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasComplex} =
    BLAS.hemv!(A.uplo, alpha, A.data, x, beta, y)
## Matmat
@inline mul!(C::StridedMatrix{T}, A::Symmetric{T,<:StridedMatrix}, B::StridedMatrix{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasFloat} =
    BLAS.symm!('L', A.uplo, alpha, A.data, B, beta, C)
@inline mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Symmetric{T,<:StridedMatrix},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasFloat} =
    BLAS.symm!('R', B.uplo, alpha, B.data, A, beta, C)
@inline mul!(C::StridedMatrix{T}, A::Hermitian{T,<:StridedMatrix}, B::StridedMatrix{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasReal} =
    BLAS.symm!('L', A.uplo, alpha, A.data, B, beta, C)
@inline mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Hermitian{T,<:StridedMatrix},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasReal} =
    BLAS.symm!('R', B.uplo, alpha, B.data, A, beta, C)
@inline mul!(C::StridedMatrix{T}, A::Hermitian{T,<:StridedMatrix}, B::StridedMatrix{T},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasComplex} =
    BLAS.hemm!('L', A.uplo, alpha, A.data, B, beta, C)
@inline mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Hermitian{T,<:StridedMatrix},
             alpha::Union{T, Bool}, beta::Union{T, Bool}) where {T<:BlasComplex} =
    BLAS.hemm!('R', B.uplo, alpha, B.data, A, beta, C)

*(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)

function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
    require_one_based_indexing(x, y)
    (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
    data = A.data
    r = zero(eltype(x)) * zero(eltype(A)) * zero(eltype(y))
    if A.uplo == 'U'
        @inbounds for j = 1:length(y)
            r += dot(x[j], real(data[j,j]), y[j])
            @simd for i = 1:j-1
                Aij = data[i,j]
                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
            end
        end
    else # A.uplo == 'L'
        @inbounds for j = 1:length(y)
            r += dot(x[j], real(data[j,j]), y[j])
            @simd for i = j+1:length(y)
                Aij = data[i,j]
                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
            end
        end
    end
    return r
end

# Fallbacks to avoid generic_matvecmul!/generic_matmatmul!
## Symmetric{<:Number} and Hermitian{<:Real} are invariant to transpose; peel off the t
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractVector) = transA.parent * B
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractMatrix) = transA.parent * B
*(A::AbstractMatrix, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = A * transB.parent
## Hermitian{<:Number} and Symmetric{<:Real} are invariant to adjoint; peel off the c
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractVector) = adjA.parent * B
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractMatrix) = adjA.parent * B
*(A::AbstractMatrix, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * adjB.parent

# ambiguities with transposed AbstractMatrix methods in linalg/matmul.jl
*(transA::Transpose{<:Any,<:RealHermSym}, transB::Transpose{<:Any,<:RealHermSym}) = transA * transB.parent
*(transA::Transpose{<:Any,<:RealHermSym}, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = transA * transB.parent
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = transA.parent * transB.parent
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transB::Transpose{<:Any,<:RealHermSym}) = transA.parent * transB
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transB::Transpose{<:Any,<:RealHermSymComplexHerm}) = transA.parent * transB
*(transA::Transpose{<:Any,<:RealHermSymComplexHerm}, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = transA * transB.parent
*(adjA::Adjoint{<:Any,<:RealHermSym}, adjB::Adjoint{<:Any,<:RealHermSym}) = adjA * adjB.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjA.parent * adjB.parent
*(adjA::Adjoint{<:Any,<:RealHermSym}, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjA * adjB.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexSym}, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjA * adjB.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjB::Adjoint{<:Any,<:RealHermSym}) = adjA.parent * adjB
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjB::Adjoint{<:Any,<:RealHermSymComplexSym}) = adjA.parent * adjB

# ambiguities with AbstractTriangular
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractTriangular) = transA.parent * B
*(A::AbstractTriangular, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = A * transB.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractTriangular) = adjA.parent * B
*(A::AbstractTriangular, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * adjB.parent

# Scaling with Number
*(A::Symmetric, x::Number) = Symmetric(A.data*x, sym_uplo(A.uplo))
*(x::Number, A::Symmetric) = Symmetric(x*A.data, sym_uplo(A.uplo))
*(A::Hermitian, x::Real) = Hermitian(A.data*x, sym_uplo(A.uplo))
*(x::Real, A::Hermitian) = Hermitian(x*A.data, sym_uplo(A.uplo))
/(A::Symmetric, x::Number) = Symmetric(A.data/x, sym_uplo(A.uplo))
/(A::Hermitian, x::Real) = Hermitian(A.data/x, sym_uplo(A.uplo))

function factorize(A::HermOrSym{T}) where T
    TT = typeof(sqrt(oneunit(T)))
    if TT <: BlasFloat
        return bunchkaufman(A)
    else # fallback
        return lu(A)
    end
end

det(A::RealHermSymComplexHerm) = real(det(factorize(A)))
det(A::Symmetric{<:Real}) = det(factorize(A))
det(A::Symmetric) = det(factorize(A))

\(A::HermOrSym{<:Any,<:StridedMatrix}, B::AbstractVector) = \(factorize(A), B)
# Bunch-Kaufman solves can not utilize BLAS-3 for multiple right hand sides
# so using LU is faster for AbstractMatrix right hand side
\(A::HermOrSym{<:Any,<:StridedMatrix}, B::AbstractMatrix) = \(lu(A), B)

function _inv(A::HermOrSym)
    n = checksquare(A)
    B = inv!(lu(A))
    conjugate = isa(A, Hermitian)
    # symmetrize
    if A.uplo == 'U' # add to upper triangle
        @inbounds for i = 1:n, j = i:n
            B[i,j] = conjugate ? (B[i,j] + conj(B[j,i])) / 2 : (B[i,j] + B[j,i]) / 2
        end
    else # A.uplo == 'L', add to lower triangle
        @inbounds for i = 1:n, j = i:n
            B[j,i] = conjugate ? (B[j,i] + conj(B[i,j])) / 2 : (B[j,i] + B[i,j]) / 2
        end
    end
    B
end
inv(A::Hermitian{<:Any,<:StridedMatrix}) = Hermitian(_inv(A), sym_uplo(A.uplo))
inv(A::Symmetric{<:Any,<:StridedMatrix}) = Symmetric(_inv(A), sym_uplo(A.uplo))

eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) = Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)

function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
    T = eltype(A)
    S = eigtype(T)
    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby=sortby)
end

eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)

"""
    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen

Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).

The `UnitRange` `irange` specifies indices of the sorted eigenvalues to search for.

!!! note
    If `irange` is not `1:n`, where `n` is the dimension of `A`, then the returned factorization
    will be a *truncated* factorization.
"""
function eigen(A::RealHermSymComplexHerm, irange::UnitRange)
    T = eltype(A)
    S = eigtype(T)
    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
end

eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where {T<:BlasReal} =
    Eigen(LAPACK.syevr!('V', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)...)

"""
    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen

Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)

Iterating the decomposition produces the components `F.values` and `F.vectors`.

The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).

`vl` is the lower bound of the window of eigenvalues to search for, and `vu` is the upper bound.

!!! note
    If [`vl`, `vu`] does not contain all eigenvalues of `A`, then the returned factorization
    will be a *truncated* factorization.
"""
function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
    T = eltype(A)
    S = eigtype(T)
    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
end

eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) =
    LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]

function eigvals(A::RealHermSymComplexHerm)
    T = eltype(A)
    S = eigtype(T)
    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A))
end

"""
    eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values

Same as [`eigvals`](@ref), but saves space by overwriting the input `A`, instead of creating a copy.
`irange` is a range of eigenvalue *indices* to search for - for instance, the 2nd to 8th eigenvalues.
"""
eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) =
    LAPACK.syevr!('N', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)[1]

"""
    eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values

Returns the eigenvalues of `A`. It is possible to calculate only a subset of the
eigenvalues by specifying a `UnitRange` `irange` covering indices of the sorted eigenvalues,
e.g. the 2nd to 8th eigenvalues.

```jldoctest
julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0  2.0   ⋅
 2.0  2.0  3.0
  ⋅   3.0  1.0

julia> eigvals(A, 2:2)
1-element Array{Float64,1}:
 0.9999999999999996

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.1400549446402604
  1.0000000000000002
  5.140054944640259
```
"""
function eigvals(A::RealHermSymComplexHerm, irange::UnitRange)
    T = eltype(A)
    S = eigtype(T)
    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
end

"""
    eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> values

Same as [`eigvals`](@ref), but saves space by overwriting the input `A`, instead of creating a copy.
`vl` is the lower bound of the interval to search for eigenvalues, and `vu` is the upper bound.
"""
eigvals!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where {T<:BlasReal} =
    LAPACK.syevr!('N', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)[1]

"""
    eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> values

Returns the eigenvalues of `A`. It is possible to calculate only a subset of the eigenvalues
by specifying a pair `vl` and `vu` for the lower and upper boundaries of the eigenvalues.

```jldoctest
julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0  2.0   ⋅
 2.0  2.0  3.0
  ⋅   3.0  1.0

julia> eigvals(A, -1, 2)
1-element Array{Float64,1}:
 1.0000000000000009

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.1400549446402604
  1.0000000000000002
  5.140054944640259
```
"""
function eigvals(A::RealHermSymComplexHerm, vl::Real, vh::Real)
    T = eltype(A)
    S = eigtype(T)
    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
end

eigmax(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, size(A, 1):size(A, 1))[1]
eigmin(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, 1:1)[1]

function eigen!(A::HermOrSym{T,S}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
end
function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasComplex,S<:StridedMatrix}
    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
    GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
end

eigvals!(A::HermOrSym{T,S}, B::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatrix} =
    LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))[1]
eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}) where {T<:BlasComplex,S<:StridedMatrix} =
    LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))[1]

eigvecs(A::HermOrSym) = eigvecs(eigen(A))

function svd(A::RealHermSymComplexHerm, full::Bool=false)
    vals, vecs = eigen(A)
    I = sortperm(vals; by=abs, rev=true)
    permute!(vals, I)
    Base.permutecols!!(vecs, I)         # left-singular vectors
    V = copy(vecs)                      # right-singular vectors
    # shifting -1 from singular values to right-singular vectors
    @inbounds for i = 1:length(vals)
        if vals[i] < 0
            vals[i] = -vals[i]
            for j = 1:size(V,1); V[j,i] = -V[j,i]; end
        end
    end
    return SVD(vecs, vals, V')
end

function svdvals!(A::RealHermSymComplexHerm)
    vals = eigvals!(A)
    for i = 1:length(vals)
        vals[i] = abs(vals[i])
    end
    return sort!(vals, rev = true)
end

# Matrix functions
^(A::Symmetric{<:Real}, p::Integer) = sympow(A, p)
^(A::Symmetric{<:Complex}, p::Integer) = sympow(A, p)
function sympow(A::Symmetric, p::Integer)
    if p < 0
        return Symmetric(Base.power_by_squaring(inv(A), -p))
    else
        return Symmetric(Base.power_by_squaring(A, p))
    end
end
function ^(A::Symmetric{<:Real}, p::Real)
    isinteger(p) && return integerpow(A, p)
    F = eigen(A)
    if all(λ -> λ ≥ 0, F.values)
        return Symmetric((F.vectors * Diagonal((F.values).^p)) * F.vectors')
    else
        return Symmetric((F.vectors * Diagonal((complex(F.values)).^p)) * F.vectors')
    end
end
function ^(A::Symmetric{<:Complex}, p::Real)
    isinteger(p) && return integerpow(A, p)
    return Symmetric(schurpow(A, p))
end
function ^(A::Hermitian, p::Integer)
    if p < 0
        retmat = Base.power_by_squaring(inv(A), -p)
    else
        retmat = Base.power_by_squaring(A, p)
    end
    for i = 1:size(A,1)
        retmat[i,i] = real(retmat[i,i])
    end
    return Hermitian(retmat)
end
function ^(A::Hermitian{T}, p::Real) where T
    isinteger(p) && return integerpow(A, p)
    F = eigen(A)
    if all(λ -> λ ≥ 0, F.values)
        retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors'
        if T <: Real
            return Hermitian(retmat)
        else
            for i = 1:size(A,1)
                retmat[i,i] = real(retmat[i,i])
            end
            return Hermitian(retmat)
        end
    else
        return (F.vectors * Diagonal((complex(F.values).^p))) * F.vectors'
    end
end

for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh)
    @eval begin
        function ($func)(A::HermOrSym{<:Real})
            F = eigen(A)
            return Symmetric((F.vectors * Diagonal(($func).(F.values))) * F.vectors')
        end
        function ($func)(A::Hermitian{<:Complex})
            n = checksquare(A)
            F = eigen(A)
            retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
            for i = 1:n
                retmat[i,i] = real(retmat[i,i])
            end
            return Hermitian(retmat)
        end
    end
end

for func in (:acos, :asin)
    @eval begin
        function ($func)(A::HermOrSym{<:Real})
            F = eigen(A)
            if all(λ -> -1 ≤ λ ≤ 1, F.values)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
            else
                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
            end
            return Symmetric(retmat)
        end
        function ($func)(A::Hermitian{<:Complex})
            n = checksquare(A)
            F = eigen(A)
            if all(λ -> -1 ≤ λ ≤ 1, F.values)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
                for i = 1:n
                    retmat[i,i] = real(retmat[i,i])
                end
                return Hermitian(retmat)
            else
                return (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
            end
        end
    end
end

function acosh(A::HermOrSym{<:Real})
    F = eigen(A)
    if all(λ -> λ ≥ 1, F.values)
        retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors'
    else
        retmat = (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors'
    end
    return Symmetric(retmat)
end
function acosh(A::Hermitian{<:Complex})
    n = checksquare(A)
    F = eigen(A)
    if all(λ -> λ ≥ 1, F.values)
        retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors'
        for i = 1:n
            retmat[i,i] = real(retmat[i,i])
        end
        return Hermitian(retmat)
    else
        return (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors'
    end
end

function sincos(A::HermOrSym{<:Real})
    n = checksquare(A)
    F = eigen(A)
    S, C = Diagonal(similar(A, (n,))), Diagonal(similar(A, (n,)))
    for i in 1:n
        S.diag[i], C.diag[i] = sincos(F.values[i])
    end
    return Symmetric((F.vectors * S) * F.vectors'), Symmetric((F.vectors * C) * F.vectors')
end
function sincos(A::Hermitian{<:Complex})
    n = checksquare(A)
    F = eigen(A)
    S, C = Diagonal(similar(A, (n,))), Diagonal(similar(A, (n,)))
    for i in 1:n
        S.diag[i], C.diag[i] = sincos(F.values[i])
    end
    retmatS, retmatC = (F.vectors * S) * F.vectors', (F.vectors * C) * F.vectors'
    for i = 1:n
        retmatS[i,i] = real(retmatS[i,i])
        retmatC[i,i] = real(retmatC[i,i])
    end
    return Hermitian(retmatS), Hermitian(retmatC)
end

for func in (:log, :sqrt)
    @eval begin
        function ($func)(A::HermOrSym{<:Real})
            F = eigen(A)
            if all(λ -> λ ≥ 0, F.values)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
            else
                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
            end
            return Symmetric(retmat)
        end

        function ($func)(A::Hermitian{<:Complex})
            n = checksquare(A)
            F = eigen(A)
            if all(λ -> λ ≥ 0, F.values)
                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
                for i = 1:n
                    retmat[i,i] = real(retmat[i,i])
                end
                return Hermitian(retmat)
            else
                retmat = (F.vectors * Diagonal(($func).(complex(F.values)))) * F.vectors'
                return retmat
            end
        end
    end
end

# disambiguation methods: *(Adj of RealHermSymComplexHerm, Trans of RealHermSymComplexSym) and symmetric partner
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A.parent * B.parent
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A.parent * B.parent
# disambiguation methods: *(Adj/Trans of AbsVec/AbsMat, Adj/Trans of RealHermSymComplex{Herm|Sym})
*(A::Adjoint{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
# disambiguation methods: *(Adj/Trans of RealHermSymComplex{Herm|Sym}, Adj/Trans of AbsVec/AbsMat)
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractVector}) = A.parent * B
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractMatrix}) = A.parent * B
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractVector}) = A.parent * B
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractMatrix}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractVector}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractMatrix}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractVector}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractMatrix}) = A.parent * B

# disambiguation methods: *(Adj/Trans of AbsTri or RealHermSymComplex{Herm|Sym}, Adj/Trans of other)
*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent
*(A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractTriangular}) = A.parent * B
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractTriangular}) = A.parent * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B

max_stars_repo_name: smallnamespace/julia

max_issues_repo_name: smallnamespace/julia

max_forks_repo_name: smallnamespace/julia


text:
let

midi = readMIDIFile("serenade_full.mid")
piano = midi.tracks[4]
notes = getnotes(piano, midi.tpq)

@testset "MIDI IO" begin
    @test midi.tpq == 960
    @test length(midi.tracks) == 4

    notes = getnotes(midi.tracks[2])

    @test length(notes.notes) > 1
    @test notes.tpq == 960

    @testset "extract numbers" begin
        for f in (velocities, positions, pitches, durations)
            @test length(f(notes)) == length(notes)
        end
    end
end

@testset "Replace Notes" begin

    midi = readMIDIFile("serenade_full.mid")
    notes2 = getnotes(midi.tracks[2])
    newtrack = deepcopy(midi.tracks[2])
    newnotes = getnotes(newtrack)

    for f in (velocities, positions, pitches, durations)
        @test f(notes2) == f(newnotes)
    end

    for note in newnotes
        note.velocity = rand(1:126)
    end

    newtrack = replace_notes(newtrack, newnotes)
    newnotes = getnotes(newtrack)

    for f in (positions, pitches, durations)
        @test f(notes2) == f(newnotes)
    end

    @test velocities(notes) !== velocities(newnotes)
end
end


text:
# IDEA: GradientLaplacianKernel?
struct HessianKernel{T, K} <: MultiKernel{T}
    k::K
    # input_type::IT  IT<:InputTrait
    # could A temporary storage for gradient calculation here
end
input_trait(G::HessianKernel) = input_trait(G.k)
HessianKernel(k::AbstractKernel{T}) where {T} = HessianKernel{T, typeof(k)}(k)
HessianKernel(k) = HessianKernel{Float64, typeof(k)}(k) # use fieldtype here?

function elsize(G::Gramian{<:AbstractMatrix, <:HessianKernel}, i::Int)
    i ≤ 2 ? length(G.x[1])^2 : 1
end

# IDEA: input_trait(G.k) could be pre-computed
function (G::HessianKernel)(x::AbstractVector, y::AbstractVector)
    hessian_kernel(G.k, x, y, input_trait(G.k))
end

function hessian_kernel(k, x::AbstractVector, y::AbstractVector, T::InputTrait)#::GenericInput = GenericInput())
    K = allocate_hessian_kernel(k, x, y, T)
    hessian_kernel!(K, k, x, y, T)
end

########################## Generic Hessian Kernel ##############################
# allocates space for gradient kernel evaluation but does not evaluate
# separation from evaluation useful for ValueHessianKernel
function allocate_hessian_kernel(k, x, y, ::InputTrait)
    d = length(x)
    zeros(d^2, d^2)
end

function hessian_kernel!(K::AbstractMatrix, k, x::AbstractVector, y::AbstractVector, ::InputTrait)
    d = length(x)
    value = similar(x, d^2)
    derivs = K # Hessian of Hessian (higher order terms)
    result = DiffResults.DiffResult(value, derivs)
    g(y) = vec(ForwardDiff.hessian(z->k(z, y), x)) # covariance between vectorized Hessian and value
    ForwardDiff.jacobian!(result, z->ForwardDiff.jacobian(g, z), y) # hessian of vector valued function
    return K
end

############################ Lazy Kernel Element ###############################
# struct HessianKernelElement{T, K, X<:AbstractVector{T}, Y} <: Factorization{T}
#     k::K
#     x::X
#     y::Y
#     # input_trait::IT
# end
#
# # HessianKernelElement(k, x, y) = HessianKernelElement(k, x, y, input_trait(k))
# function Base.size(K::HessianKernelElement)
#     d = length(K.x)
#     (d^2, d^2)
# end
# Base.size(K::HessianKernelElement, i::Int) = 0 < i ≤ 2 ? size(K)[i] : 1

# function allocate_hessian_kernel(k, x, y, T::Union{IsotropicInput, DotProductInput})
#     HessianKernelElement(k, copy(x), copy(y))
# end

# function hessian_kernel!(K::HessianKernelElement, k, x::AbstractVector, y::AbstractVector, T::Union{IsotropicInput, DotProductInput})
#     K.x .= x; K.y .= y
#     return K
# end

# function evaluate_block!(K::HessianKernelElement, G::HessianKernel, x::AbstractVector, y::AbstractVector, T::InputTrait = input_trait(G))
#     hessian_kernel!(K, G.k, x, y, T)
# end

########################## Isotropic Hessian Kernel ############################
struct IsotropicHessianKernelElement{T, K, R<:AbstractVector{T}, RR<:AbstractMatrix,
                                     UT, CT, KD} <: Factorization{T}
    k::K
    r::R
    rr::RR # r*r'
    U::UT # for [vId, rr]
    C::CT
    kd::KD

    # these could have different element type depending on input
    # Ua::
    # CUa::
    # Ar::AT
    # BA
end
function IsotropicHessianKernelElement(k, x, y)
    # universal variables
    d = length(x)
    Id = I(d)
    vId = vec(Id)

    r = x - y
    rr = r * r'
    vrr = vec(rr)
    T = eltype(r)
    U = hcat(vId, vrr) # micro optimization IDEA: could link second column to rr via view
    C = zeros(T, (2, 2))

    r² = sum(abs2, r)
    kd = derivatives(k, r², 4)[2:end] # get first to fourth derivative
    @. kd = kd * 2^(1:4) # constant adjustment since kernels are function of r², not r²/2

    # Ua = zeros(T, 2)
    # CUa = zeros(T, 2)
    # Ar = zeros(T, d)
    # BA = zeros(T, (d, d))
    IsotropicHessianKernelElement(k, r, rr, U, C, kd) #, Ua, CUa, Ar, BA)
end
function Base.size(K::IsotropicHessianKernelElement)
    d = length(K.r)
    (d^2, d^2)
end
Base.size(K::IsotropicHessianKernelElement, i::Int) = 0 < i ≤ 2 ? size(K)[i] : 1

function allocate_hessian_kernel(k, x, y, ::IsotropicInput)
    IsotropicHessianKernelElement(k, copy(x), copy(y))
end

function evaluate_block!(K::IsotropicHessianKernelElement, G::HessianKernel,
    x::AbstractVector, y::AbstractVector, ::IsotropicInput = IsotropicInput())
    hessian_kernel!(K, G.k, x, y)
end

function hessian_kernel!(K::IsotropicHessianKernelElement, k, x::AbstractVector, y::AbstractVector, ::IsotropicInput = IsotropicInput())
    @. K.r = x-y
    @. K.rr = K.r*K.r'
    r² = sum(abs2, K.r)
    K.kd .= derivatives(k, r², 4)[2:end] # get first to fourth derivative
    @. K.kd = K.kd * 2^(1:4) # constant adjustment since kernels are function of r², not r²/2

    @. K.U[:, 2] = $vec(K.rr) # [vec(Id) vec(rr)], only rr changes
    K.C[1, 1] = K.kd[2]
    K.C[2, 1] = K.kd[3]
    K.C[1, 2] = K.kd[3]
    K.C[2, 2] = K.kd[4]
    return K
end

function Base.Matrix(K::IsotropicHessianKernelElement)
    r = K.r
    rr = vec(K.rr)
    r² = sum(abs2, r)

    d = length(r)
    S = perfect_shuffle(d)

    B = kron(I(d), r*r') + kron(r*r', I(d))
    HH = K.U * K.C * K.U' # Hessian-Hessian
    HH .+= (S + I(d^2)) * (K.kd[2] * I(d^2) + K.kd[3] * B)
end

function LinearAlgebra.mul!(b::AbstractVector, K::IsotropicHessianKernelElement,
                            a::AbstractVector, α::Real = 1, β::Real = 0)
    d = length(K.r)
    d^2 == length(a) == length(b) || throw(DimensionMismatch())
    # below: multiplying by
    # [vec(Id) vec(rr)] * [kd[2] kd[3]; kd[3] kd[4]] * [vec(Id) rr]' i.e. U * C * U'
    # mul!(K.Ua, U', a) # Ua = U'a # this is tricky, because a might have different type
    # mul!(K.CUa, C, K.Ua) # CUa = C*U'a
    Ua = K.U' * a # more general, and these are 2 vectors
    CUa = K.C * Ua
    mul!(b, K.U, CUa, α, β)

    # below: multiplying with (Id2 + Shuffle)
    kd = K.kd
    A = reshape(a, d, d)
    @. b += α * kd[2] * a # I*a
    at = vec(A') # S*a where S is perfect_shuffle
    @. b += α * kd[2] * at

    # below: multiplying with
    # (Id2 + Shuffle) * f3(rn) * B, where
    # B = kronecker(Id, r*r') + kronecker(r*r', Id)
    # BA = zeros(d, d)
    # Ar = zeros(d)
    # mul!(K.Ar, A, K.r)
    # @. K.BA += (K.Ar*K.r') # A * rr' # this could in principle act on b
    # mul!(K.Ar, A', K.r)
    # @. K.BA += (K.r*K.Ar') # rr * A
    # @. b += α * kd[3] * ($vec(K.BA) + $vec(K.BA'))

    Ar = A * K.r # d vector
    BA = Ar * K.r' # A * rr' - d x d matrix
    mul!(Ar, A', K.r)
    @. BA += K.r * Ar' # rr * A
    @. b += α * kd[3] * $vec(BA)
    @. b += α * kd[3] * $vec(BA')
    return b
end

####################### Hessian DotProductInput Kernel #########################
struct DotProductHessianKernelElement{T, K, X<:AbstractVector{T}, Y<:AbstractVector,
                XX<:AbstractMatrix{T}, YY<:AbstractMatrix, KD} <: Factorization{T}
    k::K
    x::X
    y::Y
    xx::XX
    yy::YY
    kd::KD
end
function DotProductHessianKernelElement(k, x, y)
    # universal variables
    d = length(x)
    Id, Id2 = I(d), I(d^2)
    vId = vec(Id)
    xy = dot(x, y)
    xx = x*x'
    yy = y*y'
    kd = derivatives(k, xy, 4)[2:end]
    DotProductHessianKernelElement(k, x, y, xx, yy, kd)
end
function Base.size(K::DotProductHessianKernelElement)
    d = length(K.x)
    (d^2, d^2)
end
Base.size(K::DotProductHessianKernelElement, i::Int) = 0 < i ≤ 2 ? size(K)[i] : 1

function allocate_hessian_kernel(k, x, y, ::DotProductInput)
    DotProductHessianKernelElement(k, copy(x), copy(y))
end
# necessary for BlockFactorization
function evaluate_block!(K::DotProductHessianKernelElement, G::HessianKernel,
    x::AbstractVector, y::AbstractVector, ::DotProductInput = DotProductInput())
    hessian_kernel!(K, G.k, x, y)
end
function hessian_kernel!(K::DotProductHessianKernelElement, k, x::AbstractVector, y::AbstractVector, ::DotProductInput = DotProductInput())
    @. K.x = x
    @. K.y = y
    @. K.xx = x*x'
    @. K.yy = y*y'
    xy = dot(x, y)
    K.kd .= derivatives(k, xy, 4)[2:end] # get first to fourth derivative
    return K
end
function LinearAlgebra.mul!(b::AbstractVector, K::DotProductHessianKernelElement, a::AbstractVector,
                            α::Real, β::Real)
    d = size(K.xx, 1)
    d^2 == length(a) == length(b) || throw(DimensionMismatch())
    xx, yy = vec(K.xx), vec(K.yy)

    # below: multiplying by yy * xx'
    @. b *= β
    xxa = dot(xx, a)
    @. b += α * K.kd[4] * yy * xxa

    # below: multiplying with (Id2 + Shuffle)
    A = reshape(a, d, d)
    @. b += α * K.kd[2] * a # I*a
    @. b += α * K.kd[2] * $vec(A') # + S*a where S is perfect_shuffle

    # below: multiplying with kron(I(d), yx) + kron(yx, I(d))
    BA = zeros(d, d)
    Av = zeros(d)
    mul!(Av, A, K.x)
    @. BA += (Av*K.y') # A * y*x' # this could in principle act on b
    mul!(Av, A', K.x)
    @. BA += (K.y*Av') # y*x' * A
    @. b += α * K.kd[3] * $vec(BA)
    @. b += α * K.kd[3] * $vec(BA')
    return b
end
function Base.Matrix(K::DotProductHessianKernelElement)
    d = length(K.x)
    # universal variables
    Id2 = I(d^2)
    S = perfect_shuffle(d)
    # specific variables
    yx = K.y*K.x'
    xx = vec(K.xx)
    yy = vec(K.yy)
    kd = K.kd # get first to fourth derivative
    B = kron(I(d), yx) + kron(yx, I(d))
    HH = (Id2 + S) * (kd[2] * Id2 + kd[3] * B) + kd[4] * (yy * xx')
end

################################################################################
# kernel for value, gradient, and Hessian observations
struct ValueGradientHessianKernel{T, K} <: MultiKernel{T}
    k::K
    # input_type::IT  IT<:InputTrait
end
input_trait(G::ValueGradientHessianKernel) = input_trait(G.k)
ValueGradientHessianKernel(k::AbstractKernel{T}) where {T} = ValueGradientHessianKernel{T, typeof(k)}(k)
ValueGradientHessianKernel(k) = ValueGradientHessianKernel{Float64, typeof(k)}(k) # use fieldtype here?

# IDEA: input_trait(G.k) could be pre-computed
function (G::ValueGradientHessianKernel)(x::AbstractVector, y::AbstractVector)
    value_gradient_hessian_kernel(G.k, x, y, input_trait(G.k))
end
function value_gradient_hessian_kernel(k, x, y, T::InputTrait)
    K = allocate_value_gradient_hessian_kernel(k, x, y, T)
    value_gradient_hessian_kernel!(K, k, x, y, T)
end
# allocating dense element (this should probably never be done in practice, only comparison)
function allocate_value_gradient_hessian_kernel(k, x, y, ::InputTrait)
    d = length(x)
    n = d^2 + d + 1
    zeros(n, n)
end
function value_gradient_hessian_kernel!(K::AbstractMatrix, k, x::AbstractVector, y::AbstractVector, ::InputTrait)
    d = length(x)
    # block indices
    vi = 1 # value index
    gi = (1+1):(d+1) # gradient index
    hi = (d+1+1):(1+d+d^2) # hessian index
    # value-value
    K[1, 1] = k(x, y)
    # value-gradient
    K[1, gi] = ForwardDiff.gradient(z->k(x, z), y)'
    K[gi, 1] = ForwardDiff.gradient(z->k(z, y), x)
    # gradient-gradient
    K[gi, gi] = ForwardDiff.jacobian(z2->ForwardDiff.gradient(z1->k(z1, z2), x), y)
    # value-hessian
    K[1, hi] = vec(ForwardDiff.hessian(z->k(x, z), y))'
    K[hi, 1] = vec(ForwardDiff.hessian(z->k(z, y), x))
    # gradient-hessian
    K[gi, hi] = ForwardDiff.jacobian(z1->vec(ForwardDiff.hessian(z2->k(z1, z2), y)), x)'
    K[hi, gi] = ForwardDiff.jacobian(z2->vec(ForwardDiff.hessian(z1->k(z1, z2), x)), y)
    # hessian-hessian
    g(y) = vec(ForwardDiff.hessian(z->k(z, y), x)) # covariance between vectorized Hessians
    H = ForwardDiff.jacobian(z->ForwardDiff.jacobian(g, z), y)
    K[hi, hi] .= reshape(H, d^2, d^2)
    return K
end

# this is neccesary when evaluate!(::ValueGradientHessianKernel) is called
# defines mul!
# TODO: increase pre-allocation
# TODO: unify this with DerivativeKernelElement
struct ValueGradientHessianKernelElement{T, K, X<:AbstractVector{T}, Y} <: Factorization{T}
   k::K
   x::X
   y::Y
   # r
   # U::UT # for [vId, rr]
   # Ua
   # CUa
   # C::CT
   # Ar::AT
   # BA
end

function value_gradient_hessian_kernel!(K::ValueGradientHessianKernelElement, k, x::AbstractVector, y::AbstractVector, T::Union{IsotropicInput, DotProductInput})
    K.x .= x; K.y .= y
    return K
end

function evaluate_block!(K::ValueGradientHessianKernelElement, G::ValueGradientHessianKernel,
        x::AbstractVector, y::AbstractVector, T::InputTrait = input_trait(G.k))
    value_gradient_hessian_kernel!(K, G.k, x, y, T)
end

function Base.size(K::ValueGradientHessianKernelElement)
    d = length(K.x)
    n = d^2 + d + 1
    (n, n)
end
Base.size(K::ValueGradientHessianKernelElement, i::Int) = 0 < i ≤ 2 ? size(K)[i] : 1

function Base.Matrix(K::ValueGradientHessianKernelElement) #, ::DotProductInput = input_trait(K.k))
    k, x, y = K.k, K.x, K.y
    d = length(x)
    n = d^2 + d + 1
    A = zeros(n, n)
    value_gradient_hessian_kernel!(A, k, x, y, input_trait(k))
    # d = length(x)
    # # universal variables
    # Id2 = I(d^2)
    # S = perfect_shuffle(d)
    #
    # # specific variables
    # xy = dot(x, y)
    # yx = y*x'
    # xx = vec(x*x')
    # yy = vec(y*y')
    #
    # kd = derivatives(k, xy, 4)[2:end] # get first to fourth derivative
    #
    # B = kron(I(d), yx) + kron(yx, I(d))
    # HH = (Id2 + S) * (kd[2] * Id2 + kd[3] * B) +  kd[4] * (yy * xx')
end
function allocate_value_gradient_hessian_kernel(k, x, y, ::IsotropicInput)
    ValueGradientHessianKernelElement(k, copy(x), copy(y))
end

function LinearAlgebra.mul!(b::AbstractVector, K::ValueGradientHessianKernelElement,
                            a::AbstractVector, α::Real = 1, β::Real = 0)
    mul!(b, K, a, α, β, input_trait(K.k))
end
# fast multiplication with isotropic second-order derivative kernel
function LinearAlgebra.mul!(b::AbstractVector, K::ValueGradientHessianKernelElement,
                            a::AbstractVector, α::Real, β::Real, T::IsotropicInput)
    k, x, y = K.k, K.x, K.y
    d = length(x)
    # block indices
    vi = 1 # value index
    gi = 2:(d+1) # gradient index
    hi = (d+2):(1+d+d^2) # hessian index

    avi, agi, ahi = @views a[vi], a[gi], a[hi]

    # universal variables
    Id, Id2 = I(d), I(d^2)
    vId = vec(Id)

    # specific variables
    r = x - y
    rr = vec(r*r')
    r² = norm(r)^2

    # evaluate kernel and derivatives
    k = K.k
    kd = derivatives(k, r², 4)
    @. kd *= 2^(0:4)
    k0, k1, k2, k3, k4 = kd

    @. b *= β # scaling of target

    # important inner products
    ra = dot(r, agi)
    ida = dot(vId, ahi)
    rra = dot(rr, ahi)

    # value
    b[vi] += α * k0 * a[vi] # value
    b[vi] -= α * k1 * ra # gradient
    b[vi] += α * (k1 * ida + k2 * rra) # hessian

    # gradient
    @. b[gi] += α * k1 * r * a[vi]
    @. b[gi] -= α * (k1 * agi + k2 * r * ra)
    @. b[gi] += α * (k2 * ida + k3 * rra) * r # this sign seems correct

    BA = zeros(d)
    Ahi = reshape(ahi, d, d)
    BA += Ahi * r + vec(r' * Ahi) # can be combined with the below
    @. b[gi] += α * k2 * BA

    # hessian
    @. b[hi] += α * (k1 * vId + k2 * rr) * a[vi] # value part of gradient
    @. b[hi] -= α * (k2 * vId + k3 * rr) * ra # this can be collected with the above
    BA = zeros(d, d)
    @. BA += agi * r' + r * agi' # can be combined with the below
    @. b[hi] -= α * k2 * $vec(BA)

    # hessian-hessian
    # below: multiplying by
    # [vec(Id) rr] * [kd[2] kd[3]; kd[3] kd[4]] * [vec(Id) rr]'
    U = hcat(vId, rr)
    C = zeros(2, 2)
    C[1, 1] = k2
    C[2, 1] = k3
    C[1, 2] = k3
    C[2, 2] = k4
    # Ua = U'a
    Ua = [ida, rra] # (U'*ahi)
    bhi = @view b[hi]
    mul!(bhi, U, C*Ua, α, 1)

    # below: multiplying with (Id2 + Shuffle) * ahi
    Ahi = reshape(ahi, d, d)
    @. b[hi] += α * k2 * ahi # I*a
    @. b[hi] += α * k2 * $vec(Ahi') # + S*a where S is perfect_shuffle

    # below: multiplying with
    # (Id2 + Shuffle) * f3(rn) * (Bi - Bj), where
    # B = kronecker(Id, r*r') + kronecker(r*r', Id)
    BA = zeros(d, d)
    Ar = zeros(d)
    mul!(Ar, Ahi, r)
    @. BA += (Ar*r') # A * rr' # this could in principle act on b
    mul!(Ar, Ahi', r)
    @. BA += (r*Ar') # rr * A

    @. b[hi] += α * k3 * $vec(BA)
    @. b[hi] += α * k3 * $vec(BA') # Adding Shuffle * A
    return b
end

# # multiplying with rank 1 kronecker sum
# # function _kron_sum_mul!(b::AbstractVector, x::AbstractVector, y::Abstractor,
# #                         a::AbstractVector, α::Real = 1, β::Real = 0)
# #
# # end
#

# block implementation probably less efficient than the Element one
# TODO: move to value_gradient_hessian_kernel!
# function value_gradient_hessian_kernel(K::BlockFactorization, k, x::AbstractVector, y::AbstractVector, T::DotProductInput)
#
#     # first, call ValueGradientKernel to calculate value-gradient covariances in upper 2x2 block of K
#     value_gradient_kernel!(K, k, x, y, T)
#
#     # universal variables
#     Id, Id2 = I(d), I(d^2)
#     vId = vec(Id)
#     S = perfect_shuffle(d)
#
#     # # specific variables
#     # r = x - y
#     # rr = vec(r*r')
#     # r² = norm(r)^2
#
#     kd = derivatives(k, r², 4)[2:end] # get first to fourth derivative
#     kd = @. kd * 2^(1:4) # constant adjustment
#
#     # value-hessian
#     K[1, 3] = kd[1]
#     K[3, 1] = -1
#
#     # gradient-hessian
#     K[2, 3] = -1
#     K[3, 2] = -1
#
#     # hessian-hessian
#     K[3, 3] = -1
# end


text:
function StatsBase.fit(::Type{Normal}, m::AbstractMoments)
    n_moments(m) >= 2 || error("Need mean and variance to estimate normal")
    return Normal(mean(m),std(m))
end

function StatsBase.fit(::Type{Normal}, lower::QuantilePoint, upper::QuantilePoint)
    # https://www.johndcook.com/quantiles_parameters.pdf
    if (upper < lower)
        lower,upper = (upper,lower)
    end
    qz1 = quantile(Normal(), lower.p)
    qz2 = quantile(Normal(), upper.p)
    dqz = (qz2 - qz1)
    σ = (upper.q - lower.q)/dqz
    μ = (lower.q*qz2 - upper.q*qz1)/dqz
    Normal(μ,σ)
end

fit_mean_quantile(D::Type{Normal}, mean::Real, qp::QuantilePoint) =
    fit(D, QuantilePoint(mean,0.5), qp)

fit_mode_quantile(D::Type{Normal}, mode::Real, qp::QuantilePoint) =
    fit(D, QuantilePoint(mode,0.5), qp)