text:
using Rotations
using Oceananigans.Grids
using Oceananigans.Grids: R_Earth, interior_indices

import Base: show, size, eltype
import Oceananigans.Grids: topology, architecture, halo_size, on_architecture

struct CubedSphereFaceConnectivityDetails{F, S}
    face :: F
    side :: S
end

short_string(deets::CubedSphereFaceConnectivityDetails) = "face $(deets.face) $(deets.side) side"

Base.show(io::IO, deets::CubedSphereFaceConnectivityDetails) =
    print(io, "CubedSphereFaceConnectivityDetails: $(short_string(deets))")

struct CubedSphereFaceConnectivity{W, E, S, N}
     west :: W
     east :: E
    south :: S
    north :: N
end

CubedSphereFaceConnectivity(; west, east, south, north) =
    CubedSphereFaceConnectivity(west, east, south, north)

function Base.show(io::IO, connectivity::CubedSphereFaceConnectivity)
    print(io, "CubedSphereFaceConnectivity:\n",
              "├── west: $(short_string(connectivity.west))\n",
              "├── east: $(short_string(connectivity.east))\n",
              "├── south: $(short_string(connectivity.south))\n",
              "└── north: $(short_string(connectivity.north))")
end

function default_face_connectivity()
    # See figure 8.4 of https://mitgcm.readthedocs.io/en/latest/phys_pkgs/exch2.html?highlight=cube%20sphere#fig-6tile
    #
    #                         face  F5   face  F6
    #                       +----------+----------+
    #                       |    ↑↑    |    ↑↑    |
    #                       |    1W    |    1S    |
    #                       |←3N F5 6W→|←5E F6 2S→|
    #                       |    4N    |    4E    |
    #              face  F3 |    ↓↓    |    ↓↓    |
    #            +----------+----------+----------+
    #            |    ↑↑    |    ↑↑    |
    #            |    5W    |    5S    |
    #            |←1N F3 4W→|←3E F4 6S→|
    #            |    2N    |    2E    |
    #            |    ↓↓    |    ↓↓    |
    # +----------+----------+----------+
    # |    ↑↑    |    ↑↑    | face  F4
    # |    3W    |    3S    |
    # |←5N F1 2W→|←1E F2 4S→|
    # |    6N    |    6E    |
    # |    ↓↓    |    ↓↓    |
    # +----------+----------+
    #   face  F1   face  F2

    face1_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(5, :north),
        east  = CubedSphereFaceConnectivityDetails(2, :west),
        south = CubedSphereFaceConnectivityDetails(6, :north),
        north = CubedSphereFaceConnectivityDetails(3, :west),
    )

    face2_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(1, :east),
        east  = CubedSphereFaceConnectivityDetails(4, :south),
        south = CubedSphereFaceConnectivityDetails(6, :east),
        north = CubedSphereFaceConnectivityDetails(3, :south),
    )

    face3_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(1, :north),
        east  = CubedSphereFaceConnectivityDetails(4, :west),
        south = CubedSphereFaceConnectivityDetails(2, :north),
        north = CubedSphereFaceConnectivityDetails(5, :west),
    )

    face4_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(3, :east),
        east  = CubedSphereFaceConnectivityDetails(6, :south),
        south = CubedSphereFaceConnectivityDetails(2, :east),
        north = CubedSphereFaceConnectivityDetails(5, :south),
    )

    face5_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(3, :north),
        east  = CubedSphereFaceConnectivityDetails(6, :west),
        south = CubedSphereFaceConnectivityDetails(4, :north),
        north = CubedSphereFaceConnectivityDetails(1, :west),
    )

    face6_connectivity = CubedSphereFaceConnectivity(
        west  = CubedSphereFaceConnectivityDetails(5, :east),
        east  = CubedSphereFaceConnectivityDetails(2, :south),
        south = CubedSphereFaceConnectivityDetails(4, :east),
        north = CubedSphereFaceConnectivityDetails(1, :south),
    )

    face_connectivity = (
        face1_connectivity,
        face2_connectivity,
        face3_connectivity,
        face4_connectivity,
        face5_connectivity,
        face6_connectivity
    )

    return face_connectivity
end

# Note: I think we want to keep faces and face_connectivity tuples
# so it's easy to support an arbitrary number of faces.

struct ConformalCubedSphereGrid{FT, F, C, Arch} <: AbstractHorizontallyCurvilinearGrid{FT, FullyConnected, FullyConnected, Bounded, Arch}
         architecture :: Arch
                faces :: F
    face_connectivity :: C
end

function ConformalCubedSphereGrid(arch = CPU(), FT=Float64; face_size, z, radius=R_Earth)
    @warn "ConformalCubedSphereGrid is experimental: use with caution!"

    # +z face (face 1)
    z⁺_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=nothing)

    # +x face (face 2)
    x⁺_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=RotX(π/2))

    # +y face (face 3)
    y⁺_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=RotY(π/2))

    # -x face (face 4)
    x⁻_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=RotX(-π/2))

    # -y face (face 5)
    y⁻_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=RotY(-π/2))

    # -z face (face 6)
    z⁻_face_grid = ConformalCubedSphereFaceGrid(arch, FT, size=face_size, z=z, radius=radius, rotation=RotX(π))

    faces = (
        z⁺_face_grid,
        x⁺_face_grid,
        y⁺_face_grid,
        x⁻_face_grid,
        y⁻_face_grid,
        z⁻_face_grid
    )

    face_connectivity = default_face_connectivity()

    return ConformalCubedSphereGrid{FT, typeof(faces), typeof(face_connectivity), typeof(arch)}(arch, faces, face_connectivity)
end

function ConformalCubedSphereGrid(filepath::AbstractString, arch = CPU(), FT=Float64; Nz, z, radius = R_Earth, halo = (1, 1, 1))
    @warn "ConformalCubedSphereGrid is experimental: use with caution!"

    face_topo = (FullyConnected, FullyConnected, Bounded)
    face_kwargs = (Nz=Nz, z=z, topology=face_topo, radius=radius, halo=halo)

    faces = Tuple(ConformalCubedSphereFaceGrid(filepath, arch, FT; face=n, face_kwargs...) for n in 1:6)

    face_connectivity = default_face_connectivity()

    grid = ConformalCubedSphereGrid{FT, typeof(faces), typeof(face_connectivity), typeof(arch)}(arch, faces, face_connectivity)

    fill_grid_metric_halos!(grid)
    fill_grid_metric_halos!(grid)

    return grid
end

function Base.show(io::IO, grid::ConformalCubedSphereGrid{FT}) where FT
    Nx, Ny, Nz, Nf = size(grid)
    print(io, "ConformalCubedSphereGrid{$FT}: $Nf faces with size = ($Nx, $Ny, $Nz)")
end

#####
##### Nodes for ConformalCubedSphereFaceGrid
#####

@inline λnode(LX::Face,   LY::Face,   LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.λᶠᶠᵃ[i, j]
@inline λnode(LX::Face,   LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.λᶠᶜᵃ[i, j]
@inline λnode(LX::Center, LY::Face,   LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.λᶜᶠᵃ[i, j]
@inline λnode(LX::Center, LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.λᶜᶜᵃ[i, j]

@inline φnode(LX::Face,   LY::Face,   LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.φᶠᶠᵃ[i, j]
@inline φnode(LX::Face,   LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.φᶠᶜᵃ[i, j]
@inline φnode(LX::Center, LY::Face,   LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.φᶜᶠᵃ[i, j]
@inline φnode(LX::Center, LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.φᶜᶜᵃ[i, j]

@inline znode(LX, LY, LZ::Face,   i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.zᵃᵃᶠ[k]
@inline znode(LX, LY, LZ::Center, i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.zᵃᵃᶜ[k]

λnodes(LX::Face, LY::Face, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.λᶠᶠᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

λnodes(LX::Face, LY::Center, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.λᶠᶜᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

λnodes(LX::Center, LY::Face, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.λᶜᶠᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

λnodes(LX::Center, LY::Center, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.λᶜᶜᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

φnodes(LX::Face, LY::Face, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.φᶠᶠᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

φnodes(LX::Face, LY::Center, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.φᶠᶜᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

φnodes(LX::Center, LY::Face, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.φᶜᶠᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

φnodes(LX::Center, LY::Center, LZ, grid::ConformalCubedSphereFaceGrid{TX, TY}) where {TX, TY} =
    view(grid.φᶜᶜᵃ, interior_indices(LX, TX, grid.Nx), interior_indices(LY, TY, grid.Ny))

#####
##### Grid utils
#####

Base.size(grid::ConformalCubedSphereGrid) = (size(grid.faces[1])..., length(grid.faces))
Base.size(grid::ConformalCubedSphereGrid, i) = size(grid)[i]
halo_size(ccsg::ConformalCubedSphereGrid) = halo_size(first(ccsg.faces)) # hack

Base.eltype(grid::ConformalCubedSphereGrid{FT}) where FT = FT

topology(::ConformalCubedSphereGrid) = (Bounded, Bounded, Bounded)
architecture(grid::ConformalCubedSphereGrid) = grid.architecture

function on_architecture(arch, grid::ConformalCubedSphereGrid)

    faces = Tuple(on_architecture(arch, grid.faces[n]) for n in 1:6)
    face_connectivity = grid.face_connectivity
    FT = eltype(grid)

    return ConformalCubedSphereGrid{FT, typeof(faces), typeof(face_connectivity), typeof(arch)}(arch, faces, face_connectivity)
end

#####
##### filling grid halos
#####

function grid_metric_halo(grid_metric, grid, location, topo, side)
    LX, LY = location
    TX, TY = topo
    side == :west  && return  underlying_west_halo(grid_metric, grid, LX, TX)
    side == :east  && return  underlying_east_halo(grid_metric, grid, LX, TX)
    side == :south && return underlying_south_halo(grid_metric, grid, LY, TY)
    side == :north && return underlying_north_halo(grid_metric, grid, LY, TY)
end

function grid_metric_boundary(grid_metric, grid, location, topo, side)
    LX, LY = location
    TX, TY = topo
    side == :west  && return  underlying_west_boundary(grid_metric, grid, LX, TX)
    side == :east  && return  underlying_east_boundary(grid_metric, grid, LX, TX)
    side == :south && return underlying_south_boundary(grid_metric, grid, LY, TY)
    side == :north && return underlying_north_boundary(grid_metric, grid, LY, TY)
end

function fill_grid_metric_halos!(grid)

    topo_bb = (Bounded, Bounded)

    loc_cc = (Center, Center)
    loc_cf = (Center, Face  )
    loc_fc = (Face,   Center)
    loc_ff = (Face,   Face  )

    for face_number in 1:6, side in (:west, :east, :south, :north)

        connectivity_info = getproperty(grid.face_connectivity[face_number], side)
        src_face_number = connectivity_info.face
        src_side = connectivity_info.side

        grid_face = grid.faces[face_number]
        src_grid_face = grid.faces[src_face_number]

        if sides_in_the_same_dimension(side, src_side)
            grid_metric_halo(grid_face.Δxᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= grid_metric_boundary(grid_face.Δxᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side)
            grid_metric_halo(grid_face.Δyᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= grid_metric_boundary(grid_face.Δyᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side)
            grid_metric_halo(grid_face.Azᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= grid_metric_boundary(grid_face.Azᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side)

            grid_metric_halo(grid_face.Δxᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= grid_metric_boundary(grid_face.Δxᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side)
            grid_metric_halo(grid_face.Δyᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= grid_metric_boundary(grid_face.Δyᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side)
            grid_metric_halo(grid_face.Azᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= grid_metric_boundary(grid_face.Azᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side)

            grid_metric_halo(grid_face.Δxᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= grid_metric_boundary(grid_face.Δxᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side)
            grid_metric_halo(grid_face.Δyᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= grid_metric_boundary(grid_face.Δyᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side)
            grid_metric_halo(grid_face.Azᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= grid_metric_boundary(grid_face.Azᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side)

            grid_metric_halo(grid_face.Δxᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= grid_metric_boundary(grid_face.Δxᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side)
            grid_metric_halo(grid_face.Δyᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= grid_metric_boundary(grid_face.Δyᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side)
            grid_metric_halo(grid_face.Azᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= grid_metric_boundary(grid_face.Azᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side)
        else
            reverse_dim = src_side in (:west, :east) ? 1 : 2
            grid_metric_halo(grid_face.Δxᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δyᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Δyᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δxᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Azᶜᶜᵃ, grid_face, loc_cc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Azᶜᶜᵃ, src_grid_face, loc_cc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)

            grid_metric_halo(grid_face.Δxᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δyᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Δyᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δxᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Azᶜᶠᵃ, grid_face, loc_cf, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Azᶠᶜᵃ, src_grid_face, loc_fc, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)

            grid_metric_halo(grid_face.Δxᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δyᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Δyᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δxᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Azᶠᶜᵃ, grid_face, loc_fc, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Azᶜᶠᵃ, src_grid_face, loc_cf, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)

            grid_metric_halo(grid_face.Δxᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δyᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Δyᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Δxᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
            grid_metric_halo(grid_face.Azᶠᶠᵃ, grid_face, loc_ff, topo_bb, side) .= reverse(permutedims(grid_metric_boundary(grid_face.Azᶠᶠᵃ, src_grid_face, loc_ff, topo_bb, src_side), (2, 1, 3)), dims=reverse_dim)
        end
    end

    return nothing
end


text:
defaultfigure(;kwargs...) = Figure(
    ;resolution = (800, 800),
    background = RGBA(0, 0, 0, 0),
    kwargs...)

# ## Plotting interface definition

"""
    plotsample(method, sample)
"""
function plotsample(method, sample)
    f = defaultfigure(resolution = (300, 150))
    plotsample!(f, method, sample)
    return f
end

"""
    plotsample!(f, method, sample)
"""
function plotsample! end

"""
    plotxy(method, (x, y))
"""
function plotxy(method, xy)
    f = defaultfigure(resolution = (300, 150))
    plotxy!(f, method, xy)
    return f
end
"""
    plotxy!(f, method, (x, y))
"""
function plotxy! end

function plotbatch(method, (xs, ys))
    n = size(xs)[end]
    nrows = Int(ceil(sqrt(n)))
    f = Figure()
    is = Iterators.product(1:nrows, 1:nrows)
    for (i, (x, y)) in zip(is, DataLoaders.obsslices((xs, ys)))
        plotxy!(f[i...], method, (x, y))
    end
    return f
end

# ## Utilities

function imageaxis(f; kwargs...)
    ax = AbstractPlotting.Axis(f; kwargs...)
    ax.aspect = DataAspect()
    ax.xzoomlock = true
    ax.yzoomlock = true
    ax.xrectzoom = false
    ax.yrectzoom = false
    ax.panbutton = nothing
    ax.xpanlock = true
    ax.ypanlock = true
    ax.bottomspinevisible = false
    ax.leftspinevisible = false
    ax.rightspinevisible = false
    ax.topspinevisible = false
    MakieLayout.tightlimits!(ax)
    hidedecorations!(ax)

    return ax
end

imageaxis(f::AbstractPlotting.FigurePosition; kwargs...) = imageaxis(f.fig; kwargs...)

# ## Plot recipes

@recipe(PlotImage, image) do scene
    Attributes()
end

function AbstractPlotting.plot!(plot::PlotImage)
    im = plot[:image]
    rim = @lift rotr90($im)
    image!(plot, rim; plot.attributes...)
    return plot
end

@recipe(PlotMask, mask, classes) do scene
    Attributes()
end

function AbstractPlotting.plot!(plot::PlotMask; kwargs...)
    mask = plot[:mask]
    classes = try
        classes = plot[:classes]
    catch
        classes = @lift unique($mask)
    end
    im = @lift maskimage($mask, $classes)
    plotimage!(plot, im; plot.attributes...)
    return plot
end

function maskimage(mask, classes)
    colors = distinguishable_colors(length(classes), transform = deuteranopic)
    im = map(c -> colors[c], mask)
end

maskimage(mask::AbstractArray{<:Gray{T}}, args...) where T =
    maskimage(reinterpret(T, mask), args...)
maskimage(mask::AbstractArray{<:Normed{T}}, args...) where T =
    maskimage(reinterpret(T, mask), args...)

ext: jl


text:
_init_field(t::Type) = _init_field(Val{t}())
_init_field(::Val{T}) where {T <: AbstractMatrix} = zeros(eltype(T), 0, 0)
_init_field(::Val{T}) where {T <: SparseMatrixCSC} = sparse(zeros(eltype(T), 0, 0))
_init_field(::Val{T}) where {T <: AbstractVector} = zeros(eltype(T), 0)
_init_field(::Val{T}) where {T <: SparseVector} = sparse(zeros(eltype(T), 0))
_init_field(::Val{BigFloat}) = BigFloat(NaN)
_init_field(::Val{Float64}) = NaN
_init_field(::Val{Float32}) = NaN32
_init_field(::Val{Float16}) = NaN16
_init_field(::Val{Missing}) = missing
_init_field(::Val{Nothing}) = nothing
_init_field(::Val{Bool}) = false #unfortunately no way to differentiate
_init_field(::Val{T}) where {T <: Number} = typemin(T)
_init_field(::Val{Counters}) = Counters()

"""
Type: GenericState

Methods: update!, reinit!

A generic State to describe the state of a problem at a point x.

Tracked data include:
 - x             : current iterate
 - d [opt]       : search direction
 - res [opt]     : residual
 - current_time  : time
 - current_score : score

Constructors:
 `GenericState(:: T, :: S; d :: T = _init_field(T), res :: T = _init_field(T), current_time :: Float64 = NaN) where {S, T <:AbstractVector}`

 `GenericState(:: T; d :: T = _init_field(T), res :: T = _init_field(T), current_time :: Float64 = NaN, current_score :: Union{T,eltype(T)} = _init_field(eltype(T))) where T <:AbstractVector`

Note:
 - By default, unknown entries are set using `_init_field`.
 - By default the type of `current_score` is `eltype(x)` and cannot be changed once the State is created.
   To have a vectorized `current_score` of length n, try something like `GenericState(x, Array{eltype(x),1}(undef, n))`.

Examples:
  GenericState(x)
  GenericState(x, Array{eltype(x),1}(undef, length(x)))
  GenericState(x, current\\_time = 1.0)
  GenericState(x, current\\_score = 1.0)

See also: Stopping, NLPAtX
"""
mutable struct GenericState{S, T <: Union{AbstractFloat, AbstractVector}} <: AbstractState{S, T}
  x::T

  d::T
  res::T

  #Current time
  current_time::Float64
  #Current score
  current_score::S

  function GenericState(
    x::T,
    current_score::S;
    d::T = _init_field(T),
    res::T = _init_field(T),
    current_time::Float64 = NaN,
  ) where {S, T <: AbstractVector}
    return new{S, T}(x, d, res, current_time, current_score)
  end
end

function GenericState(
  x::T;
  d::T = _init_field(T),
  res::T = _init_field(T),
  current_time::Float64 = NaN,
  current_score::Union{T, eltype(T)} = _init_field(eltype(T)),
) where {T <: AbstractVector}
  return GenericState(x, current_score, d = d, res = res, current_time = current_time)
end

scoretype(typestate::AbstractState{S, T}) where {S, T} = S
xtype(typestate::AbstractState{S, T}) where {S, T} = T

"""
update!: generic update function for the State

`update!(:: AbstractState; convert = false, kwargs...)`

The function compares the kwargs and the entries of the State.
If the type of the kwargs is the same as the entry, then
it is updated.

Set kargs `convert` to true to update even incompatible types.

Examples:
update!(state1)
update!(state1, current\\_time = 2.0)
update!(state1, convert = true, current\\_time = 2.0)

See also: GenericState, reinit!, update\\_and\\_start!, update\\_and\\_stop!
"""
function update!(stateatx::T; convert::Bool = false, kwargs...) where {T <: AbstractState}
  fnames = fieldnames(T)
  for k ∈ keys(kwargs)
    #check if k is in fieldnames and type compatibility
    if (k ∈ fnames) && typeof(kwargs[k]) <: typeof(getfield(stateatx, k))
      setfield!(stateatx, k, kwargs[k])
    end
  end

  return stateatx
end

#Ca serait cool d'avoir un shortcut en repérant certains keywords
#ou si il n'y a aucun keyword!
#function update!(stateatx :: T; x :: TT = stateatx.x) where {TS, TT, T <: AbstractState{TS,TT}}
# setfield!(stateatx, :x, x)
# return stateatx
#end

"""
\\_smart\\_update!: generic update function for the State without Type verification.

`_smart_update!(:: AbstractState; kwargs...)`

The function works exactly as update! without type and field verifications.
So, affecting a value to nothing or a different type will return an error.

See also: update!, GenericState, reinit!, update\\_and\\_start!, update\\_and\\_stop!
"""
function _smart_update!(stateatx::T; kwargs...) where {T <: AbstractState}
  for k ∈ keys(kwargs)
    setfield!(stateatx, k, kwargs[k])
  end

  return stateatx
end
#https://github.com/JuliaLang/julia/blob/f3252bf50599ba16640ef08eb1e43c632eacf264/base/Base.jl#L34

function _update_time!(stateatx::T, current_time::Float64) where {T <: AbstractState}
  setfield!(stateatx, :current_time, current_time)

  return stateatx
end

"""
reinit!: function that set all the entries at `_init_field` except the mandatory `x`.

`reinit!(:: AbstractState, :: T; kwargs...)`

Note: If `x` is given as a kargs it will be prioritized over
the second argument.

Examples:
reinit!(state2, zeros(2))
reinit!(state2, zeros(2), current_time = 1.0)

There is a shorter version of reinit! reusing the `x` in the state

`reinit!(:: AbstractState; kwargs...)`

Examples:
reinit!(state2)
reinit!(state2, current_time = 1.0)
"""
function reinit!(stateatx::St, x::T; kwargs...) where {S, T, St <: AbstractState{S, T}}

  #for k not in the kwargs
  for k ∈ setdiff(fieldnames(St), keys(kwargs))
    if k != :x
      setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
    end
  end

  setfield!(stateatx, :x, x)

  if length(kwargs) == 0
    return stateatx #save the update! call if no other kwargs than x
  end

  return update!(stateatx; kwargs...)
end

function reinit!(stateatx::T; kwargs...) where {T <: AbstractState}
  for k ∈ setdiff(fieldnames(T), keys(kwargs))
    if k != :x
      setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
    end
  end

  return update!(stateatx; kwargs...)
end

"""
\\_domain\\_check: returns true if there is a `NaN` or a `Missing` in the state entries (short-circuiting), false otherwise.

`_domain_check(:: AbstractState; kwargs...)`

Note:
- The fields given as keys in kwargs are not checked.

Examples:
\\_domain\\_check(state1)
\\_domain\\_check(state1, x = true)
"""
function _domain_check(stateatx::T; kwargs...) where {T <: AbstractState}
  for k ∈ setdiff(fieldnames(T), keys(kwargs))
    gf = getfield(stateatx, k)
    if _check_nan_miss(gf)
      return true
    end
  end

  return false
end

_check_nan_miss(field::Any) = false #Nothing or Counters
_check_nan_miss(field::SparseMatrixCSC) = any(isnan, field.nzval) #because checking in sparse matrices is too slow
_check_nan_miss(field::Union{AbstractVector, AbstractMatrix}) = any(isnan, field)
#We don't check for NaN's in Float as they are the _init_field
_check_nan_miss(field::AbstractFloat) = ismissing(field)

import Base.copy
ex = :(_genobj(typ) = $(Expr(:new, :typ)));
eval(ex);
function copy(state::T) where {T <: AbstractState}

  #ex=:(_genobj(typ)=$(Expr(:new, :typ))); eval(ex)
  cstate = _genobj(T)
  #cstate = $(Expr(:new, typeof(state)))

  for k ∈ fieldnames(T)
    setfield!(cstate, k, deepcopy(getfield(state, k)))
  end

  return cstate
end

"""
compress_state!: compress State with the following rules.
- If it contains matrices and save_matrix is false, then the corresponding entries
are set to _init_field(typeof(getfield(stateatx, k)).
- If it contains vectors with length greater than max_vector_size, then the
corresponding entries are replaced by a vector of size 1 containing its pnorm-norm.
- If keep is true, then only the entries given in kwargs will be saved (the others are set to _init_field(typeof(getfield(stateatx, k))).
- If keep is false and an entry in the State is in the kwargs list, then it is put as _init_field(typeof(getfield(stateatx, k)) if possible.

`compress_state!(:: AbstractState; save_matrix :: Bool = false, max_vector_size :: Int = length(stateatx.x), pnorm :: Real = Inf, keep :: Bool = false, kwargs...)`

see also: copy, copy\\_compress\\_state, ListofStates
"""
function compress_state!(
  stateatx::T;
  save_matrix::Bool = false,
  max_vector_size::Int = length(stateatx.x),
  pnorm::Real = Inf,
  keep::Bool = false,
  kwargs...,
) where {T <: AbstractState}
  for k ∈ fieldnames(T)
    if k ∈ keys(kwargs) && !keep
      try
        setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
      catch
        #nothing happens
      end
    end
    if k ∉ keys(kwargs) && keep
      try
        setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
      catch
        #nothing happens
      end
    end
    if typeof(getfield(stateatx, k)) <: AbstractVector
      katt = getfield(stateatx, k)
      if (length(katt) > max_vector_size)
        setfield!(stateatx, k, [norm(katt, pnorm)])
      end
    elseif typeof(getfield(stateatx, k)) <: Union{AbstractArray, AbstractMatrix}
      if save_matrix
        katt = getfield(stateatx, k)
        if maximum(size(katt)) > max_vector_size
          setfield!(stateatx, k, norm(getfield(stateatx, k)) * ones(1, 1))
        end
      else #save_matrix is false
        setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
      end
    else
      #nothing happens
    end
  end

  return stateatx
end

"""
copy\\_compress\\_state: copy the State and then compress it.

`copy_compress_state(:: AbstractState; save_matrix :: Bool = false, max_vector_size :: Int = length(stateatx.x), pnorm :: Real = Inf, kwargs...)`

see also: copy, compress_state!, ListofStates
"""
function copy_compress_state(
  stateatx::AbstractState;
  save_matrix::Bool = false,
  max_vector_size::Int = length(stateatx.x),
  pnorm::Real = Inf,
  kwargs...,
)
  cstate = copy(stateatx)
  return compress_state!(
    cstate;
    save_matrix = save_matrix,
    max_vector_size = max_vector_size,
    pnorm = pnorm,
    kwargs...,
  )
end

import Base.show
function show(io::IO, state::AbstractState)
  varlines = "$(typeof(state)) with an iterate of type $(xtype(state)) and a score of type $(scoretype(state))."
  println(io, varlines)
end


text:
function TransCart(ksi,eta,z,F,Topo,Param)
X=zeros(3);
X[1:2]=0.25*((1-ksi)*(1-eta)*F.P[1:2,1]+
  (1+ksi)*(1-eta)*F.P[1:2,2]+
  (1+ksi)*(1+eta)*F.P[1:2,3]+
  (1-ksi)*(1+eta)*F.P[1:2,4]);
X[3]=first(Topo(X[1],X[2],z,Param));
return X
end


text:
using Nettle
using Base.Test

# MD5 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,digest) in {
    (
        "",
        "",
        "74e6f7298a9c2d168935f58c001bad88"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "80070713463e7749b90c2dc24911e275"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b"),
        "Hi There",
        "9294727a3638bb1c13f48ef8158bfc9d"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "750c783e6ab0b503eaa86e310a5db738"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "56be34521d144c88dbb8c733f0e8b3f6"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "697eaf0aca3a3aea3a75164746ffaa79"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "6b1ab7fe4bd7bf8f0b62e6ce61b9d0cd"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key and Larger " *
        "Than One Block-Size Data",
        "6f630fad67cda0ee1fb1f562db3aa53e"
    )
}
    @test md5_hmac(key,text) == hex2bytes(digest)
end

# SHA1 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,digest) in {
    (
        "",
        "",
        "fbdb1d1b18aa6c08324b7d64b71fb76370690e1d"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "de7c9b85b8b78aa6bc8a7a36f70a90701c9db4d9"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b"),
        "Hi There",
        "b617318655057264e28bc0b6fb378c8ef146be00"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "effcdf6ae5eb2fa2d27416d5f184df9c259a7c79"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "125d7342b9ac11cd91a39af48aa17b4f63f175d3"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "4c9007f4026250c6bc8414f9bf50c86c2d7235da"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "aa4ae5e15272d00e95705637ce8a3b55ed402112"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key and Larger " *
        "Than One Block-Size Data",
        "e8e99d0f45237d786d6bbaa7965c7808bbff1a91"
    )
}
    @test sha1_hmac(key,text) == hex2bytes(digest)
end

# SHA256 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,digest) in {
   (
        "",
        "",
        "b613679a0814d9ec772f95d778c35fc5ff1697c493715653c6c712144292c5ad"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "f7bc83f430538424b13298e6aa6fb143ef4d59a14946175997479dbc2d1a3cd8"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b" *
        "0b0b0b0b"),
        "Hi There",
        "b0344c61d8db38535ca8afceaf0bf12b" *
        "881dc200c9833da726e9376c2e32cff7"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "5bdcc146bf60754e6a042426089575c7" *
        "5a003f089d2739839dec58b964ec3843"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "773ea91e36800e46854db8ebd09181a7" *
        "2959098b3ef8c122d9635514ced565fe"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "82558a389a443c0ea4cc819899f2083a" *
        "85f0faa3e578f8077a2e3ff46729665b"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "60e431591ee0b67f0d8a26aacbf5b77f" *
        "8e0bc6213728c5140546040f0ee37f54"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaa"),
        "This is a test using a larger than block-size ke" *
        "y and a larger than block-size data. The key nee" *
        "ds to be hashed before being used by the HMAC al" *
        "gorithm.",
        "9b09ffa71b942fcb27635fbcd5b0e944" *
        "bfdc63644f0713938a7f51535c3a35e2"
    )
}
    @test sha256_hmac(key,text) == hex2bytes(digest)
end

# MD5 hash tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/md5-test.c
for (text, hash) in {
    (
        "",
        "d41d8cd98f00b204e9800998ecf8427e"
    ),(
        "a",
        "0cc175b9c0f1b6a831c399e269772661"
    ),(
        "abc",
        "900150983cd24fb0d6963f7d28e17f72"
    ),(
        "message digest",
        "f96b697d7cb7938d525a2f31aaf161d0"
    ),(
        "abcdefghijklmnopqrstuvwxyz",
        "c3fcd3d76192e4007dfb496cca67e13b"
    ),(
        "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdef" *
        "ghijklmnopqrstuvwxyz0123456789",
        "d174ab98d277d9f5a5611c2c9f419d9f"
    ),(
        "12345678901234567890123456789012" *
        "34567890123456789012345678901234" *
        "5678901234567890",
        "57edf4a22be3c955ac49da2e2107b67a"
    )
}
    @test md5_hash(text) == hex2bytes(hash)
end

# SHA1 hash tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/sha1-test.c
for (text, hash) in {
    (
        "",
        "da39a3ee5e6b4b0d3255bfef95601890afd80709",
    ),(
        "a",
        "86f7e437faa5a7fce15d1ddcb9eaeaea377667b8",
    ),(
        "abc",
        "a9993e364706816aba3e25717850c26c9cd0d89d",
    ),(
        "abcdefghijklmnopqrstuvwxyz",
        "32d10c7b8cf96570ca04ce37f2a19d84240d3a89",
    ),(
        "message digest",
        "c12252ceda8be8994d5fa0290a47231c1d16aae3",
    ),(
        "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmn" *
        "opqrstuvwxyz0123456789",
        "761c457bf73b14d27e9e9265c46f4b4dda11f940",
    ),(
        "1234567890123456789012345678901234567890" *
        "1234567890123456789012345678901234567890",
        "50abf5706a150990a08b2c5ea40fa0e585554732",
    ),(
        "38",
        "5b384ce32d8cdef02bc3a139d4cac0a22bb029e8"
    )
}
    @test sha1_hash(text) == hex2bytes(hash)
end

# SHA256 hash tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/sha256-test.c
for (text,hash) in {
    (
        "abc",
        "ba7816bf8f01cfea414140de5dae2223b00361a396177a9cb410ff61f20015ad"
    ),(
        "abcdbcdecdefdefgefghfghighijhijkijkljklmklmnlmnomnopnopq",
        "248d6a61d20638b8e5c026930c3e6039a33ce45964ff2167f6ecedd419db06c1"
    ),(
        "abcdefghbcdefghicdefghijdefghijkefghijklfghijklmghijklmnhijklmno" *
        "ijklmnopjklmnopqklmnopqrlmnopqrsmnopqrstnopqrstu",
        "cf5b16a778af8380036ce59e7b0492370b249b11e8f07a51afac45037afee9d1"
    ),(
        "",
        "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"
    ),(
        "a",
        "ca978112ca1bbdcafac231b39a23dc4da786eff8147c4e72b9807785afee48bb"
    ),(
        "38",
        "aea92132c4cbeb263e6ac2bf6c183b5d81737f179f21efdc5863739672f0f470"
    ),(
        "message digest",
        "f7846f55cf23e14eebeab5b4e1550cad5b509e3348fbc4efa3a1413d393cb650"
    ),(
        "abcdefghijklmnopqrstuvwxyz",
        "71c480df93d6ae2f1efad1447c66c9525e316218cf51fc8d9ed832f2daf18b73"
    ),(
        "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789",
        "db4bfcbd4da0cd85a60c3c37d3fbd8805c77f15fc6b1fdfe614ee0a7c8fdb4c0"
    ),(
        "1234567890123456789012345678901234567890123456789012345678901234" *
        "5678901234567890",
        "f371bc4a311f2b009eef952dd83ca80e2b60026c8e935592d0f9c308453c813e"
    )
}
    @test sha256_hash(text) == hex2bytes(hash)
end

# Test that digest!() actually resets the HashState
h = HashState(SHA1)
update!(h,"")
@test hexdigest!(h) == "da39a3ee5e6b4b0d3255bfef95601890afd80709"
update!(h,"The quick brown fox jumps over the lazy dog")
@test hexdigest!(h) == "2fd4e1c67a2d28fced849ee1bb76e7391b93eb12"
h = HashState(SHA256)
update!(h,"The quick brown fox jumps over the lazy dog")
@test hexdigest!(h) == "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"


text:
module MangroveBase

using Statistics
import Base.convert, Base.hash, Statistics.mean

Sub = String
Kid = String
struct ReviewInfo
  sub::Sub
  kid::Kid
end
hash(info::ReviewInfo) = hash(info.sub, hash(info.kid))

struct ReviewContent
  rating::Int
  opinion::Union{Missing, String}
end
normalize(rating::Int)::Float64 = (rating + 0.5) / 101

"Normalized version of `ReviewContent`"
ReviewSummary = Tuple{Float64, Bool}
rating(s::ReviewSummary) = s[1]

convert(::Type{ReviewSummary}, c::ReviewContent) =
  (normalize(c.rating), !ismissing(c.opinion))

Reviews = Dict{ReviewInfo, ReviewContent}
ReviewsSummary = Dict{ReviewInfo, ReviewSummary}

convert(::Type{ReviewsSummary}, r::Reviews) =
  Dict(i => convert(ReviewSummary, c) for (i, c) in r)

subs(infos::AbstractSet{ReviewInfo}) = Set(info.sub for info in infos)
kids(infos::AbstractSet{ReviewInfo}) = Set(info.kid for info in infos)
subs(infos::AbstractDict{ReviewInfo, T}) where T = subs(keys(infos))
kids(infos::AbstractDict{ReviewInfo, T}) where T = kids(keys(infos))

RATINGS = 100
check_rating(rating::Int) = @assert(0 <= rating <= RATINGS)

generate_data(subs::Int, ratings::Int)::Reviews =
  Dict(
    ReviewInfo(string(i % subs), string(i % 100)) =>
      ReviewContent((i % subs % 11) * 10, missing)
    for i in 0:ratings - 1
  )

function mean(ratings::ReviewsSummary)::Dict{Sub, Int}
  acc = Dict()
  for (info, s) in ratings
    push!(get!(acc, info.sub, []), rating(s))
  end
  out = Dict()
  for (sub, v) in acc
    out[sub] = mean(v)
  end
  out
end

mean(rs::Base.ValueIterator{ReviewsSummary}) = mean(rating(s) for s in rs)

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
"""
Bridges
`CP.Knapsack{CP.FIXED_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK}`
to `CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK}`
by creating
capacity variables.
"""
struct Knapsack2VariableCapacityKnapsackBridge{T} <: MOIBC.AbstractBridge
    capa_var::MOI.VariableIndex
    capa_con::Union{MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer}, Nothing}
    kp::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T}}
end

function MOIBC.bridge_constraint(
    ::Type{Knapsack2VariableCapacityKnapsackBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Knapsack{CP.FIXED_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T},
) where {T}
    return MOIBC.bridge_constraint(
        Knapsack2VariableCapacityKnapsackBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Knapsack2VariableCapacityKnapsackBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Knapsack{CP.FIXED_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T},
) where {T <: Integer}
    # Add the capacity variables.
    capa_var, capa_con = MOI.add_constrained_variable(model, MOI.Integer())

    # Add the capacity constraints.
    f_scalars = MOIU.scalarize(f)
    f_capa = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), capa_var)], zero(T))
    new_f = MOIU.vectorize([f_scalars..., f_capa])
    kp_set = CP.Knapsack(s.weights)
    kp = MOI.add_constraint(model, new_f, kp_set)

    return Knapsack2VariableCapacityKnapsackBridge(capa_var, capa_con, kp)
end

function MOIBC.bridge_constraint(
    ::Type{Knapsack2VariableCapacityKnapsackBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Knapsack{CP.FIXED_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T},
) where {T <: Real}
    # Add the capacity variables.
    capa_var = MOI.add_variable(model)
    capa_con = nothing

    # Add the capacity constraints.
    f_scalars = MOIU.scalarize(f)
    f_capa = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), capa_var)], zero(T))
    new_f = MOIU.vectorize([f_scalars..., f_capa])
    kp_set = CP.Knapsack(s.weights)
    kp = MOI.add_constraint(model, new_f, kp_set)

    return Knapsack2VariableCapacityKnapsackBridge(capa_var, capa_con, kp)
end

function MOI.supports_constraint(
    ::Type{Knapsack2VariableCapacityKnapsackBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Knapsack{CP.FIXED_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:Knapsack2VariableCapacityKnapsackBridge{T}}) where {T <: Integer}
    return [(MOI.Integer,)]
end

function MOIB.added_constrained_variable_types(::Type{<:Knapsack2VariableCapacityKnapsackBridge{T}}) where {T <: Real}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{Knapsack2VariableCapacityKnapsackBridge{T}}) where {T <: Integer}
    return [
        (MOI.VectorAffineFunction{T}, CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T}),
    ]
end

function MOIB.added_constraint_types(::Type{Knapsack2VariableCapacityKnapsackBridge{T}}) where {T <: Real}
    return [
        (MOI.VectorAffineFunction{T}, CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T}),
    ]
end

MOI.get(b::Knapsack2VariableCapacityKnapsackBridge, ::MOI.NumberOfVariables) = 1

function MOI.get(
    ::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex,
        MOI.Integer,
    },
) where {T <: Integer}
    return 1
end

function MOI.get(
    ::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex,
        MOI.Integer,
    },
) where {T <: Real}
    return 0
end

function MOI.get(
    ::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T},
    },
) where {T}
    return 1
end

function MOI.get(
    b::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return [b.capa_var]
end

function MOI.get(
    b::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        CP.Knapsack{CP.VARIABLE_CAPACITY_KNAPSACK, CP.UNVALUED_KNAPSACK, T},
    },
) where {T}
    return [b.kp]
end

function MOI.get(
    b::Knapsack2VariableCapacityKnapsackBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex,
        MOI.Integer,
    },
) where {T <: Integer}
    return [b.capa_con]
end
# Undefined for <: Real.


text:
# indexnotation/tensormacro.jl
#
# Defines the @tensor macro which switches to an index-notation environment.
"""
    @tensor(block)

Specify one or more tensor operations using Einstein's index notation. Indices can
be chosen to be arbitrary Julia variable names, or integers. When contracting several
tensors together, this will be evaluated as pairwise contractions in left to right
order, unless the so-called NCON style is used (positive integers for contracted
indices and negative indices for open indices).
"""
macro tensor(ex::Expr)
    tensorify(ex)
end

"""
    @tensoropt(optex, block)
    @tensoropt(block)

Specify one or more tensor operations using Einstein's index notation. Indices can
be chosen to be arbitrary Julia variable names, or integers. When contracting several
tensors together, the macro will determine (at compile time) the optimal contraction
order depending on the cost associated to the individual indices. If no `optex` is
provided, all indices are assumed to have an abstract scaling `χ` which is optimized
in the asympotic limit of large `χ`.

The cost can be specified in the following ways:

```julia
@tensoropt (a=>χ,b=>χ^2,c=>2*χ,e=>5) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost as specified for listed indices, unlisted indices have cost 1 (any symbol for χ can be used)
@tensoropt (a,b,c,e) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost χ for indices a,b,c,e, other indices (d,f) have cost 1
@tensoropt !(a,b,c,e) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# cost 1 for indices a,b,c,e; other indices (d,f) have asymptotic cost χ
@tensoropt C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost χ for all indices (a,b,c,d,e,f)
```

Note that `@tensoropt` will optimize any tensor contraction sequence it encounters
in the (block of) expressions. It will however not break apart expressions that have
been explicitly grouped with parenthesis, i.e. in
```julia
@tensoroptC[a,b,c,d] := A[a,e,c,f,h]*(B[f,g,e,b]*C[g,d,h])
```
it will always contract `B` and `C` first. For a single tensor contraction sequence,
the optimal contraction order and associated (asymptotic) cost can be obtained using
`@optimalcontractiontree`.
"""
macro tensoropt(ex::Expr)
    tensorify(ex, optdata(ex))
end
macro tensoropt(optex::Expr, ex::Expr)
    tensorify(ex, optdata(optex, ex))
end
macro optimalcontractiontree(ex::Expr)
    if isassignment(ex) || isdefinition(ex)
        _, ex = getlhsrhs(ex::Expr)
    end
    if !(ex.head == :call && ex.args[1] == :*)
        error("cannot compute optimal contraction tree for this expression")
    end
    network = [getindices(ex.args[k]) for k = 2:length(ex.args)]
    tree, cost = optimaltree(network, optdata(ex))
    return tree, cost
end
macro optimalcontractiontree(optex::Expr, ex::Expr)
    if isassignment(ex) || isdefinition(ex)
        _, ex = getlhsrhs(ex::Expr)
    end
    if !(ex.head == :call && ex.args[1] == :*)
        error("cannot compute optimal contraction tree for this expression")
    end
    network = [getindices(ex.args[k]) for k = 2:length(ex.args)]
    tree, cost = optimaltree(network, optdata(optex, ex))
    return tree, cost
end

# Process data for @tensoropt and @optimalcontractiontree
function parsecost(ex::Expr)
    if ex.head == :call && ex.args[1] == :*
        return *(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :+
        return +(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :-
        return -(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :^
        return ^(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :/
        return /(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :big
        return big(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :float
        return float(map(parsecost, ex.args[2:end])...)
    elseif ex.head == :call && ex.args[1] == :Int128
        return Int128(map(parsecost, ex.args[2:end])...)
    else
        error("invalid index cost specification: $ex")
    end
end
parsecost(ex::Number) = ex
parsecost(ex::Symbol) = Power{ex}(1,1)

function optdata(ex::Expr)
    allindices = getallindices(ex)
    cost = Power{:χ}(1,1)
    return Dict{Any, typeof(cost)}(i=>cost for i in allindices)
end

function optdata(optex::Expr, ex::Expr)
    if optex.head == :tuple
        isempty(optex.args) && return nothing
        args = optex.args
        if isa(args[1], Expr) && args[1].head == :call && args[1].args[1] == :(=>)
            indices = Vector{Any}(undef, length(args))
            costs = Vector{Any}(undef, length(args))
            costtype = typeof(parsecost(args[1].args[3]))
            for k = 1:length(args)
                if isa(args[k], Expr) && args[k].head == :call && args[k].args[1] == :(=>)
                    indices[k] = makeindex(args[k].args[2])
                    costs[k] = parsecost(args[k].args[3])
                    costtype = promote_type(costtype, typeof(costs[k]))
                else
                    error("invalid index cost specification")
                end
            end
            costs = convert(Vector{costtype}, costs)
        else
            indices = map(makeindex, args)
            costtype = Power{:χ,Int}
            costs = fill(Power{:χ,Int}(1,1), length(args))
        end
        return Dict{Any, costtype}(indices[k]=>costs[k] for k = 1:length(args))
    elseif optex.head == :call && optex.args[1] == :!
        allindices = unique(getallindices(ex))
        excludeind = map(makeindex, optex.args[2:end])
        cost = Power{:χ}(1,1)
        d = Dict{Any, typeof(cost)}(i=>cost for i in allindices)
        for i in excludeind
            d[i] = 1
        end
        return d
    else
         error("invalid index cost specification")
     end
 end

# functions for parsing and processing tensor expressions
function tensorify(ex::Expr, optdata = nothing)
    ex = expandconj(ex)
    ex = processcontractorder(ex, optdata)
    # assignment case
    if isassignment(ex) || isdefinition(ex)
        lhs, rhs = getlhsrhs(ex)
        if isa(rhs, Expr) && rhs.head == :call && rhs.args[1] == :throw
            return rhs
        end

        # process left hand side
        if istensor(lhs) && istensorexpr(rhs)
            indices = getindices(rhs)

            if lhs.head == :ref && length(lhs.args) == 2 && lhs.args[2] == :(:)
                if all(isa(i, Integer) && i < 0 for i in indices)
                    lhs = Expr(:ref, lhs.args[1], sort(indices, rev=true)...)
                else
                    error("cannot automatically infer index order of left hand side")
                end
            end

            if hastraceindices(lhs)
                err = "left hand side of an assignment should have unique indices: $lhs"
                return :(throw(IndexError($err)))
            end
            dst, leftind, rightind = maketensor(lhs)
            if Set(vcat(leftind,rightind)) != Set(indices)
                err = "non-matching indices between left and right hand side: $ex"
                return :(throw(IndexError($err)))
            end
            if isassignment(ex)
                if ex.head == :(=)
                    return deindexify(dst, false, rhs, true, leftind, rightind)
                elseif ex.head == :(+=)
                    return deindexify(dst, true, rhs, 1, leftind, rightind)
                else
                    return deindexify(dst, true, rhs, -1, leftind, rightind)
                end
            else
                return Expr(:(=), dst, deindexify(nothing, false, rhs, true, leftind, rightind, false))
            end
        elseif isassignment(ex) && isscalarexpr(lhs)
            if istensorexpr(rhs) && isempty(getindices(rhs))
                return Expr(ex.head, makescalar(lhs), Expr(:call, :scalar, deindexify(nothing, false, rhs, true, [], [], true)))
            elseif isscalarexpr(rhs)
                return Expr(ex.head, makescalar(lhs), makescalar(rhs))
            end
        else
            return ex # likely an error
        end
    end

    if ex.head == :block
        return Expr(ex.head, map(x->tensorify(x, optdata), ex.args)...)
    end
    if ex.head == :for
        return Expr(ex.head, esc(ex.args[1]), tensorify(ex.args[2], optdata))
    end
    # constructions of the form: a = @tensor ...
    if isscalarexpr(ex)
        return makescalar(ex)
    end
    if istensorexpr(ex)
        if !isempty(getindices(ex))
            err = "cannot evaluate $ex to a scalar: uncontracted indices"
            return :(throw(IndexError($err)))
        end
        ex = processcontractorder(ex, optdata)
        return Expr(:call, :scalar, deindexify(nothing, false, ex, true, [], [], true))
    end
    error("invalid syntax in @tensor macro: $ex")
end
tensorify(ex::Symbol, optdata = nothing) = esc(ex)
tensorify(ex, optdata = nothing) = ex

# expandconj: conjugate individual terms or factors instead of a whole expression
function expandconj(ex::Expr)
    if isgeneraltensor(ex) || isscalarexpr(ex)
        return ex
    elseif ex.head == :call && ex.args[1] == :conj
        @assert length(ex.args) == 2
        return conjexpr(expandconj(ex.args[2]))
    else
        return Expr(ex.head, map(expandconj, ex.args)...)
    end
end
expandconj(ex) = ex

function conjexpr(ex::Expr)
    if ex.head == :call && ex.args[1] == :conj
        return ex.args[2]
    elseif isgeneraltensor(ex) || isscalarexpr(ex)
        return Expr(:call, :conj, ex)
    elseif ex.head == :call && (ex.args[1] == :* || ex.args[1] == :+ || ex.args[1] == :-)
        return Expr(ex.head, ex.args[1], map(conjexpr, ex.args[2:end])...)
    elseif ex.head == :call && (ex.args[1] == :/ || ex.args[1] == :\)
        return Expr(ex.head, ex.args[1], map(conjexpr, ex.args[2:end])...)
    else
        error("cannot conjugate expression: $ex")
    end
end
conjexpr(ex::Number) = conj(ex)
conjexpr(ex::Symbol) = Expr(:call, :conj, ex)
conjexpr(ex) = ex

# processcontractorder: convert multi-argument multiplication into tree of pairwise multiplication
function processcontractorder(ex::Expr, optdata)
    ex = Expr(ex.head, map(e->processcontractorder(e, optdata), ex.args)...)
    if ex.head == :call && ex.args[1] == :* && length(ex.args) > 3
        args = ex.args[2:end]
        network = map(getindices, args)
        err = "invalid tensor contraction: $ex"
        for a in getallindices(ex)
            count(a in n for n in network) <= 2 || return :(throw(IndexError($err)))
        end
        if optdata === nothing
            if isnconstyle(network)
                tree = ncontree(network)
                ex = tree2expr(args, tree)
            else
                ex = Expr(:call, :*, args[1], args[2])
                for k = 3:length(args)
                    ex = Expr(:call, :*, ex, args[k])
                end
            end
        else
            tree, = optimaltree(network, optdata)
            ex = tree2expr(args, tree)
        end
    end
    return ex
end
processcontractorder(ex, optdata) = ex

function tree2expr(args, tree)
    if isa(tree, Int)
        return args[tree]
    else
        return Expr(:call, :*, tree2expr(args, tree[1]), tree2expr(args, tree[2]))
    end
end

# deindexify: parse tensor operations
function deindexify(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    if isgeneraltensor(ex)
        return deindexify_generaltensor(dst, β, ex, α, leftind, rightind, istemporary)
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-) # linear combination
        return deindexify_linearcombination(dst, β, ex, α, leftind, rightind, istemporary)
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 # multiplication: should be pairwise by now
        if isscalarexpr(ex.args[2])
            return deindexify(dst, β, ex.args[3], Expr(:call, :*, makescalar(ex.args[2]), α), leftind, rightind, istemporary)
        elseif isscalarexpr(ex.args[3])
            return deindexify(dst, β, ex.args[2], Expr(:call, :*, α, makescalar(ex.args[3])), leftind, rightind, istemporary)
        else
            return deindexify_contraction(dst, β, ex, α, leftind, rightind, istemporary)
        end
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        return deindexify(dst, β, ex.args[2], Expr(:call, :/, α, makescalar(ex.args[3])), leftind, rightind, istemporary)
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        return deindexify(dst, β, ex.args[3], Expr(:call, :\, makescalar(ex.args[2]), α), leftind, rightind, istemporary)
    end
    throw(ArgumentError("problem with parsing $ex"))
end

function deindexify_generaltensor(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    src, srcleftind, srcrightind, α2, conj = makegeneraltensor(ex)
    srcind = vcat(srcleftind, srcrightind)
    conjarg = conj ? :(:C) : :(:N)

    p1 = (map(l->_findfirst(isequal(l), srcind), leftind)...,)
    p2 = (map(l->_findfirst(isequal(l), srcind), rightind)...,)

    αsym = gensym()
    if dst === nothing
        dst = gensym()
        if istemporary
            initex = quote
                $αsym = $α*$α2
                $dst = cached_similar_from_indices($(QuoteNode(dst)), promote_type(eltype($src),typeof($αsym)), $p1, $p2, $src, $conjarg)
            end
        else
            initex = quote
                $αsym = $α*$α2
                $dst = similar_from_indices(promote_type(eltype($src),typeof($αsym)), $p1, $p2, $src, $conjarg)
            end
        end
    else
        initex = :($αsym = $α*$α2)
    end

    if hastraceindices(ex)
        traceind = unique(setdiff(setdiff(srcind,leftind), rightind))
        q1 = (map(l->_findfirst(isequal(l), srcind), traceind)...,)
        q2 = (map(l->_findlast(isequal(l), srcind), traceind)...,)
        if !isperm((p1...,p2...,q1...,q2...)) ||
            length(srcind) != length(leftind) + length(rightind) + 2*length(traceind)
            err = "trace: $(tuple(srcleftind..., srcrightind...)) to $(tuple(leftind..., rightind...)))"
            return :(throw(IndexError($err)))
        end
        return quote
            $initex
            trace!($α * $α2, $src, $conjarg, $β, $dst, $p1, $p2, $q1, $q2)
            # $dst
        end
    else
        if !isperm((p1...,p2...)) ||
            length(srcind) != length(leftind) + length(rightind)
            err = "add: $(tuple(srcleftind..., srcrightind...)) to $(tuple(leftind..., rightind...)))"
            return :(throw(IndexError($err)))
        end
        return quote
            $initex
            add!($α * $α2, $src, $conjarg, $β, $dst, $p1, $p2)
            # $dst
        end
    end
end
function deindexify_linearcombination(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    if ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-) # addition: add one by one
        if dst === nothing
            αnew = Expr(:call, :*, α, Expr(:call, :one, geteltype(ex)))
            ex1 = deindexify(dst, β, ex.args[2], αnew, leftind, rightind, istemporary)
            dst = gensym()
            returnex = :($dst = $ex1)
        else
            returnex = deindexify(dst, β, ex.args[2], α, leftind, rightind, istemporary)
        end
        αnew = (ex.args[1] == :-) ? Expr(:call, :-, α) : α
        for k = 3:length(ex.args)
            ex1 = deindexify(dst, true, ex.args[k], αnew, leftind, rightind)
            returnex = quote
                $returnex
                $ex1
            end
        end
        return quote
            $returnex
            # $dst
        end
    else
        throw(ArgumentError("unable to deindexify linear combination: $ex"))
    end
end
function deindexify_contraction(dst, β, ex::Expr, α, leftind::Vector{Any}, rightind::Vector{Any}, istemporary = false)
    @assert ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 &&
        istensorexpr(ex.args[2]) && istensorexpr(ex.args[3])
    exA = ex.args[2]
    exB = ex.args[3]

    indA = getindices(exA)
    indB = getindices(exB)
    cind = intersect(indA, indB)
    indC = vcat(leftind, rightind)
    oindA = intersect(indC, indA) # in the order they appear in C
    oindB = intersect(indC, indB) # in the order they appear in C

    if indC == vcat(oindB, oindA) # reorder
        exA, exB = exB, exA
        indA, indB = indB, indA
        oindA, oindB = oindB, oindA
    end

    symA = gensym()
    symB = gensym()
    symC = gensym()
    symTC = gensym()

    # prepare tensors or tensor expressions
    if dst === nothing
        TA = geteltype(exA)
        TB = geteltype(exB)
        TC = Expr(:call, :promote_type, TA, TB, :(typeof($α)))
    else
        TC = Expr(:call, :eltype, dst)
    end

    if !isgeneraltensor(exA) || hastraceindices(exA)
        initA = deindexify(nothing, false, exA, true, oindA, cind, true)
        poA = ((1:length(oindA))...,)
        pcA = length(oindA) .+ ((1:length(cind))...,)
        conjA = :(:N)
        initA = Expr(:(=), symA, initA)
        αA = 1
    else
        A, indlA, indrA, αA, conj = makegeneraltensor(exA)
        indA = vcat(indlA, indrA)
        poA = (map(l->_findfirst(isequal(l), indA), oindA)...,)
        pcA = (map(l->_findfirst(isequal(l), indA), cind)...,)
        TA = dst === nothing ? :(float(eltype($A))) : :(eltype($dst))
        conjA = conj ? :(:C) : :(:N)
        initA = Expr(:(=), symA, A)
    end

    if !isgeneraltensor(exB) || hastraceindices(exB)
        initB = deindexify(nothing, false, exB, true, oindB, cind, true)
        poB = ((1:length(oindB))...,)
        pcB = length(oindB) .+ ((1:length(cind))...,)
        conjB = :(:N)
        initB = Expr(:(=), symB, initB)
        αB = 1
    else
        B, indlB, indrB, αB, conj = makegeneraltensor(exB)
        indB = vcat(indlB, indrB)
        poB = (map(l->_findfirst(isequal(l), indB), oindB)...,)
        pcB = (map(l->_findfirst(isequal(l), indB), cind)...,)
        conjB = conj ? :(:C) : :(:N)
        initB = Expr(:(=), symB, B)
    end

    oindAB = vcat(oindA, oindB)
    p1 = (map(l->_findfirst(isequal(l), oindAB), leftind)...,)
    p2 = (map(l->_findfirst(isequal(l), oindAB), rightind)...,)
    if !(isperm((poA...,pcA...)) && length(indA) == length(poA)+length(pcA)) ||
        !(isperm((pcB...,poB...)) && length(indB) == length(poB)+length(pcB)) ||
        !(isperm((p1...,p2...)) && length(oindAB) == length(p1)+length(p2))
        err = "contraction: $(tuple(leftind..., rightind...)) from $(tuple(indA...,)) and $(tuple(indB...,)))"
        return :(throw(IndexError($err)))
    end
    if dst === nothing
        if istemporary
            initC = :($symC = cached_similar_from_indices($(QuoteNode(symC)), $symTC, $poA, $poB, $p1, $p2, $symA, $symB, $conjA, $conjB))
        else
            initC = :($symC = similar_from_indices($symTC, $poA, $poB, $p1, $p2, $symA, $symB, $conjA, $conjB))
        end
    else
        initC = :($symC = $dst)
    end

    return quote
        $symTC = $TC
        $initA
        $initB
        $initC
        contract!($α*$αA*$αB, $symA, $conjA, $symB, $conjB, $β, $symC,
                    $poA, $pcA, $poB, $pcB, $p1, $p2,
                    $((gensym(),gensym(),gensym())))
        # $symC
    end
end


text:
"""
    make_optimization_model(
        model::MetabolicModel,
        optimizer;
        sense = MOI.MAX_SENSE,
    )

Convert `MetabolicModel`s to a JuMP model, place objectives and the equality
constraint.
"""
function make_optimization_model(model::MetabolicModel, optimizer; sense = MOI.MAX_SENSE)
    m, n = size(stoichiometry(model))
    xl, xu = bounds(model)

    optimization_model = COBREXA.JuMP.Model(optimizer)
    @variable(optimization_model, x[i = 1:n])
    @objective(optimization_model, sense, objective(model)' * x)
    @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
    @constraint(optimization_model, lbs, xl .<= x) # lower bounds
    @constraint(optimization_model, ubs, x .<= xu) # upper bounds

    C = coupling(model) # empty if no coupling
    cl, cu = coupling_bounds(model)
    isempty(C) || @constraint(optimization_model, c_lbs, cl .<= coupling(model) * x) # coupling lower bounds
    isempty(C) || @constraint(optimization_model, c_ubs, coupling(model) * x .<= cu) # coupling upper bounds

    return optimization_model
end

"""
    optimize_model(
        model::MetabolicModel,
        optimizer;
        sense = MOI.MIN_SENSE,
    )

Use JuMP to solve an instance of CoreModel
"""
function optimize_model(model::MetabolicModel, optimizer; sense = MOI.MIN_SENSE)
    optimization_model = make_optimization_model(model, optimizer; sense = sense)
    COBREXA.JuMP.optimize!(optimization_model)
    return optimization_model
end

"""
    is_solved(optmodel)

Return `true` if `optmodel` solved successfully (solution is optimal or locally optimal).
Return `false` if any other termination status is reached. Termination status is defined
in the documentation of `JuMP`.
"""
function is_solved(optmodel)
    COBREXA.JuMP.termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ? true :
    false
end


text:
using Test
using PrettyTables
using Tables
using Markdown

println("Text backend")
println("============")
println()
include("./text_backend.jl")
println()

println("HTML backend")
println("============")
println()
include("./html_backend.jl")
println()

println("LaTeX backend")
println("=============")
println()
include("./latex_backend.jl")
println()

println("General")
println("=======")
println()
include("./general.jl")
println()


text:
module Exceptions

using SearchLight

export NotPersistedException, MissingDatabaseConfigurationException, DatabaseAdapterException
export UnretrievedModelException, InvalidModelException
export TypeConversionException, MissingConversionMethodException

abstract type SearchLightException <: Exception end

Base.showerror(io::IO, ex::T) where {T<:SearchLightException} = print(io, ex.msg)

struct NotPersistedException{T} <: Exception
  model::T
  msg
end
NotPersistedException(model) = NotPersistedException(model, "NotPersistedException: Model is not persisted \n$model")

struct MissingDatabaseConfigurationException <: Exception
  msg::String
end
MissingDatabaseConfigurationException() = MissingDatabaseConfigurationException("The database configuration can not be found")

struct NotConnectedException <: Exception
  msg::String
end
NotConnectedException() = NotConnectedException("SearchLight is not connected to the database")

struct DatabaseAdapterException <: Exception
  msg::String
end
DatabaseAdapterException() = DatabaseAdapterException("The SearchLight database adapter has thrown an unexpected exception")

struct UnretrievedModelException{T} <: Exception
  model::T
  id
  msg
end
UnretrievedModelException(model, id) = UnretrievedModelException(model, id, "UnretrievedModelException: the $(typeof(model)) data could not be retrieved for id $id. \nModel: \n$model")

struct InvalidModelException{T} <: Exception
  model::T
  errors
  msg
end
InvalidModelException(model, errors) = InvalidModelException("The $(typeof(model)) model has validation errors:\n$errors")

struct TypeConversionException{T,R,X<:Exception} <: Exception
  model::T
  field::Symbol
  value::R
  exception::X
end

struct MissingConversionMethodException{T,R} <: Exception
  type::T
  value::R
end

end


text:
macro check_if_loaded()
  quote
    if ! MKL_PARDISO_LOADED[]
      error("MKL pardiso is not properly loaded")
    end
  end
end

function pardisoinit!(
  pt::Vector{Int},
  mtype::Integer,
  iparm::Vector{Int32})

  @check_if_loaded

  ccall(
    pardisoinit_sym[],
    Cvoid, (
      Ptr{Int},
      Ptr{Int32},
      Ptr{Int32}),
    pt,
    Ref(Int32(mtype)),
    iparm)

end

function pardiso!(
  pt::Vector{Int},
  maxfct::Integer,
  mnum::Integer,
  mtype::Integer,
  phase::Integer,
  n::Integer,
  a::Vector{T},
  ia::Vector{Int32},
  ja::Vector{Int32},
  perm::Vector{Int32},
  nrhs::Integer,
  iparm::Vector{Int32},
  msglvl::Integer,
  b::Vector{T},
  x::Vector{T}) where T

  @check_if_loaded

  @assert T == pardiso_data_type(mtype,iparm)

  err = Ref(zero(Int32))

  ccall(
    pardiso_sym[],
    Cvoid, (
      Ptr{Int},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Cvoid},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Int32},
      Ptr{Cvoid},
      Ptr{Cvoid},
      Ptr{Int32}),
    pt,
    Ref(Int32(maxfct)),
    Ref(Int32(mnum)),
    Ref(Int32(mtype)),
    Ref(Int32(phase)),
    Ref(Int32(n)),
    a,
    ia,
    ja,
    perm,
    Ref(Int32(nrhs)),
    iparm,
    Ref(Int32(msglvl)),
    b,
    x,
    err)

  return Int(err[])

end

function pardiso_64!(
  pt::Vector{Int},
  maxfct::Integer,
  mnum::Integer,
  mtype::Integer,
  phase::Integer,
  n::Integer,
  a::Vector{T},
  ia::Vector{Int64},
  ja::Vector{Int64},
  perm::Vector{Int64},
  nrhs::Integer,
  iparm::Vector{Int64},
  msglvl::Integer,
  b::Vector{T},
  x::Vector{T}) where T

  @check_if_loaded

  @assert T == pardiso_data_type(mtype,iparm)

  err = Ref(zero(Int64))

  ccall(
    pardiso_64_sym[],
    Cvoid, (
      Ptr{Int},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Cvoid},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Int64},
      Ptr{Cvoid},
      Ptr{Cvoid},
      Ptr{Int64}),
    pt,
    Ref(Int64(maxfct)),
    Ref(Int64(mnum)),
    Ref(Int64(mtype)),
    Ref(Int64(phase)),
    Ref(Int64(n)),
    a,
    ia,
    ja,
    perm,
    Ref(Int64(nrhs)),
    iparm,
    Ref(Int64(msglvl)),
    b,
    x,
    err)

  return Int(err[])

end

function pardiso_getdiag!(
  pt::Vector{Int},
  df::Vector{T},
  da::Vector{T},
  mnum::Integer,
  mtype::Integer,
  iparm::Vector{<:Integer}) where T

  @assert T == pardiso_data_type(mtype,iparm)

  pardiso_getdiag!(pt,df,da,mnum)

end

function pardiso_getdiag!(
  pt::Vector{Int},
  df::Vector{T},
  da::Vector{T},
  mnum::Integer) where T

  @check_if_loaded

  err = Ref(zero(Int32))

  ccall(
   pardiso_getdiag_sym[],
   Cvoid,(
     Ptr{Int},
     Ptr{Cvoid},
     Ptr{Cvoid},
     Ptr{Int32},
     Ptr{Int32}),
   pt,
   df,
   da,
   Ref(Int32(mnum)),
   err)

  return Int(err[])

end

function pardiso_data_type(mtype::Integer,iparm::Vector{<:Integer})

  # Rules taken from
  # https://software.intel.com/en-us/mkl-developer-reference-fortran-pardiso-data-type

  T::DataType = Any

  if mtype in (1,2,-2,11)
    if iparm[28] == 0
      T = Float64
    else
      T = Float32
    end
  elseif mtype in (3,6,13,4,-4)
    if iparm[28] == 0
      T = Complex{Float64}
    else
      T = Complex{Float32}
    end
  else
    error("Unknown matrix type: mtype = $mtype")
  end

  T

end