text:
using Statistics
using Printf
using Parameters

using ..UnitCells: UnitCell
using ..Lattices: Lattice, sorted_neighbor_table_perm!, loc_to_site, calc_neighbor_table
using ..Checkerboard: checkerboard_mul!, checkerboard_transpose_mul!, checkerboard_order, checkerboard_groups
using ..IterativeSolvers: IterativeSolver, GMRES, ConjugateGradient, BiCGStab, solve!
using ..Utilities: get_index

export HolsteinModel
export assign_μ!, assign_ω!, assign_λ!, assign_ω₄!
export assign_t!, assign_ωᵢⱼ!
export get_index, get_site, get_τ
export setup_checkerboard!, update_model!
export write_phonons, read_phonons

import Base: eltype, size, length, *
import LinearAlgebra: mul!, ldiv!, transpose!
import Random: randn!

mutable struct HolsteinModel{T1,T2,T3,T4} <: AbstractModel{T1,T2,T3,T4}

    ################################
    ## COMPLETE MODEL HAMILTONIAN ##
    ################################

    # H =  ∑ Pᵢ²/2 + ∑ (ωᵢ²/2) xᵢ² [Einstein Phonons]
    #   +  ∑ λᵢ xᵢ nᵢ              [El-Ph Coupling]
    #   +  ∑ ωᵢⱼ(xᵢ ± xⱼ)²         [Phonon Dispersion]
    #   -  ∑ μᵢ nᵢ                 [Chemical Potential]
    #   -  ∑ tᵢⱼ(c⁺ᵢcⱼ + h.c.)     [Electron Kinetic Energy]

    ####################################################################
    ## CHARACTERIZING TEMPERATURE AND SIZE OF D+1 DIMENSIONAL LATTICE ##
    ####################################################################

    "inverse temperature"
    β::T1

    "discretization in imaginary time direction"
    Δτ::T1

    "length of imaginary time axis"
    Lτ::Int

    "number of sites in physical lattice"
    Nsites::Int

    "Number of bonds/hoppings in lattice"
    Nbonds::Int

    "number of phonons in lattice"
    Nph::Int

    "dimension of M matrix"
    Ndim::Int

    "number of degrees of freedom"
    Ndof::Int

    "number of type of bonds in lattice."
    nbonds::Int

    "number of type of phonons in lattice."
    nph::Int

    #######################################
    ## FOR REPRESENTING LATTICE GEOMETRY ##
    #######################################

    "represents lattice"
    lattice::Lattice{T1}

    ############################
    ## HOLSTEIN PHONON FIELDS ##
    ############################

    "phonon fields stored in the order `[x₁(1),...,x₁(Lτ),...,xₙ(1),...,xₙ(Lτ)]`
    where `n` is the number of sites in the lattice and `Lτ` is the length of the imaginary time axis."
    x::Vector{T1}

    ############################################
    ## VECTORS REPRESENTING DIAGONAL MATRICES ##
    ############################################

    "vector representing the diagonal matrix exp(-Δτ⋅V[x])"
    expnΔτV::Vector{T2}

    #############################################################
    ## SPECIFIES ELECTRON KINETIC ENERGY (TIGHT BINDING MODEL) ##
    #############################################################

    "definitions of bonds in lattice."
    bond_definitions::Vector{Bond{T2}}

    "electron hopping energies in tight binding model"
    t::Vector{T2}

    "checkerboard permutation"
    checkerboard_perm::Vector{Int}

    "cosh of electron hopping parameters in t. Sorted according to checkerboard decomposition."
    cosht::Vector{T2}

    "sinh of electron hopping parameters in t. Sorted according to checkerboard decomposition."
    sinht::Vector{T2}

    "neighboring sites in tight binding model. Sorted according to checkerboard decomposition."
    neighbor_table::Matrix{Int}

    #####################################################
    ## SPECIFIES ON-SITE ENERGIES (CHEMICAL POTENTIAL) ##
    #####################################################

    "chemical potential for each site in lattice"
    μ::Vector{T1}

    ##############################
    ## SPECIFIES HOLSTEIN MODEL ##
    ##############################

    "frequency of each phonon"
    ω::Vector{T1}

    "coefficient for anharmonic term X^4"
    ω₄::Vector{T1}

    "local electron-phonon coupling"
    λ::Vector{T1}

    "local non-linear electron-phonon coupling"
    λ₂::Vector{T1}

    ###################################
    ## SPECIFIES HOLSTEIN DISPERSION ##
    ###################################

    "extended holstein frequency of the form ωᵢⱼ(xᵢ±xⱼ)²"
    ωᵢⱼ::Vector{T1}

    "specifies which two sites i,j that are coupled in ωᵢⱼ(xᵢ±xⱼ)²"
    neighbor_table_ωᵢⱼ::Matrix{Int}

    "specifies the sign: ωᵢⱼ(xᵢ+xⱼ)² or ωᵢⱼ(xᵢ-xⱼ)²"
    sign_ωᵢⱼ::Vector{Int}

    #################################################
    ## VARIBALES FOR SOLVING M⋅x=g VIA ITERATIVELY ##
    #################################################

    "If true the default mul! routine multiplies by the M (or Mᵀ) matrix.
    If false the default mul! routine multiplies by the symmetric matrix MᵀM (or MMᵀ) instead."
    mul_by_M::Bool

    "If true multiply by Mᵀ (or MMᵀ) instead of M (or MᵀM)."
    transposed::Bool

    "A vector of length `Ndim` to temporarily store data.
    Used in multiplication by MᵀM or MMᵀ."
    v′::Vector{T2}

    "A vector of length `Ndim` to temporarily store data.
    Used to store Mᵀ⋅b result when solving M⋅x=b via MᵀM⋅x=Mᵀ⋅b."
    v″::Vector{T2}

    "A vector of length `Ndim` to temporarily store data.
    Used for calculating true residual error after linear solve."
    v‴::Vector{T2}

    """
    Iterative Solver
    """
    solver::T3

    """
    Random Number Generator
    """
    rng::T4

    """
    data folder
    """
    datafolder::String

    #######################
    ## INNER CONSTRUCTOR ##
    #######################

    """
    Constructor for Holstein type.
    """
    function HolsteinModel(lattice::Lattice{T1}, β::T1, Δτ::T1;
                           is_complex::Bool=false, iterativesolver::String="cg",rng::AbstractRNG,
                           tol::T1=1e-4, maxiter::Int=10000, restart::Int=-1) where {T1<:AbstractFloat}

        # define data type of matrix elements of M matrix
        if is_complex
            T2 = Complex{T1}
        else
            T2 = T1
        end

        # calculating length of imaginary time axis
        Lτ = round(Int,β/Δτ)

        # number of sites in lattice
        Nsites = lattice.nsites

        # number of phonons on lattice
        Nph = Nsites

        # number of degrees of freedom
        Ndof = Nph*Lτ

        # dimension of M matrix
        Ndim = Ndof

        # number of bonds in lattice
        Nbonds = 0

        # zero types of bonds initially
        nbonds = 0

        # number of types of phonons in lattice
        nph = lattice.unit_cell.norbits

        # initialize ineraction matrix, which is diagonal so stored as vector
        expnΔτV = zeros(T1,Ndim)

        # intializing phonon fields to zero
        x = zeros(T1,Ndof)

        # initializing all on-site energies to zero
        μ = zeros(T1,Nsites)

        # initializing vector to contain bond definitions
        bond_definitions = Vector{Bond{T2}}(undef,0)

        # initialize hopping parameters to empty vector
        t = Vector{T2}(undef,0)

        # initialize hopping parameters to empty vector
        cosht = Vector{T2}(undef,0)

        # initialize hopping parameters to empty vector
        sinht = Vector{T2}(undef,0)

        # initializing empty matrix to contain tight binding model neighbor_table
        neighbor_table = Matrix{Int}(undef,2,0)

        # checkerboard perumtation
        checkerboard_perm = Vector{Int}(undef,0)

        # intializing phonon frequencies to one
        ω = zeros(T1,Nph)

        # initialize electron-phonon coupling to zero
        λ = zeros(T1,Nph)

        # initialize non-linear electron-phonon coupling to zero
        λ₂ = zeros(T1,Nph)

        # initialize anharmonic X^4 term coefficient
        ω₄ = zeros(T1,Nph)

        # initizlize empty vector for inter-site phonon frequencies
        ωᵢⱼ = Vector{T1}(undef,0)

        # intialize empty matrix for storing inter-site phonon frequency neighbor_table
        neighbor_table_ωᵢⱼ = Matrix{Int}(undef,2,0)

        # intialize empty vector for sign_ωᵢⱼ
        sign_ωᵢⱼ = Vector{Int}(undef,0)

        # temporary vectors
        v′ = zeros(T2,Ndim)

        # temporary vector
        v″ = zeros(T2,Ndim)

        # temporary vector
        v‴ = zeros(T2,Ndim)

        # if true multiply by Mᵀ instead of M
        transposed = false

        # construct solver
        if lowercase(iterativesolver)=="cg"
            mul_by_M = false
            solver = ConjugateGradient(v″,tol=tol,maxiter=maxiter)
        elseif lowercase(iterativesolver)=="gmres"
            mul_by_M = true
            solver = GMRES(v″,tol=tol,restart=restart,maxiter=maxiter)
        elseif lowercase(iterativesolver)=="bicgstab"
            mul_by_M = true
            solver = BiCGStab(v″,tol=tol,maxiter=maxiter)
        end

        # get data type of iterative solver used
        T3 = typeof(solver)

        # declare random number generator
        T4  = typeof(rng)

        # constructing holstein model
        return new{T1,T2,T3,T4}(β,Δτ,Lτ,Nsites,Nbonds,Nph,Ndim,Ndof,nbonds,nph,
                                lattice, x, expnΔτV,
                                bond_definitions, t, checkerboard_perm, cosht, sinht, neighbor_table,
                                μ, ω, ω₄, λ, λ₂, ωᵢⱼ, neighbor_table_ωᵢⱼ, sign_ωᵢⱼ,
                                mul_by_M, transposed, v′, v″, v‴, solver,rng,"")
    end

end

#############################################################################
## DEFINING METHODS TO INCREMENTALLY SPECIFY THE HOLSTEIN MODEL PARAMETERS ##
#############################################################################

"""
Assign chemical potential in model for given type of orbital.
"""
function assign_μ!(holstein::HolsteinModel,μ0::T,σ0::T,orbit::Int=0) where {T<:AbstractFloat}

    if orbit==0
        R = randn(holstein.rng,holstein.Nsites)
        @. holstein.μ = μ0 + σ0 * R
    else
        for i in 1:holstein.Nsites
            if holstein.lattice.site_to_orbit[i]==orbit
                holstein.μ[i] = μ0 + σ0 * randn(holstein.rng)
            end
        end
    end

    return nothing
end

"""
Assign electron-phonon coupling in model for given type of orbital.
"""
function assign_λ!(holstein::HolsteinModel,μ0::T,σ0::T,orbit::Int=0) where {T<:AbstractFloat}

    if orbit==0
        R = randn(holstein.rng,holstein.Nph)
        @. holstein.λ = μ0 + σ0 * R
    else
        for i in 1:holstein.Nph
            if holstein.lattice.site_to_orbit[i]==orbit
                holstein.λ[i] = μ0 + σ0 * randn(holstein.rng)
            end
        end
    end

    return nothing
end

"""
Assign electron-phonon coupling in model for given type of orbital.
"""
function assign_λ₂!(holstein::HolsteinModel,μ0::T,σ0::T,orbit::Int=0) where {T<:AbstractFloat}

    if orbit==0
        R = randn(holstein.rng,holstein.Nph)
        @. holstein.λ₂ = μ0 + σ0 * R
    else
        for i in 1:holstein.Nph
            if holstein.lattice.site_to_orbit[i]==orbit
                holstein.λ₂[i] = μ0 + σ0 * randn(holstein.rng)
            end
        end
    end

    return nothing
end

"""
Assign phonon frequency in model.
"""
function assign_ω!(holstein::HolsteinModel,μ0::T,σ0::T,orbit::Int=0) where {T<:AbstractFloat}

    if orbit==0
        R = randn(holstein.rng,holstein.Nph)
        @. holstein.ω = μ0 + σ0 * R
    else
        for i in 1:holstein.Nph
            if holstein.lattice.site_to_orbit[i]==orbit
                holstein.ω[i] = μ0 + σ0 * randn(holstein.rng)
            end
        end
    end

    return nothing
end

"""
Assign coefficient for quartic phonon term.
"""
function assign_ω₄!(holstein::HolsteinModel,μ0::T,σ0::T,orbit::Int=0) where {T<:AbstractFloat}

    if orbit==0
        R = randn(holstein.rng,holstein.Nsites)
        @. holstein.ω₄ = μ0 + σ0 * R
    else
        for i in 1:holstein.Nsites
            if holstein.lattice.site_to_orbit[i]==orbit
                holstein.ω₄[i] = μ0 + σ0 * randn(holstein.rng)
            end
        end
    end

    return nothing
end

"""
Define a type of hopping/bond in the lattice.
"""
function assign_t!(holstein::HolsteinModel{T1,T2,T3}, t::T2, σt::T2, o₁::Int, o₂::Int, v::AbstractVector{Int}) where {T1,T2,T3}

    # increment number of bond definitions
    holstein.nbonds += 1

    # bond definition
    bond = Bond(t,σt,o₁,o₂,v)
    push!(holstein.bond_definitions,bond)

    # getting new neighbors accounting for periodic boundary conditions
    newneighbors = calc_neighbor_table(holstein.lattice,o₁,o₂,v)

    # getting number of new neighbors
    Nnewneighbors = size(newneighbors,2)

    # adding new neighbors to neighbor tables
    holstein.neighbor_table = hcat(holstein.neighbor_table,newneighbors)

    # phase of hopping
    phase = t/abs(t)

    # adding new hopping values
    t_new = @. phase * ( fill(abs(t),Nnewneighbors) + σt*randn(holstein.rng,Nnewneighbors) )
    append!(holstein.t, t_new)

    return nothing
end

"""
Define phonon dispersion relationship between phonons.
"""
function assign_ωᵢⱼ!(holstein::HolsteinModel{T1,T2,T3}, μω::T1, σω::T1, sgn::Int, o1::Int, o2::Int, v::AbstractVector{Int}) where {T1,T2,T3}

    # getting new neighbors accounting for periodic boundary conditions
    newneighbors = calc_neighbor_table(holstein.lattice,o1,o2,v)

    # getting number of new neighbors
    Nnewneighbors = size(newneighbors,2)

    # adding new neighbors to neighbor tables
    holstein.neighbor_table = hcat(holstein.neighbor_table,newneighbors)

    # adding new hopping values
    append!(holstein.ωᵢⱼ, fill(μω,Nnewneighbors) + σω*randn(holstein.rng,Nnewneighbors) )

    # updating neighbor table and parameter values
    assign_ωᵢⱼ!(holstein,μω,σω,o1,o2,v)

    # modifying holsteinmodel.sign_ωᵢⱼ array
    @assert abs(sgn)==1
    append!( holstein.sign_ωᵢⱼ , fill(Int,sgn,Nnewneighbors) )

    return nothing
end

##########################################################
## FUNCTIONS FOR SPECIFYING AND UPDATING HOLSTEIN MODEL ##
##########################################################

"""
    function setup_checkerboard!(holstein::HolsteinModel{T1,T2}) where {T1<:AbstractFloat,T2<:Number}

Function for sorting the hopping parameters in HolsteinModel and calculating the
cosh and sinh of the hopping parameters so that the checkerboard decomposition
is ready to go.
"""
function initialize_model!(holstein::HolsteinModel)

    if length(holstein.t)>0

        # set number of bonds in lattice
        holstein.Nbonds = length(holstein.t)

        # calculate cosh and sinh of hopping parameters
        holstein.cosht = @. cosh(holstein.Δτ*holstein.t)
        holstein.sinht = @. sinh(holstein.Δτ*holstein.t)

        # sort neighbor_table
        perm = sorted_neighbor_table_perm!(holstein.neighbor_table)
        holstein.neighbor_table .= holstein.neighbor_table[:,perm]
        holstein.cosht          .= holstein.cosht[perm]
        holstein.sinht          .= holstein.sinht[perm]

        # get checkerboard groups
        groups = checkerboard_groups(holstein.neighbor_table)

        # get checkerboard permutation
        new_perm = checkerboard_order(groups)

        # applying checkerboard permutation
        holstein.neighbor_table .= holstein.neighbor_table[:,new_perm]
        holstein.cosht          .= holstein.cosht[new_perm]
        holstein.sinht          .= holstein.sinht[new_perm]

        # record checkerboard permutation
        holstein.checkerboard_perm = sortperm(perm[new_perm])
    end

    return nothing
end

"""
    function update_model!(holstein::HolsteinModel)

Constructs the exponentiated interaction matrix for the Holstein Model
exp(-Δτ⋅V[x]) based on the current phonon fields x. Note that the matrix
exp(-Δτ⋅V[x]) is stored as a vector as it is a diagonal matrix.
"""
function update_model!(holstein::HolsteinModel)

    expnΔτV  = holstein.expnΔτV
    λ        = holstein.λ
    λ₂       = holstein.λ₂
    μ        = holstein.μ
    x        = holstein.x
    Δτ       = holstein.Δτ
    Lτ       = holstein.Lτ
    Nsites   = holstein.Nsites

    # iterating over time slices
    @inbounds @fastmath for i in 1:Nsites
        # iterating over sites in lattice
        for τ in 1:Lτ
            # getting index in vector
            index = get_index(τ,i,Lτ)
            # updating matrix element exp{-Δτ⋅Vᵢᵢ(τ)} = exp{-Δτ⋅(λᵢ⋅xᵢ(τ)-μᵢ)}
            expnΔτV[index] = exp( -Δτ * ( λ[i] * x[index] + λ₂[i] * x[index]^2 +  - μ[i] ) )
        end
    end

    return nothing
end

"""
Fill vector that will ultimate be used to update fields in model.
"""
function randn!(v::AbstractVector{T},holstein::HolsteinModel{T}) where {T}

    # fill with random numbers
    randn!(holstein.rng,v)

    return nothing
end

####################################
## MATRIX MULTIPLICATION ROUTINES ##
####################################

"""
Multiply vector by M.
"""
function mulM!(y::AbstractVector{T2},holstein::HolsteinModel{T1,T3},v::AbstractVector{T2}) where {T1<:AbstractFloat,T2<:Number,T3<:Number}

    ####################################
    ## PERFORM MULTIPLICATION y = M⋅v ##
    ####################################

    #           EXAMPLE OF M MATRIX CONVENTION
    #     |   1     0      0      0      0   +B(1) |
    #     | -B(2)   1      0      0      0     0   |
    # M = |   0   -B(3)    1      0      0     0   |
    #     |   0     0    -B(4)    1      0     0   |
    #     |   0     0      0    -B(5)    1     0   |
    #     |   0     0      0      0    -B(6)   1   |

    # Notes:
    # • y(1) = [M⋅v](1) = v(1) + B(1)⋅v(Lτ)  for τ = 1
    # • y(τ) = [M⋅v](τ) = v(τ) - B(τ)⋅v(τ-1) for τ > 1
    # • B(τ) = exp{-Δτ⋅K}⋅exp{-Δτ⋅V[x(τ)]}
    # • exp{-Δτ⋅V[x(τ)]} is the exponentiated interaction matrix and is diagonal,
    #   and as such is stored as a vector
    # • exp{-Δτ⋅K} is given by the checkerboard approximation matrix.

    @unpack Nsites, Lτ, cosht, sinht, neighbor_table, expnΔτV = holstein

    # y(τ) = exp{-Δτ⋅V[x(τ)]}⋅v(τ-1)
    for i in 1:Nsites
        for τ in 1:Lτ
            τm1  = mod1(τ-1,Lτ)
            iτ   = get_index(τ,i,Lτ)
            iτm1 = get_index(τm1,i,Lτ)
            y[iτ] = expnΔτV[iτ] * v[iτm1]
        end
    end

    # y(τ) = B(τ)⋅v(τ-1) = exp{-Δτ⋅K}⋅exp{-Δτ⋅V[x(τ)]}⋅v(τ-1)

    checkerboard_mul!(y, holstein.neighbor_table, holstein.cosht, holstein.sinht, holstein.Lτ)

    # y(τ) = v(τ) - B(τ)⋅v(τ-1)

    # iterate over sites in lattice
    @fastmath @inbounds for i in 1:Nsites

        # y(1) = v(1) + B(1)⋅v(Lτ)
        i1     = get_index(1,i,Lτ)
        y[i1] = v[i1] + y[i1]

        # iterate of time slices
        for τ in 2:Lτ

            # y(τ) = v(τ) - B(τ)⋅v(τ-1)
            iτ     = get_index(τ,i,Lτ)
            y[iτ] = v[iτ] - y[iτ]
        end
    end

    return nothing
end

"""
Multiply vector by Mᵀ.
"""
function mulMᵀ!(y::AbstractVector{T2},holstein::HolsteinModel{T1,T3},v::AbstractVector{T2}) where {T1<:AbstractFloat,T2<:Number,T3<:Number}

    #####################################
    ## PERFORM MULTIPLICATION y = Mᵀ⋅v ##
    #####################################

    #           EXAMPLE OF Mᵀ MATRIX CONVENTION
    #      |    1   -Bᵀ(2)   0      0      0      0    |
    #      |    0     1    -Bᵀ(3)   0      0      0    |
    # Mᵀ = |    0     0      1    -Bᵀ(4)   0      0    |
    #      |    0     0      0      1    -Bᵀ(5)   0    |
    #      |    0     0      0      0      1    -Bᵀ(6) |
    #      | +Bᵀ(1)   0      0      0      0      1    |

    # Notes:
    # • y(τ)  = [Mᵀ⋅v](τ)  = v(τ)  - Bᵀ(τ+1)⋅v(τ+1) for τ < Lτ
    # • y(Lτ) = [Mᵀ⋅v](Lτ) = v(Lτ) + Bᵀ(1)⋅v(1)     for τ = Lτ
    # • Bᵀ(τ) = exp{-Δτ⋅V[x(τ)]}ᵀ⋅exp{-Δτ⋅K}ᵀ
    # • exp{-Δτ⋅V[x(τ)]}ᵀ =exp{-Δτ⋅V[x(τ)]} is the exponentiated diagonal interaction
    #   matrix and is diagonal, and as such is stored as a vector
    # • [exp{-Δτ⋅K}]ᵀ is given by adjoint of the checkerboard approximation matrix.

    @unpack Nsites, Lτ, cosht, sinht, neighbor_table, expnΔτV = holstein

    # y(τ+1) = exp{-Δτ⋅K}ᵀ⋅v(τ+1)

    copyto!(y,v)
    checkerboard_transpose_mul!(y, neighbor_table, cosht, sinht, Lτ)

    # y(τ) = v(τ) - Bᵀ(τ+1)⋅v(τ+1) = v(τ) - exp{-Δτ⋅V[x(τ+1)]}⋅exp{-Δτ⋅K}ᵀ⋅v(τ+1)

    # iterate over sites in lattice
    for i in 1:Nsites

        # record y(L) = v(L) + Bᵀ(1)⋅v(1) = v(L) + exp{-Δτ⋅V[x(1)]}⋅exp{-Δτ⋅K}ᵀ⋅v(1)
        i1   = get_index(1 ,i,Lτ)
        iL   = get_index(Lτ,i,Lτ)
        y_iL = v[iL] + expnΔτV[i1] * y[i1]

        # iterate over time slices
        for τ in 1:Lτ-1

            # y(τ) = v(τ) - Bᵀ(τ+1)⋅v(τ+1) = v(τ) - exp{-Δτ⋅V[x(τ+1)]}⋅exp{-Δτ⋅K}ᵀ⋅v(τ+1)
            iτ    = get_index(τ,i,Lτ)
            iτp1  = get_index(τ+1,i,Lτ)
            y[iτ] = v[iτ] - expnΔτV[iτp1] * y[iτp1]
        end

        # y(L) = v(L) + Bᵀ(1)⋅v(1)
        y[iL] = y_iL
    end

    return nothing
end

"""
Calculates ⟨∂M/∂xᵢⱼ(τ)⟩ = uᵀ⋅[∂M/∂xᵢⱼ(τ)]⋅v for each xᵢⱼ(τ) degree of freedom,
writing each of the ⟨∂M/∂xᵢⱼ(τ)⟩ values to a vector dMdx.
"""
function muldMdx!(dMdx::AbstractVector{T2},u::AbstractVector{T2},holstein::HolsteinModel{T1,T2,T3},v::AbstractVector{T2}) where {T1,T2,T3}

    ########################################
    ## PERFORM MULTIPLICATION y = ∂M/∂x⋅v ##
    ########################################

    #           EXAMPLE OF dM/dx MATRIX CONVENTION
    #         |     0        0        0        0        0   +dB₁/dx₂ |
    #         | -dB₂/dx₂     0        0        0        0       0    |
    # dM/dx = |     0    -dB₃/dx₃     0        0        0       0    |
    #         |     0        0    -dB₄/dx₄     0        0       0    |
    #         |     0        0        0    -dB₅/dx₅     0       0    |
    #         |     0        0        0        0    -dB₆/dx₆    0    |

    # Notes:
    # • Consider y = ∂M/∂xᵢ(τ)⋅v ==>
    #
    # • ⟨∂M/∂xᵢ(τ)⟩ = uᵢ(τ)⋅-∂B/∂xᵢ(τ)⋅vᵢ(τ-1) for τ > 1
    # • ⟨∂M/∂xᵢ(1)⟩ = uᵢ(1)⋅+∂B/∂xᵢ(1)⋅vᵢ(L)   for τ = 1
    #
    # • B(τ) = exp{-Δτ⋅K}⋅exp{-Δτ⋅V[x(τ)]}
    # • ∂B/∂xᵢ(τ) = -Δτ⋅exp{-Δτ⋅K}⋅dV/dxᵢ(τ)⋅exp{-Δτ⋅V[x(τ)]}
    # • ∂B/∂xᵢ(τ) = -Δτ⋅exp{-Δτ⋅K}⋅   λᵢ    ⋅exp{-Δτ⋅V[x(τ)]}
    #
    # • Therefore the final expression is:
    # • ⟨∂M/∂xᵢ(τ)⟩ = uᵢ(τ)⋅+Δτ⋅exp{-Δτ⋅K}⋅λᵢ⋅exp{-Δτ⋅V[x(τ)]}⋅vᵢ(τ-1) for τ > 1
    # • ⟨∂M/∂xᵢ(1)⟩ = uᵢ(1)⋅-Δτ⋅exp{-Δτ⋅K}⋅λᵢ⋅exp{-Δτ⋅V[x(1)]}⋅vᵢ(L)   for τ = 1
    #
    # • Simplifying a little bit:
    # • ⟨∂M/∂xᵢ(τ)⟩ = uᵢ(τ)⋅+Δτ⋅λᵢ⋅B(τ)⋅vᵢ(τ-1) for τ > 1
    # • ⟨∂M/∂xᵢ(1)⟩ = uᵢ(1)⋅-Δτ⋅λᵢ⋅B(1)⋅vᵢ(L)   for τ = 1

    @unpack Nsites, Lτ, Δτ, cosht, sinht, neighbor_table, expnΔτV, λ, λ₂, x = holstein

    y = holstein.v′

    # ⟨∂M/∂xᵢ(τ)⟩ = +Δτ⋅λᵢ⋅exp{-Δτ⋅V[x(τ)]}⋅v(τ-1)

    # iterating over sites in lattice
    @fastmath @inbounds for i in 1:Nsites

        # ⟨∂M/∂xᵢⱼ(τ)⟩ = -Δτ⋅λᵢ⋅exp{-Δτ⋅V[x(τ)]}⋅v(1) for τ=1
        idx_1    = get_index(1,  i,  Lτ)
        idx_L    = get_index(Lτ, i,  Lτ)
        dMdx[idx_1] = -Δτ * (λ[i] + 2*λ₂[i]*x[idx_1]) * expnΔτV[idx_1] * v[idx_L]

        # iterating over time slices
        for τ in 2:Lτ

            # ⟨∂M/∂xᵢⱼ(τ)⟩ = +Δτ⋅λᵢ⋅exp{-Δτ⋅V[x(τ)]}⋅v(τ-1) for τ>1
            idx_τ   = get_index(τ,   i, Lτ)
            idx_τm1 = get_index(τ-1, i, Lτ)
            dMdx[idx_τ] = Δτ * (λ[i] + 2*λ₂[i]*x[idx_τ]) * expnΔτV[idx_τ] * v[idx_τm1]
        end
    end

    # y = exp{-Δτ⋅K}ᵀ⋅u ⟹ yᵀ = uᵀ⋅exp{-Δτ⋅K}
    copyto!(y,u)
    checkerboard_transpose_mul!(y, neighbor_table, cosht, sinht, Lτ)

    # finalize calculation of ⟨∂M/∂xᵢⱼ(τ)⟩ = uᵀ⋅[∂M/∂xᵢⱼ(τ)]⋅v
    @. dMdx = y * dMdx

    return nothing
end

##############################
## PHONON FIELD IO ROUTINES ##
##############################

"""
Writes the current phonon field configuration for a HolsteinModel to file.
"""
function write_phonons!(holstein::HolsteinModel,filename::String)

    # get lattice associated with holstein model
    lattice = holstein.lattice

    # get phonon fields
    x = holstein.x

    # open file
    open(filename,"w") do file

        # write header to file
        write(file, "L3 L2 L1 orbit tau x\n")

        # iterate over unit cells
        for l3 in 0:lattice.L3-1
            for l2 in 0:lattice.L2-1
                for l1 in 0:lattice.L1-1

                    # iterate over orbitals/sites in each unit cell
                    for orbit in 1:lattice.unit_cell.norbits

                        # get site in lattice
                        site = loc_to_site(lattice,orbit,l1,l2,l3)

                        # iterate of time slice
                        for τ in 1:holstein.Lτ

                            # get index
                            i = get_index(τ, site, holstein.Lτ)

                            # get phonon field for site
                            xi = x[i]

                            # write to file
                            write(file, @sprintf("%d %d %d %d %d %.6f\n",l3,l2,l1,orbit,τ,xi))
                        end
                    end
                end
            end
        end
    end

    return nothing
end

"""
Read phonon config from file.
"""
function read_phonons!(holstein::HolsteinModel{T1,T2,T3},filename::String) where {T1<:AbstractFloat,T2<:Number,T3<:IterativeSolver}

    # open file
    open(filename,"r") do file

        # read in header
        header = readline(file)

        # iterate over lines in file
        for line in eachline(file)

            # split line at white space
            atoms = split(line," ")

            # extract info about location
            τ     = parse(Int,atoms[5])
            orbit = parse(Int,atoms[4])
            l1    = parse(Int,atoms[3])
            l2    = parse(Int,atoms[2])
            l3    = parse(Int,atoms[1])

            # get phonon field value
            xi = parse(T1,atoms[6])

            # map location on site
            site = loc_to_site(holstein.lattice, orbit, l1, l2, l3)

            # map site and τ onto index
            i = get_index(τ, site, holstein.Lτ)

            # assign phonon value
            holstein.x[i] = xi
        end

    end

    # construct exponentiated interaction matrix
    update_model!(holstein)

    return nothing
end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
module Game2048Core

export Dirs, left, right, up, down, move!, move, Bitboard

@enum Dirs left up right down

const DIRS = (left, up, right, down)

# Vector{Dirs}(init_len) = [left for _ in 1:init_len]

# Do NOT delete; this is used to generate the LEFT RIGHT movement lookup tables
include("uint_board/move-up.jl")
include("uint_board/move-board.jl")
include("uint_board/sim-game.jl")

include("rotate-mirror.jl")

include("bitboard/make_lookups.jl")
include("bitboard/bitboard.jl")
include("bitboard/simulate_bb.jl")
# include("bitboard/playahead.jl")
# include("bitboard/fingerprint.jl")

end

ext: jl


text:
# udat.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files udat.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
"""
module UDAT
const NONE     = Int32(-1)
const FULL     = Int32(0)
const LONG     = Int32(1)
const MEDIUM   = Int32(2)
const SHORT    = Int32(3)
const RELATIVE = Int32(128)
end # module UDAT

export UDAT, UDateFormat

macro libdat(s)     ; _libicu(s, iculibi18n, "udat_")     ; end

mutable struct UDateFormat
    ptr::Ptr{Cvoid}

    function UDateFormat(tstyle::Integer, dstyle::Integer, tz::Ptr{UInt16}, len)
        err = Ref{UErrorCode}(0)
        p = ccall(@libdat(open), Ptr{Cvoid},
                  (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32,
                   Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  tstyle, dstyle, locale[], tz, len, C_NULL, 0, err)
        FAILURE(err[]) && error("ICU: $(err[]), bad date format")
        self = new(p)
        finalizer(self, close)
        self
    end

    function UDateFormat(pattern::Ptr{UInt16}, patlen, tz::Ptr{UInt16}, tzlen)
        err = Ref{UErrorCode}(0)
        p = ccall(@libdat(open), Ptr{Cvoid},
                  (Int32, Int32, Ptr{UInt8}, Ptr{UChar}, Int32,
                   Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  -2, -2, locale[], tz, tzlen,
                  pattern, patlen, err)
        FAILURE(err[]) && error("ICU: $(err[]), bad date format")
        self = new(p)
        finalizer(self, close)
        self
    end
end

UDateFormat(tstyle::Integer, dstyle::Integer, tz::WordStrings) =
    @preserve tz UDateFormat(tstyle, dstyle, pointer(tz), ncodeunits(tz))
UDateFormat(pattern::WordStrings, tz::WordStrings) =
    @preserve pattern tz UDateFormat(pointer(pattern), ncodeunits(pattern),
                                          pointer(tz), ncodeunits(tz))

UDateFormat(tstyle::Integer, dstyle::Integer, tz::Vector{UInt16}) =
    @preserve tz UDateFormat(tstyle, dstyle, pointer(tz), length(tz))
UDateFormat(pattern::Vector{UInt16}, tz::Vector{UInt16}) =
    @preserve pattern tz UDateFormat(pointer(pattern), length(pattern),
                                          pointer(tz), length(tz))

UDateFormat(pattern::AbstractString, tz::AbstractString) =
    UDateFormat(cvt_utf16(pattern), cvt_utf16(tz))

UDateFormat(tstyle::Integer, dstyle::Integer, tz::AbstractString) =
    UDateFormat(tstyle, dstyle, cvt_utf16(tz))

close(df::UDateFormat) =
    df.ptr == C_NULL || (ccall(@libdat(close), Cvoid, (Ptr{Cvoid},), df.ptr) ; df.ptr = C_NULL)

function format(df::UDateFormat, millis::Float64)
    err = Ref{UErrorCode}(0)
    buflen = 64
    buf, pnt = _allocate(UInt16, buflen)
    len = ccall(@libdat(format), Int32,
                (Ptr{Cvoid}, Float64, Ptr{UChar}, Int32, Ptr{Cvoid}, Ptr{UErrorCode}),
                df.ptr, millis, pnt, buflen, C_NULL, err)
    FAILURE(err[]) && error("ICU: $(err[]), failed to format time")
    Str(UTF16CSE, buf[1:len])
end

parse(df::UDateFormat, s::AbstractString) = parse(df, cvt_utf16(s))
parse(df::UDateFormat, s::Vector{UInt16}) = @preserve s parse(df, pointer(s), length(s))
parse(df::UDateFormat, s::WordStrings)    = @preserve s parse(df, pointer(s), ncodeunits(s))

function parse(df::UDateFormat, s16::Ptr{UInt16}, slen)
    err = Ref{UErrorCode}(0)
    ret = ccall(@libdat(parse), Float64,
                (Ptr{Cvoid}, Ptr{UChar}, Int32, Ptr{Int32}, Ptr{UErrorCode}),
                df.ptr, s16, slen, C_NULL, err)
    FAILURE(err[]) && error("ICU: $(err[]), failed to parse string")
    ret
end


text:
using JuLIP, SKTB

at = bulk(:Si)
tb = SKTB.Kwon.KwonTBModel(; potential = SKTB.ZeroT(),
                             bzquad = SKTB.MPGrid(at, (10,10,10)) )
energy(tb, at)
F0 = defm(at)
ss = range(1.0, 1.4, length = 100)
EE = Float64[]
for s in ss
   F = s * F0
   set_defm!(at, F, updatepositions=true)
   push!(EE, energy(tb, at))
end

using Plots
plot(ss, EE)

findmin(EE)
ss[67]


text:
module BallroomSkatingSystem

using Reexport
@reexport using DataFrames
using Statistics

include("helper_functions.jl")
include("skating_single_dance.jl")
include("skating_combined.jl")

export skating_single_dance, skating_combined

end


text:
using EngEconomics, Roots, Plots

# Given
aInit = -350
aAnnual = -90
aLife = 10

bInit = -120
bAnnual = -100
bLife = 4

MARR = 0.05

# Find IRR
ΔInit_eq(x) = aInit * (1 + presentWorthFactor(x, aLife)) - bInit * (1 + presentWorthFactor(x, bLife) + presentWorthFactor(x, 2 * bLife) + presentWorthFactor(x, 3 * bLife) + presentWorthFactor(x, 4 * bLife))
ΔAnnual_eq(x) = aAnnual - bAnnual
IRR_eq(x) = ΔInit_eq(x) + ΔAnnual_eq(x) * seriesPresentAmountFactor(x, 20)
xVec = collect(0:0.01:1)
yVec = zeros(size(xVec))
plot(xVec, IRR_eq.(xVec))
plot!(xVec, yVec)

find_zero(IRR_eq, MARR)


text:
##using Markdown # hide

# ## Using the `SimilaritySearch` package
# By Eric S. Téllez
#
# This is a small tutorial showing a minimum example for working with `SimilaritySearch`
# it accepts several options that are let to defaults. While this should be enough for many
# purposes, you are invited to see the rest of the tutorials to take advantage of other
# features.
ENV["GKSwstype"] = "nul" # hide

using SimilaritySearch
# MatrixDatabase is a required wrapper that tells `SimilaritySearch` how to access underlying objects
# since it can support different kinds of objects. In this setup, each column is an object
# and will be accessed through views using the MatrixDatabase. Since the backend doesn't support
# appends or pushes, the index can be seen as an static index.

db = MatrixDatabase(randn(2, 10^5))
dist = L2Distance() # squared L2

# it can use any distance function described in `SimilaritySearch` and `Distances.jl`,
# and in fact any `SemiMetric` as described in the later package.
# The index construction is made as follows
G = SearchGraph(; dist, db)
index!(G)
#- callbacks = SearchGraphCallbacks(hyperparameters=OptimizeParameters(MinRecall(0.9)))
#- optimize!(G; callbacks)

# this will display a lot of information in the console, since as construction advances
# the hyperparameters of the index are adjusted. The default optimization takes into account
# both quality and speed, and tries to adjust to take the best of both.
# See `ParetoRecall` in docs.
#
# Once the index is created, the index can solve nearest neighbor queries

Q = MatrixDatabase(randn(2, 100))
k = 30
I, D = searchbatch(G, Q, k)
@show typeof(I) => size(I)
@show typeof(D) => size(D)
# `I` is a matrix of identifiers in `db`. Each column stores the $k$ nearest neighbors (approx.)
# for the $i$-th colum (object) in `Q`. The matrix `D` stores the the corresponding distances for
# each identifier in `I`.
#
# ## Visualizing what we just did

using Plots
@views scatter(db.matrix[1, :], db.matrix[2, :], fmt=:png, size=(600, 600), color=:cyan, ma=0.3, a=0.3, ms=1, msw=0, label="")
for c in eachcol(I)
    R = db.matrix[:, c]
    @views scatter!(R[1, :], R[2, :], m=:diamond, ma=0.3, a=0.3, color=:auto, ms=2, msw=0, label="")
end

@views scatter!(Q.matrix[1, :], Q.matrix[2, :], color=:black, m=:star, ma=0.5, a=0.5, ms=4, msw=0, label="")

savefig(joinpath(@OUTPUT, "fig-2d-t1.png")) # hide

# \fig{fig-2d-t1.png}
#
# Cyan points identify the dataset while starts are query points. The nearest neighbor points
# are colored automatically and can repeat, but they come quite close to query points,
# in dense areas they are even hidding them.

#- Using BeutifulMakie style
# ## Dependencies
using Pkg # hide
Pkg.status(["SimilaritySearch", "Plots"]) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# This file was generated, do not modify it. # hide
X = reshape([1,2,3,4], 2, 2) |> permutedims


text:
using VectorEngine, Test

@testset "vadd float64" begin

    dims = (1123,)             # intentionally an odd number
    a = round.(rand(Float64, dims) * 100)
    b = round.(rand(Float64, dims) * 100)
    c = similar(a)

    d_a = VEArray(a)
    d_b = VEArray(b)
    d_c = VEArray(c)

    function add_float64(a, b, c)
        @vectorize for i = 1:prod(size(a))
            @inbounds c[i] = a[i] + b[i]
        end
        return
    end

    @synchronize @veda add_float64(d_a, d_b, d_c)
    c = Array(d_c)
    @test a+b ≈ c

    out = sprint(io->(@device_code_native io=io @veda add_float64(d_a, d_b, d_c)))

    # do we see vld, vst, vfadd in the output?
    @test occursin("vld", out)
    @test occursin("vst", out)
    @test occursin("vfadd", out)

end

@testset "vadd float32" begin

    dims = (1123,)             # intentionally an odd number
    a = round.(rand(Float32, dims) * 100)
    b = round.(rand(Float32, dims) * 100)
    c = similar(a)

    d_a = VEArray(a)
    d_b = VEArray(b)
    d_c = VEArray(c)

    function add_float32(a, b, c)
        @vectorize length=512 for i = 1:prod(size(a))
            @inbounds c[i] = a[i] + b[i]
        end
        return
    end

    @synchronize @veda add_float32(d_a, d_b, d_c)
    c = Array(d_c)
    @test a+b ≈ c

    out = sprint(io->(@device_code_native io=io @veda add_float32(d_a, d_b, d_c)))

    # do we see vld, vst, pvfadd in the output?
    @test occursin("vld", out)
    @test occursin("vst", out)
    @test occursin("pvfadd", out)
    @test occursin("512", out)

end


text:
"""
Initialize options
"""
function create_result(data, fixedvariables, outsample, criteria, ttest, modelavg, residualtest, orderresults)

    if outsample == nothing
        outsample = 0
    end

    if (outsample isa Array && size(outsample, 1) > 0) || (!(outsample isa Array) && outsample > 0)
        push!(criteria, :rmseout)
    end

    criteria = unique(criteria)

    datanames = create_datanames(data, criteria, ttest, modelavg, residualtest)

    if modelavg
        modelavg_datanames = []
    else
        modelavg_datanames = nothing
    end

    outsample_max = data.nobs - INSAMPLE_MIN - size(data.expvars, 1) - if (data.intercept) 1 else 0 end

    if isa(outsample, Int) && outsample_max <= outsample
        outsample = 0
    end

    return AllSubsetRegressionResult(
        datanames,
        modelavg_datanames,
        fixedvariables,
        outsample,
        criteria,
        modelavg,
        ttest,
        residualtest,
        orderresults
    )
end

"""
Constructs the datanames array for results based on this structure.
    - Index
    - Covariates
        * b
        * bstd
        * T-test
    - Equation general information merged with criteria user-defined options
    - Order index from user combined criteria
    - Weight
"""
function create_datanames(data, criteria, ttest, modelavg, residualtest)

    datanames = []

    push!(datanames, INDEX)

    for expvar in data.expvars
        push!(datanames, Symbol(string(expvar, "_b")))
        if ttest
            push!(datanames, Symbol(string(expvar, "_bstd")))
            push!(datanames, Symbol(string(expvar, "_t")))
        end
    end

    testfields = (residualtest != nothing && residualtest) ? ((data.time != nothing) ? RESIDUAL_TESTS_TIME : RESIDUAL_TESTS_CROSS) : []
    general_information_criteria = unique([ EQUATION_GENERAL_INFORMATION; criteria; testfields ])
    datanames = vcat(datanames, general_information_criteria)

    push!(datanames, ORDER)
    if modelavg
        push!(datanames, WEIGHT)
    end

    return datanames
end

"""
Get insample data view
"""
function get_insample_subset(depvar_data, expvars_data, outsample, selected_variables_index)
    depvar_view = nothing
    expvars_view = nothing
    if isa(outsample, Array)
        insample = setdiff(1:size(depvar_data, 1), outsample)
        depvar_view = depvar_data[insample, 1]
        expvars_view = expvars_data[insample, selected_variables_index]
    else
        depvar_view = depvar_data[1:end-outsample, 1]
        expvars_view = expvars_data[1:end-outsample, selected_variables_index]
    end
    return depvar_view, expvars_view
end

"""
Get outsample data view
"""
function get_outsample_subset(depvar_data, expvars_data, outsample, selected_variables_index)
    depvar_view = nothing
    expvars_view = nothing
    if isa(outsample, Array)
        depvar_view = depvar_data[outsample, 1]
        expvars_view = expvars_data[outsample, selected_variables_index]
    else
        depvar_view = depvar_data[end-outsample+1:end, 1]
        expvars_view = expvars_data[end-outsample+1:end, selected_variables_index]
    end
    return depvar_view, expvars_view
end

"""
Sort rows
"""
function sortrows(B::AbstractMatrix,cols::Array; kws...)
    for i = 1:length(cols)
        if i == 1
            p = sortperm(B[:,cols[i]]; kws...)
            B = B[p,:]
        else
            i0_old = 0
            i1_old = 0
            i0_new = 0
            i1_new = 0
            for j = 1:size(B,1)-1
                if B[j,cols[1:i-1]] == B[j+1,cols[1:i-1]] && i0_old == i0_new
                    i0_new = j
                elseif B[j,cols[1:i-1]] != B[j+1,cols[1:i-1]] && i0_old != i0_new && i1_new == i1_old
                    i1_new = j
                elseif i0_old != i0_new && j == size(B,1)-1
                    i1_new = j+1
                end
                if i0_new != i0_old && i1_new != i1_old
                    p = sortperm(B[i0_new:i1_new,cols[i]]; kws...)
                    B[i0_new:i1_new,:] = B[i0_new:i1_new,:][p,:]
                    i0_old = i0_new
                    i1_old = i1_new
                end
            end
        end
    end
    return B
end

function get_varnames(datanames)
    map(h -> chop("$h", tail=2), filter(s -> endswith("$s", "_b"), datanames))
end

"""
Add values to extras
"""
function addextras(data, result)
    data.extras[ModelSelection.generate_extra_key(ALLSUBSETREGRESSION_EXTRAKEY, data.extras)] = Dict(
        :datanames => result.datanames,
        :depvar => data.depvar,
        :expvars => data.expvars,
        :nobs => data.nobs,
        :time => data.time,
        :residualtest => result.residualtest,
        :criteria => result.criteria,
        :intercept => data.intercept,
        :ttest => result.ttest,
        :outsample => result.outsample,
        :modelavg => result.modelavg,
        :fixedvariables => result.fixedvariables,
        :orderresults => result.orderresults
    )
    return data
end


text:
mutable struct Optimizer <: AbstractSampledOptimizer
    instance_optimizer
    optimizer_params::Dict{MOI.AbstractOptimizerAttribute, Any}
    parameters::SAAParameters{Float64}
    status::MOI.TerminationStatusCode
    algorithm::Union{SampleAverageApproximation, Nothing}

    function Optimizer(; optimizer = nothing,
                       kw...)
        return new(optimizer,
                   Dict{MOI.AbstractOptimizerAttribute, Any}(),
                   SAAParameters{Float64}(; kw...),
                   MOI.OPTIMIZE_NOT_CALLED,
                   nothing)
    end
end

function load_model!(optimizer::Optimizer,
                     model::StochasticModel,
                     sampler::AbstractSampler,
                     x₀::AbstractVector)
    instance_optimizer = MOI.get(optimizer, InstanceOptimizer())
    # Create new SAA algorithm
    optimizer.algorithm = SampleAverageApproximation(model,
                                                     sampler,
                                                     x₀;
                                                     optimizer = instance_optimizer,
                                                     type2dict(optimizer.parameters)...)
    return nothing
end

function MOI.optimize!(optimizer::Optimizer)
    if optimizer.algorithm === nothing
        error("SAA algorithm has not been loaded properly")
    end
    optimizer.status = optimizer.algorithm()
    return nothing
end

function optimizer_name(optimizer::Optimizer)
    return "Sample average approximation"
end

function optimal_instance(optimizer::Optimizer)
    if optimizer.algorithm === nothing
        error("SAA algorithm has not been loaded properly")
    end
    return instance(optimizer.algorithm)
end

# MOI #
# ========================== #
function MOI.get(optimizer::Optimizer, ::MOI.Silent)
    return !MOI.get(optimizer, MOI.RawParameter("log"))
end

function MOI.set(optimizer::Optimizer, attr::MOI.Silent, flag::Bool)
    MOI.set(optimizer, MOI.RawParameter("log"), !flag)
    optimizer.optimizer_params[attr] = flag
    return nothing
end

function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(SAAParameters))
        error("Unrecognized parameter name: $(name).")
    end
    return getfield(optimizer.parameters, name)
end

function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(SAAParameters))
        error("Unrecognized parameter name: $(name).")
    end
    setfield!(optimizer.parameters, name, value)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::RelativeTolerance)
    return MOI.get(optimizer, MOI.RawParameter("tolerance"))
end

function MOI.set(optimizer::Optimizer, ::RelativeTolerance, limit::Real)
    MOI.set(optimizer, MOI.RawParameter("tolerance"), limit)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::Confidence)
    return MOI.get(optimizer, MOI.RawParameter("confidence"))
end

function MOI.set(optimizer::Optimizer, ::Confidence, confidence::Real)
    MOI.set(optimizer, MOI.RawParameter("confidence"), confidence)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumSamples)
    if optimizer.algorithm != nothing
        return optimizer.algorithm.data.sample_size
    end
    return MOI.get(optimizer, MOI.RawParameter("num_samples"))
end

function MOI.set(optimizer::Optimizer, ::NumSamples, num_samples::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_samples"), num_samples)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumEvalSamples)
    return MOI.get(optimizer, MOI.RawParameter("num_eval_samples"))
end

function MOI.set(optimizer::Optimizer, ::NumEvalSamples, num_eval_samples::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_eval_samples"), num_eval_samples)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumEWSSamples)
    return MOI.get(optimizer, MOI.RawParameter("num_ews_samples"))
end

function MOI.set(optimizer::Optimizer, ::NumEWSSamples, num_ews_samples::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_ews_samples"), num_ews_samples)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumEEVSamples)
    return MOI.get(optimizer, MOI.RawParameter("num_eev_samples"))
end

function MOI.set(optimizer::Optimizer, ::NumEEVSamples, num_eev_samples::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_eev_samples"), num_eev_samples)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumLowerTrials)
    return MOI.get(optimizer, MOI.RawParameter("num_lower_trials"))
end

function MOI.set(optimizer::Optimizer, ::NumLowerTrials, num_lower_trials::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_lower_trials"), num_lower_trials)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::NumUpperTrials)
    return MOI.get(optimizer, MOI.RawParameter("num_upper_trials"))
end

function MOI.set(optimizer::Optimizer, ::NumUpperTrials, num_upper_trials::Integer)
    MOI.set(optimizer, MOI.RawParameter("num_upper_trials"), num_upper_trials)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::InstanceOptimizer)
    if optimizer.instance_optimizer === nothing
        error("Instance optimizer not set. Consider setting `InstanceOptimizer` attribute.")
    end
    return () -> begin
        instance_opt = optimizer.instance_optimizer()
        for (attr, value) in optimizer.optimizer_params
            MOI.set(instance_opt, attr, value)
        end
        try
            MOI.set(instance_opt, MOI.RawParameter("log"), optimizer.parameters.log)
            MOI.set(instance_opt, MOI.RawParameter("keep"), optimizer.parameters.keep)
            MOI.set(instance_opt, MOI.RawParameter("offset"), optimizer.parameters.offset + 1)
            MOI.set(instance_opt, MOI.RawParameter("indent"), 2 * optimizer.parameters.indent)
        catch
        end
        return instance_opt
    end
end

function MOI.set(optimizer::Optimizer, ::InstanceOptimizer, optimizer_constructor)
    optimizer.instance_optimizer = optimizer_constructor
    return nothing
end

function MOI.get(optimizer::Optimizer, ::InstanceOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute)
    if !haskey(optimizer.optimizer_params, attr)
        error("Instance optimizer attribute $(typeof(attr)) has not been set.")
    end
    return optimizer.optimizer_params[attr]
end

function MOI.set(optimizer::Optimizer, ::InstanceOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute, value)
    optimizer.optimizer_params[attr] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, param::RawInstanceOptimizerParameter)
    moi_param = MOI.RawParameter(param.name)
    if !haskey(optimizer.optimizer_params, moi_param)
        error("Instance optimizer attribute $(param.name) has not been set.")
    end
    return optimizer.optimizer_params[moi_param]
end

function MOI.set(optimizer::Optimizer, param::RawInstanceOptimizerParameter, value)
    # Only set parameter if instance optimizer supports it
    if !MOI.supports(MOI.get(optimizer, InstanceOptimizer())(), param)
        return nothing
    end
    moi_param = MOI.RawParameter(param.name)
    optimizer.optimizer_params[moi_param] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    return optimizer.status
end

function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, index::MOI.VariableIndex)
    if optimizer.algorithm === nothing
        error("SAA algorithm has not been loaded properly")
    end
    return decision(optimizer.algorithm, index)
end

function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue)
    if optimizer.algorithm === nothing
        error("SAA algorithm has not been loaded properly")
    end
    return optimizer.algorithm.data.interval
end

function MOI.is_empty(optimizer::Optimizer)
    return optimizer.algorithm === nothing
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.supports(::Optimizer, ::MOI.RawParameter) = true
MOI.supports(::Optimizer, ::AbstractSampledOptimizerAttribute) = true


text:
"""
    HintikkaXOData()

Retrieves the HintikkaXO MultiAssayExperiment, consisting of the microbiota,
metabolite and biomarker SummarizedExperiment data, by Hintikka et al.
(2021), https://doi.org/10.3390/ijerph18084049.
"""
function HintikkaXOData()

    # ##### MICROBIOTA EXPERIMENT #####

    # # read in csv file on microbiota as DataFrame
    # counts = CSV.File(joinpath(@__DIR__, "assets/counts.csv")) |> DataFrame
    # # store rowdata
    # rowdata_microbiota = DataFrame(
    #     name = counts[:, 1]
    # )
    # # store coldata
    # coldata_microbiota = DataFrame(
    #     name = names(counts)[2:end]
    # )
    # # create counts assay
    # assay_microbiota = OrderedDict{String, AbstractArray}("counts" => counts[:, 2:end] |> Matrix)

    # ##### METABOLITES EXPERIMENT #####

    # # read in csv file on metabolites as DataFrame
    # nmr = CSV.File(joinpath(@__DIR__, "assets/nmr.csv")) |> DataFrame
    # # store rowdata
    # rowdata_metabolites = DataFrame(
    #     name = nmr[:, 1]
    # )
    # # store coldata
    # coldata_metabolites = DataFrame(
    #     name = names(nmr)[2:end]
    # )
    # # create nmr assay
    # assay_metabolites = OrderedDict{String, AbstractArray}("nmr" => nmr[:, 2:end] |> Matrix)

    # ##### BIOMARKERS EXPERIMENT #####

    # # read in csv file on biomarkers as DataFrame
    # signals = CSV.File(joinpath(@__DIR__, "assets/signals.csv")) |> DataFrame
    # # store rowdata
    # rowdata_biomarkers = DataFrame(
    #     name = signals[:, 1]
    # )
    # # store coldata
    # coldata_biomarkers = DataFrame(
    #     name = names(signals)[2:end]
    # )
    # # create signals assay
    # assay_biomarkers = OrderedDict{String, AbstractArray}("signals" => signals[:, 2:end] |> Matrix)

    # ##### EXPERIMENTLIST #####

    # # create OrderedDict object to store SummarizedExperiment objects
    # expo = OrderedDict{String, SummarizedExperiment}()
    # # store first experiment as microbiota
    # expo["microbiota"] = SummarizedExperiment(assay_microbiota, rowdata_microbiota, coldata_microbiota)
    # # store second experiment as metabolites
    # expo["metabolites"] = SummarizedExperiment(assay_metabolites, rowdata_metabolites, coldata_metabolites)
    # # store third experiment as biomarkers
    # expo["biomarkers"] = SummarizedExperiment(assay_biomarkers, rowdata_biomarkers, coldata_biomarkers)

    # ##### SAMPLE DATA #####

    # # read in csv file on sample data as DataFrame
    # sample_data = CSV.File(joinpath(@__DIR__, "assets/sample_data.csv")) |> DataFrame
    # # rename cols to proper names
    # rename!(sample_data, Dict(:Sample => "name"))

    # ##### SAMPLE MAP ######

    # # read in csv file on sample map as DataFrame
    # sample_map = CSV.File(joinpath(@__DIR__, "assets/sample_map.csv")) |> DataFrame
    # # rename cols to proper names
    # rename!(sample_map, Dict(:assay => "experiment"))
    # rename!(sample_map, Dict(:primary => "sample"))
    # # put cols in proper order
    # select!(sample_map, [:sample, :experiment, :colname])

    # ##### MULTIASSAY EXPERIMENT #####

    # # construct MultiAssayExperiment object
    # MultiAssayExperiment(expo, sample_data, sample_map)

    experiment_files = [joinpath(@__DIR__, "assets/XO_microbiota_assays.csv"),
                        joinpath(@__DIR__, "assets/XO_metabolites_assays.csv"),
                        joinpath(@__DIR__, "assets/XO_biomarkers_assays.csv")]

    experiment_names = ["microbiota", "metabolites", "biomarkers"]

    import_mae_from_csv(experiment_files,
        joinpath(@__DIR__, "assets/XO_sample_data.csv"),
        joinpath(@__DIR__, "assets/XO_sample_map.csv"),
        experiment_names = experiment_names)

end

"""
    OKeefeDSData()

Retrieves the OKeefeDS SummarizedExperiment, consisting of the microbiota,
contains microbiome data from a study with African and African American groups
undergoing a two-week diet swap, by O'Keefe et al.
(2015), https://dx.doi.org/10.1038/ncomms7342.
"""
function OKeefeDSData()

    # DS_counts = CSV.File(joinpath(@__DIR__, "assets/DS_counts.csv")) |> DataFrame
    # DS_assays = OrderedDict{String, AbstractArray}("counts" => DS_counts[:, 2:end] |> Matrix)

    # DS_rowdata = CSV.File(joinpath(@__DIR__, "assets/DS_rowdata.csv")) |> DataFrame
    # DS_rowdata[!, :name] = ["strain$i" for i in 1:size(DS_rowdata, 1)]
    # select!(DS_rowdata, vcat("name", names(DS_rowdata)[1:end - 1]))

    # DS_coldata = CSV.File(joinpath(@__DIR__, "assets/DS_coldata.csv")) |> DataFrame
    # DS_coldata[!, :name] = DS_coldata[!, :sample]
    # select!(DS_coldata, vcat("name", names(DS_coldata)[1: end - 1]))

    # SummarizedExperiment(DS_assays, DS_rowdata, DS_coldata)

    import_se_from_csv(joinpath(@__DIR__, "assets/DS_assays.csv"),
                    joinpath(@__DIR__, "assets/DS_rowdata.csv"),
                    joinpath(@__DIR__, "assets/DS_coldata.csv"))

end


text:
using OceanStateEstimation, MeshArrays, Statistics
using Test

p=dirname(pathof(OceanStateEstimation))

@testset "OceanStateEstimation.jl" begin
    γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
    tmp=OceanStateEstimation.get_ecco_files(γ,"oceQnet")
    tmp=[mean(tmp[j][findall((!isnan).(tmp[j]))]) for j=1:5]
    ref=[19.88214831145215,47.63055475475805,-44.1122401210416,
         3.4402271721659816,30.14270126344508]
    @test tmp==ref

    lst=joinpath(p,"../examples/nctiles_climatology.csv")
    lists=dataverse_lists(lst)

    get_occa_velocity_if_needed()
    get_occa_variable_if_needed("DDuvel")
    @test isfile(joinpath(OCCAclim_path,"DDuvel.0406clim.nc"))

    get_ecco_velocity_if_needed()
    get_ecco_variable_if_needed("UVELMASS")
    @test isdir(joinpath(ECCOclim_path,"UVELMASS"))

    isdir(MITPROFclim_path)
    isdir(CBIOMESclim_path)
    @test true

end


text:
#############################################################################
# ProxAL/ExaTron Example File
# This example runs ProxAL with ExaTron as a backend and outputs a profile
# file. PProf has to be installed in the global environment.
#############################################################################

using ProxAL
using DelimitedFiles, Printf
using LinearAlgebra, JuMP, Ipopt
using CatViews
using CUDA
using MPI
using Logging

MPI.Init()

# Select one of the following
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case_ACTIVSg2000_Corrected", 20, 0.3, 1.0, 0.1, 1e5) # 290 iterations
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case_ACTIVSg2000_Corrected", 20, 0.3, 1.0, 0.01, 1e3) # 449 iterations
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case_ACTIVSg2000_Corrected", 20, 0.3, 1.0, 0.01, 1e5) # 521 iterations
(case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case_ACTIVSg2000_Corrected", 20, 0.3, 1.0, 0.1, 1e3) # 225 iterations
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case9241pegase", 2, 0.3, 0.8, 0.01, 1e3)
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case118", 168, 0.2, 1.0, 0.1, 1e5)
# (case, T, ramp_scale, load_scale, maxρ, quad_penalty) = ("case9", 6, 0.04, 1.0, 0.1, 0.1)

T = parse(Int, ARGS[1])
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
  println("$T periods")
end

# No contingencies in this example
K = 0

# Load case
DATA_DIR = joinpath(dirname(@__FILE__), "..", "data")
case_file = joinpath(DATA_DIR, "$(case).m")
# load_file = joinpath(DATA_DIR, "mp_demand", "$(case)_oneweek_168")
load_file = joinpath(DATA_DIR, "mp_demand", "$(case)_3600")

# Model/formulation settings
modelinfo = ModelInfo()
modelinfo.num_time_periods = T
modelinfo.load_scale = load_scale
modelinfo.ramp_scale = ramp_scale
modelinfo.allow_obj_gencost = true
modelinfo.allow_constr_infeas = false
modelinfo.weight_freq_ctrl = quad_penalty
modelinfo.time_link_constr_type = :penalty
modelinfo.ctgs_link_constr_type = :frequency_ctrl
modelinfo.case_name = case
modelinfo.num_ctgs = K

# Algorithm settings
algparams = AlgParams()
algparams.verbose = 1
algparams.tol = 1e-3
algparams.decompCtgs = true
algparams.iterlim = 100
algparams.device = ProxAL.CUDADevice
algparams.optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
algparams.tron_rho_pq=5*1e4
algparams.tron_rho_pa=5*1e5
algparams.tron_outer_iterlim=30
algparams.tron_inner_iterlim=2000
algparams.tron_scale=1e-5
algparams.mode = :coldstart
algparams.init_opf = false

# case_ACTIVSg2000_Corrected parameters
algparams.tron_rho_pq=5*1e4
algparams.tron_rho_pa=5*1e5
algparams.tron_outer_iterlim=30
algparams.tron_inner_iterlim=2000
algparams.tron_scale=1e-5

opt_sol, lyapunov_sol = Dict(), Dict()

ranks = MPI.Comm_size(MPI.COMM_WORLD)
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
   println("ProxAL/ExaTron $ranks ranks, $T periods")
end
cur_logger = global_logger(NullLogger())
elapsed_t = @elapsed begin
  redirect_stdout(devnull) do
    global nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, ProxAL.ExaTronBackend(), opt_sol, lyapunov_sol)
  end
end
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    global_logger(cur_logger)
    println("Creating problem: $elapsed_t")
    println("Benchmark Start")
    np = MPI.Comm_size(MPI.COMM_WORLD)
    elapsed_t = @elapsed begin
      info = ProxAL.optimize!(nlp)
    end
    println("AugLag iterations: $(info.iter) with $np ranks in $elapsed_t seconds")
else
  info = ProxAL.optimize!(nlp)
end

MPI.Finalize()


text:
module LibCImGui

using CImGui_jll
export CImGui_jll

using CSyntax.CEnum

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libcimgui_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libcimgui_api.jl"))
include(joinpath(@__DIR__, "..", "gen", "libhelper_api.jl"))

foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "ig") || startswith(string(s), "Im")
        @eval export $s
    end
end

end # module


text:
using ArrayInterface
using Pkg

function dev_subpkg(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.develop(PackageSpec(path=subpkg_path))
end
dev_subpkg("ArrayInterfaceCore")

using ArrayInterfaceCore
using Documenter

makedocs(;
    modules=[ArrayInterface, ArrayInterfaceCore],
    sitename="ArrayInterface",
    pages=[
        "Home" => "index.md",
        "API" => "api.md"
    ]
)

deploydocs(;
    repo="github.com/JuliaArrays/ArrayInterface.jl"
)


text:
module SparseRegressionTest
using SparseRegression, Test
include("datagenerator.jl")

#-----------------------------------------------------------------------# Sanity Check and show
@testset "Sanity Check" begin
    n, p = 1000, 5
    x, y, β = SparseRegression.fakedata(L2DistLoss(), n, p)
    show(SModel(x, y))
    println("\n")
end

name(o) = replace(string(typeof(o)), r"([a-zA-Z]*\.)", "")

function _test(l::Loss, p::Penalty, alg)
    x, y, β = SparseRegression.fakedata(l, 100, 5)
    o = @inferred SModel(x, y, l, p)
    a = @inferred alg(o)
    @inferred learn!(o, strategy(a, MaxIter(3)))
    @test coef(o) == o.β
    @test length(predict(o, x)) == length(y)

    w = rand(100)
    o = @inferred SModel(x, y, l, Weights(w), p)
    a = @inferred alg(o)
    learn!(o, strategy(a, MaxIter(3)))
    coef(o)
    predict(o, x)
end

@testset "Sanity Check" begin
    for a in [ProxGrad, Fista, AdaptiveProxGrad, GradientDescent]
        @info(a)
        for l in [L2DistLoss(), .5*L2DistLoss(), HuberLoss(2.0), QuantileLoss(.7),
            LogitMarginLoss(), DWDMarginLoss(2.0), L2HingeLoss(), ExpLoss()]
            println("  > $l")
            for p in [NoPenalty(), L2Penalty(), L1Penalty(), ElasticNetPenalty()]
                _test(l, p, a)
            end
        end
        println()
    end
    for a in [Sweep, LinRegCholesky]
        for l in [L2DistLoss(), .5L2DistLoss()]
            for p in [NoPenalty(), L2Penalty()]
                _test(l ,p, a)
            end
        end
    end
end

@testset "Linear Regression" begin
    n, p = 1000, 10
    x, y, β = DataGenerator.linregdata(n, p)
    s = SModel(x, y, .5 * L2DistLoss(), NoPenalty())

    for a in [ProxGrad(s, .1), Fista(s, .1), GradientDescent(s, .1), Sweep(s), LinRegCholesky(s)]
        learn!(s, strategy(a, MaxIter(25)))
        @test coef(s) ≈ x\y atol = .5
    end
end

@testset "LineSearch" begin
    n, p = 1000, 10
    x, y, β = SparseRegression.fakedata(L2DistLoss(), n, p)
    s = SModel(x, y, L1Penalty())
    strat = strategy(MaxIter(5), LineSearch(GradientDescent(s)))
    learn!(s, strat)
    @test issorted(coef(s))
    strat = strategy(MaxIter(5), LineSearch(ProxGrad(s)))
    learn!(s, strat)
    @test issorted(coef(s))
    strat = strategy(MaxIter(5), LineSearch(Fista(s)))
    learn!(s, strat)
    @test issorted(coef(s))
end

@testset "SModel" begin
    n, p = 100, 5
    x, y, β = DataGenerator.linregdata(n, p)
    @testset "Constructor type stability" begin
        @inferred SModel(x, y)

        @inferred SModel(x, y, L2DistLoss())
        @inferred SModel(x, y, L2Penalty())
        @inferred SModel(x, y, rand(5))

        @inferred SModel(x, y, L2DistLoss(), L2Penalty())
        @inferred SModel(x, y, L2DistLoss(), rand(5))
        @inferred SModel(x, y, L2Penalty(), L2DistLoss())
        @inferred SModel(x, y, L2Penalty(), rand(5))
        @inferred SModel(x, y, rand(5), L2DistLoss())
        @inferred SModel(x, y, rand(5), L2Penalty())

        @inferred SModel(x, y, L2DistLoss(), L2Penalty(), rand(5))
        @inferred SModel(x, y, L2DistLoss(), rand(5), L2Penalty())
        @inferred SModel(x, y, L2Penalty(), L2DistLoss(), rand(5))
        @inferred SModel(x, y, L2Penalty(), rand(5), L2DistLoss())
        @inferred SModel(x, y, rand(5), L2DistLoss(), L2Penalty())
        @inferred SModel(x, y, rand(5), L2Penalty(), L2DistLoss())
    end
    @testset "predict" begin
        o = SModel(x, y)
        @test predict(o, randn(5)) == 0.0
        @test predict(o, randn(10, 5)) == zeros(10)
    end
end

end


text:
#pansharpen?
#others?


text:
const VPtr = Ptr{Cvoid}

"""
    get_string(x::VPtr, fn_sym)::Maybe{String}

C-call the SBML function `fn_sym` with a single parameter `x`, interpret the
result as a string and return it, or throw exception in case the pointer is
NULL.
"""
function get_string(x::VPtr, fn_sym)::String
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        throw(DomainError(x, "Calling $fn_sym returned NULL, valid string expected."))
    end
end

"""
    get_optional_string(x::VPtr, fn_sym)::Maybe{String}

Like [`get_string`](@ref), but returns `nothing` instead of throwing an
exception.

This is used to get notes and annotations and several other things (see
`get_notes`, `get_annotations`)
"""
function get_optional_string(x::VPtr, fn_sym)::Maybe{String}
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        return nothing
    end
end

"""
    get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}

Helper for getting out boolean flags.
"""
function get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cint, (VPtr,), x) != 0
    else
        return nothing
    end
end

"""
    get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{UInt}

Helper for getting out unsigned integers.
"""
function get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{Int}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cint, (VPtr,), x)
    else
        return nothing
    end
end

"""
    get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}

Helper for getting out C doubles aka Float64s.
"""
function get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}
    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
        return ccall(sbml(get_sym), Cdouble, (VPtr,), x)
    else
        return nothing
    end
end

"""
    function readSBML(fn::String, sbml_conversion = model->nothing)::SBML.Model

Read the SBML from a XML file in `fn` and return the contained `SBML.Model`.

The `sbml_conversion` is a function that does an in-place modification of the
single parameter, which is the C pointer to the loaded SBML document (C type
`SBMLDocument*`). Several functions for doing that are prepared, including
[`set_level_and_version`](@ref), [`libsbml_convert`](@ref), and
[`convert_simplify_math`](@ref).

# Example
```
m = readSBML("my_model.xml", doc -> begin
    set_level_and_version(3, 1)(doc)
    convert_simplify_math(doc)
end)
```
"""
function readSBML(fn::String, sbml_conversion = document -> nothing)::SBML.Model
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
            msg = strip(get_string(err, :XMLError_getMessage))
            @error "SBML reported error: $msg"
        end
        if n_errs > 0
            throw(AssertionError("Opening SBML document has reported errors"))
        end

        sbml_conversion(doc)

        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
            throw(AssertionError("SBML document contains no model"))
        end

        model = ccall(sbml(:SBMLDocument_getModel), VPtr, (VPtr,), doc)

        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
    end
end

get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)

"""
    function getAssociation(x::VPtr)::GeneProductAssociation

Convert a pointer to SBML `FbcAssociation_t` to the `GeneProductAssociation`
tree structure.
"""
function getAssociation(x::VPtr)::GeneProductAssociation
    # libsbml C API is currently missing functions to check this in a normal
    # way, so we use a bit of a hack.
    typecode = ccall(sbml(:SBase_getTypeCode), Cint, (VPtr,), x)
    if typecode == 808 # SBML_FBC_GENEPRODUCTREF
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
    elseif typecode == 809 # SBML_FBC_AND
        return GPAAnd([
            getAssociation(
                ccall(sbml(:FbcAnd_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcAnd_getNumAssociations), Cuint, (VPtr,), x)
        ])
    elseif typecode == 810 # SBML_FBC_OR
        return GPAOr([
            getAssociation(
                ccall(sbml(:FbcOr_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcOr_getNumAssociations), Cuint, (VPtr,), x)
        ])
    else
        throw(ErrorException("Unsupported FbcAssociation type"))
    end
end

""""
    function extractModel(mdl::VPtr)::SBML.Model

Take the `SBMLModel_t` pointer and extract all information required to make a
valid [`SBML.Model`](@ref) structure.
"""
function extractModel(mdl::VPtr)::SBML.Model
    # get the FBC plugin pointer (FbcModelPlugin_t)
    mdl_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), mdl, "fbc")

    # get the parameters
    parameters = Dict{String,Float64}()
    for i = 1:ccall(sbml(:Model_getNumParameters), Cuint, (VPtr,), mdl)
        p = ccall(sbml(:Model_getParameter), VPtr, (VPtr, Cuint), mdl, i - 1)
        id = get_string(p, :Parameter_getId)
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

    # parse out the unit definitions
    units = Dict{String,Vector{SBML.UnitPart}}()
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
        id = get_string(ud, :UnitDefinition_getId)
        units[id] = [
            begin
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
                SBML.UnitPart(
                    unsafe_string(
                        ccall(
                            sbml(:UnitKind_toString),
                            Cstring,
                            (Cint,),
                            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
                        ),
                    ),
                    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getScale), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u),
                )
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
        ]
    end

    # parse out compartment names
    compartments = Dict{String,Compartment}()
    for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
        co = ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1)

        compartments[get_string(co, :Compartment_getId)] = Compartment(
            get_optional_string(co, :Compartment_getName),
            get_optional_bool(co, :Compartment_isSetConstant, :Compartment_getConstant),
            get_optional_int(
                co,
                :Compartment_isSetSpatialDimensions,
                :Compartment_getSpatialDimensions,
            ),
            get_optional_double(co, :Compartment_isSetSize, :Compartment_getSize),
            get_optional_string(co, :Compartment_getUnits),
            get_notes(co),
            get_annotation(co),
        )
    end

    # parse out species
    species = Dict{String,Species}()
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t

        formula = nothing
        charge = nothing
        if sp_fbc != C_NULL
            # if the FBC plugin is present, try to get the chemical formula and charge
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
                formula = get_string(sp_fbc, :FbcSpeciesPlugin_getChemicalFormula)
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
        end

        ia = nothing
        if ccall(sbml(:Species_isSetInitialAmount), Cint, (VPtr,), sp) != 0
            ia = (
                ccall(sbml(:Species_getInitialAmount), Cdouble, (VPtr,), sp),
                get_optional_string(sp, :Species_getSubstanceUnits),
            )
        end

        ic = nothing
        if ccall(sbml(:Species_isSetInitialConcentration), Cint, (VPtr,), sp) != 0
            ic = (
                ccall(sbml(:Species_getInitialConcentration), Cdouble, (VPtr,), sp),
                get_optional_string(sp, :Species_getSubstanceUnits),
            )
        end

        species[get_string(sp, :Species_getId)] = Species(
            get_optional_string(sp, :Species_getName),
            get_string(sp, :Species_getCompartment),
            get_optional_bool(
                sp,
                :Species_isSetBoundaryCondition,
                :Species_getBoundaryCondition,
            ),
            formula,
            charge,
            ia,
            ic,
            get_optional_bool(
                sp,
                :Species_isSetHasOnlySubstanceUnits,
                :Species_getHasOnlySubstanceUnits,
            ),
            get_notes(sp),
            get_annotation(sp),
        )
    end

    # parse out the flux objectives (these are complementary to the objectives
    # that appear in the reactions, see comments lower)
    objectives_fbc = Dict{String,Float64}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumObjectives), Cuint, (VPtr,), mdl_fbc)
            o = ccall(
                sbml(:FbcModelPlugin_getObjective),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )
            for j = 1:ccall(sbml(:Objective_getNumFluxObjectives), Cuint, (VPtr,), o)
                fo = ccall(sbml(:Objective_getFluxObjective), VPtr, (VPtr, Cuint), o, j - 1)
                objectives_fbc[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
            end
        end
    end

    # reactions!
    reactions = Dict{String,Reaction}()
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
        lb = (-Inf, "") # (bound value, unit id)
        ub = (Inf, "")
        oc = 0.0
        math = nothing

        # kinetic laws store a second version of the bounds and objectives
        kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
        if kl != C_NULL
            for j = 1:ccall(sbml(:KineticLaw_getNumParameters), Cuint, (VPtr,), kl)
                p = ccall(sbml(:KineticLaw_getParameter), VPtr, (VPtr, Cuint), kl, j - 1)
                id = get_string(p, :Parameter_getId)
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
                punit = () -> get_string(p, :Parameter_getUnits)
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
            end

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
            end
        end

        # TRICKY: SBML spec is completely silent about what should happen if
        # someone specifies both the above and below formats of the flux bounds
        # for one reaction. Notably, these do not really specify much
        # interaction with units. In this case, we'll just set a special
        # "[fbc]" unit that has no specification in `units`, and hope the users
        # can make something out of it.
        re_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), re, "fbc")
        if re_fbc != C_NULL
            fbcb = get_optional_string(re_fbc, :FbcReactionPlugin_getLowerFluxBound)
            if !isnothing(fbcb) && haskey(parameters, fbcb)
                lb = (parameters[fbcb], "[fbc]")
            end
            fbcb = get_optional_string(re_fbc, :FbcReactionPlugin_getUpperFluxBound)
            if !isnothing(fbcb) && haskey(parameters, fbcb)
                ub = (parameters[fbcb], "[fbc]")
            end
        end

        # extract stoichiometry
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
                stoi[get_string(sr, :SpeciesReference_getSpecies)] =
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor

        # reactants and products
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
            add_stoi(sr, -1)
        end

        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
            add_stoi(sr, 1)
        end

        # gene product associations
        association = nothing
        if re_fbc != C_NULL
            gpa = ccall(
                sbml(:FbcReactionPlugin_getGeneProductAssociation),
                VPtr,
                (VPtr,),
                re_fbc,
            )
            if gpa != C_NULL
                a = ccall(sbml(:GeneProductAssociation_getAssociation), VPtr, (VPtr,), gpa)
                a != C_NULL
                association = getAssociation(a)
            end
        end

        # explicit reversible flag (defaults to true in SBML)
        reversible = Bool(ccall(sbml(:Reaction_getReversible), Cint, (VPtr,), re))

        reid = get_string(re, :Reaction_getId)
        reactions[reid] = Reaction(
            stoi,
            lb,
            ub,
            haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
            association,
            math,
            reversible,
            get_notes(re),
            get_annotation(re),
        )
    end

    # extract gene products
    gene_products = Dict{String,GeneProduct}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumGeneProducts), Cuint, (VPtr,), mdl_fbc)
            gp = ccall(
                sbml(:FbcModelPlugin_getGeneProduct),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )

            id = get_optional_string(gp, :GeneProduct_getId) # IDs don't need to be set

            if id != nothing
                gene_products[id] = GeneProduct(
                    get_optional_string(gp, :GeneProduct_getName),
                    get_optional_string(gp, :GeneProduct_getLabel),
                    get_notes(gp),
                    get_annotation(gp),
                )
            end
        end
    end

    function_definitions = Dict{String,FunctionDefinition}()
    for i = 1:ccall(sbml(:Model_getNumFunctionDefinitions), Cuint, (VPtr,), mdl)
        fd = ccall(sbml(:Model_getFunctionDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
        def = nothing
        if ccall(sbml(:FunctionDefinition_isSetMath), Cint, (VPtr,), fd) != 0
            def = parse_math(ccall(sbml(:FunctionDefinition_getMath), VPtr, (VPtr,), fd))
        end

        function_definitions[get_string(fd, :FunctionDefinition_getId)] =
            FunctionDefinition(
                get_optional_string(fd, :FunctionDefinition_getName),
                def,
                get_notes(fd),
                get_annotation(fd),
            )
    end

    return Model(
        parameters,
        units,
        compartments,
        species,
        reactions,
        gene_products,
        function_definitions,
        get_notes(mdl),
        get_annotation(mdl),
    )
end


text:
function measData(f::Union{MDFFileV2, MDFv2InMemory}, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),
  receivers=1:rxNumChannels(f))

  _data = measDataRaw(f)
  if measIsFastFrameAxis(f)
    data = _data[frames, :, receivers, periods]
    data = reshape(data, length(frames), size(_data,2), length(receivers), length(periods))
  else
    data = _data[:, receivers, periods, frames]
    data = reshape(data, size(_data,1), length(receivers), length(periods), length(frames))
  end
  return data
end

function measDataTDPeriods(f::Union{MDFFileV2, MDFv2InMemory}, periods=1:acqNumPeriods(f),
  receivers=1:rxNumChannels(f))
  if measIsFastFrameAxis(f)
    error("measDataTDPeriods can currently not handle transposed data!")
  end

  _data = measDataRaw(f)
  data = reshape(_data,Val(3))[:, receivers, periods]

  return data
end

function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=true)
  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) ||
    !measIsFourierTransformed(f)
    return nothing
  end

  rows_ = rowsToSubsampledRows(f, rows)

  data_ = measDataRaw(f)

  data_ = reshape(data_, size(data_,1),
          size(data_,2)*size(data_,3),
          size(data_,4))[:, rows_, :]
  data = reshape(data_, Val(2))

  fgdata = data[measFGFrameIdx(f),:]

  if measIsSparsityTransformed(f)
    dataBackTrafo = similar(fgdata, prod(calibSize(f)), size(fgdata,2))
    B = linearOperator(measSparsityTransformation(f), calibSize(f))

    tmp = measSubsamplingIndices(f)
    subsamplingIndices_ = reshape(tmp, size(tmp,1),
                          size(tmp,2)*size(tmp,3),
                          size(tmp,4))[:, rows_, :]
    subsamplingIndices = reshape(subsamplingIndices_, Val(2))

    for l=1:size(fgdata,2)
      dataBackTrafo[:,l] .= 0.0
      dataBackTrafo[subsamplingIndices[:,l],l] .= fgdata[:,l]
      dataBackTrafo[:,l] .= adjoint(B) * vec(dataBackTrafo[:,l])
    end
    fgdata = dataBackTrafo
  end

  if bgCorrection # this assumes equidistant bg frames
    @debug "Applying bg correction on system matrix (MDF)"
    bgdata = data[measBGFrameIdx(f),:]
    blockLen = measBGFrameBlockLengths( invpermute!(deepcopy(measIsBGFrame(f)), measFramePermutation(f)) ) # Added deepcopy to be side-effect free in in-memory MDF
    st = 1
    for j=1:length(blockLen)
      bgdata[st:st+blockLen[j]-1,:] .=
      mean(bgdata[st:st+blockLen[j]-1,:], dims=1)
      st += blockLen[j]
    end

    bgdataInterp = interpolate(bgdata, (BSpline(Linear()), NoInterp()))
    # Cubic does not work for complex numbers
    origIndex = measFramePermutation(f)
    M = size(fgdata,1)
    K = size(bgdata,1)
    N = M + K
    for m=1:M
      alpha = (origIndex[m]-1)/(N-1)*(K-1)+1
      for k=1:size(fgdata,2)
        fgdata[m,k] -= bgdataInterp(alpha,k)
      end
    end
  end
  return fgdata
end

function systemMatrixWithBG(f::Union{MDFFileV2, MDFv2InMemory})
  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) || !measIsFourierTransformed(f)
    return nothing
  end

  data_ = measDataRaw(f)
  data = data_[:, :, :, :]
  return data
end

# This is a special variant used for matrix compression
function systemMatrixWithBG(f::Union{MDFFileV2, MDFv2InMemory}, freq)
  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) || !measIsFourierTransformed(f)
    return nothing
  end

  data_ = measDataRaw(f)
  data = data_[:, freq, :, :]
  return data
end


text:
using LinearAlgebra
using Fermi.DIIS
# Functions specific to CTF
#
############ PREPROCESSING ################
"""
    Fermi.CoupledCluster.RCCSD{T}(Alg::CTF)

Compute a RCCSD wave function using the Compiled time factorization algorithm (CTF)
"""
function RCCSD{T}(guess::RCCSD{Tb},Alg::CTF) where { T <: AbstractFloat,
                                                    Tb <: AbstractFloat }
    refwfn = Fermi.HartreeFock.RHF()
    RCCSD{T}(refwfn, guess, Alg)
end

function RCCSD{T}(refwfn::RHF, guess::RCCSD{Tb},Alg::CTF) where { T <: AbstractFloat,
                                                    Tb <: AbstractFloat }
    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    @output "Transforming Integrals..."
    tint = @elapsed begin
        refwfn.ints["OOOO"]
        refwfn.ints["OOOV"]
        refwfn.ints["OOVV"]
        refwfn.ints["OVOV"]
        refwfn.ints["OVVV"]
        refwfn.ints["VVVV"]
    end
    @output " done in {} s" tint
    RCCSD{T}(refwfn, guess, Alg)
end

function RCCSD{T}(refwfn::RHF, guess::RCCSD{Tb}, Alg::CTF) where { T <: AbstractFloat,
                                                                  Tb <: AbstractFloat }
    ints = refwfn.ints
    d = [i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
    D = [i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]
    newT1 = ints["FOV"]./d
    newT2 = ints["OOVV"]./D
    o_small,v_small = size(guess.T1)
    o = 1:o_small
    v = 1:v_small
    newT1[o,v] .= Fermi.data(guess.T1)
    newT2[o,o,v,v] .= Fermi.data(guess.T2)
    RCCSD{T}(refwfn, ints, newT1, newT2, Alg)
end
############ END PREPROCESSING ################
#
#
#
############# KERNEL FUNCTIONS ###################
function print_alg(alg::CTF)
    @output "Computing RCCSD with CTF algorithm.\n"
end

function compute_integrals(ints,alg::CTF)
    ints["OOOO"]
    ints["OOOV"]
    ints["OOVV"]
    ints["OVOV"]
    ints["OVVV"]
    ints["VVVV"]
end
function update_energy(T1::Array{T, 2}, T2::Array{T, 4}, f::Array{T,2}, ints, alg::CTF) where { T <: AbstractFloat }

    Voovv = ints["OOVV"]
    @tensoropt (k=>x, l=>x, c=>100x, d=>100x)  begin
        CC_energy = 2.0*f[k,c]*T1[k,c]
        B[l,c,k,d] := -1.0*T1[l,c]*T1[k,d]
        B[l,c,k,d] += -1.0*T2[l,k,c,d]
        B[l,c,k,d] += 2.0*T2[k,l,c,d]
        CC_energy += B[l,c,k,d]*Voovv[k,l,c,d]
        CC_energy += 2.0*T1[l,c]*T1[k,d]*Voovv[l,k,c,d]
    end

    return CC_energy
end

function compute_oovv(ints,alg::CTF)
    ints["OOVV"]
end

function update_T1(T1::Array{T,2}, T2::Array{T,4}, newT1::Array{T,2}, foo, fov, fvv, ints::IntegralHelper, alg::CTF) where T <: AbstractFloat
    fill!(newT1, 0.0)
    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = ints["OOOO"], ints["OOOV"], ints["OOVV"], ints["OVOV"], ints["OVVV"], ints["VVVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT1[i,a] += fov[i,a]
        newT1[i,a] -= foo[i,k]*T1[k,a]
        newT1[i,a] += fvv[c,a]*T1[i,c]
        newT1[i,a] -= fov[k,c]*T1[i,c]*T1[k,a]
        newT1[i,a] += 2.0*fov[k,c]*T2[i,k,a,c]
        newT1[i,a] -= fov[k,c]*T2[k,i,a,c]
        newT1[i,a] -= T1[k,c]*Vovov[i,c,k,a]
        newT1[i,a] += 2.0*T1[k,c]*Voovv[k,i,c,a]
        newT1[i,a] -= T2[k,i,c,d]*Vovvv[k,a,d,c]
        newT1[i,a] += 2.0*T2[i,k,c,d]*Vovvv[k,a,d,c]
        newT1[i,a] += -2.0*T2[k,l,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += T2[l,k,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += -2.0*T1[k,c]*T1[l,a]*Vooov[l,k,i,c]
        newT1[i,a] -= T1[k,c]*T1[i,d]*Vovvv[k,a,d,c]
        newT1[i,a] += 2.0*T1[k,c]*T1[i,d]*Vovvv[k,a,c,d]
        newT1[i,a] += T1[k,c]*T1[l,a]*Vooov[k,l,i,c]
        newT1[i,a] += -2.0*T1[k,c]*T2[i,l,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[k,c]*T2[l,i,a,d]*Voovv[k,l,c,d]
        newT1[i,a] += T1[k,c]*T2[l,i,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[i,c]*T2[l,k,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += T1[i,c]*T2[l,k,a,d]*Voovv[k,l,c,d]
        newT1[i,a] += -2.0*T1[l,a]*T2[i,k,d,c]*Voovv[k,l,c,d]
        newT1[i,a] += T1[l,a]*T2[i,k,c,d]*Voovv[k,l,c,d]
        newT1[i,a] += T1[k,c]*T1[i,d]*T1[l,a]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[k,c]*T1[i,d]*T1[l,a]*Voovv[k,l,c,d]
        newT1[i,a] += 4.0*T1[k,c]*T2[i,l,a,d]*Voovv[k,l,c,d]
    end
end

function update_T2(T1::Array{T,2},T2::Array{T,4},newT2::Array{T,4},foo,fov,fvv,ints::IntegralHelper, alg::CTF) where T <: AbstractFloat
    fill!(newT2, 0.0)
    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = ints["OOOO"], ints["OOOV"], ints["OOVV"], ints["OVOV"], ints["OVVV"], ints["VVVV"]

    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT2[i,j,a,b] += Voovv[i,j,a,b]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*Vvvvv[c,d,a,b]
        newT2[i,j,a,b] += T2[i,j,c,d]*Vvvvv[c,d,a,b]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] += T2[k,l,a,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,a]*Vovvv[k,b,c,d]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,b]*Vovvv[k,a,d,c]
        newT2[i,j,a,b] += T1[i,c]*T1[k,a]*T1[l,b]*Vooov[l,k,j,c]
        newT2[i,j,a,b] += T1[j,c]*T1[k,a]*T1[l,b]*Vooov[k,l,i,c]
        newT2[i,j,a,b] += T2[k,l,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[l,k,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[i,k,a,c]*T2[l,j,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[k,i,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,a,c]*T2[l,k,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,j,a,c]*T2[k,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,j,a,c]*T2[i,l,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += 4.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,d,c]*T2[l,k,a,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T1[k,a]*T1[l,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T2[l,k,a,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*T2[i,j,d,c]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] := -1.0*foo[i,k]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += fvv[c,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[k,b]*Vooov[j,i,k,a]
        P_OoVv[i,j,a,b] += T1[j,c]*Vovvv[i,c,a,b]
        P_OoVv[i,j,a,b] += -1.0*fov[k,c]*T1[i,c]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += -1.0*fov[k,c]*T1[k,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,i,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[i,c]*T1[k,a]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[i,c]*T1[k,b]*Vovov[j,c,k,a]
        P_OoVv[i,j,a,b] += 2.0*T2[i,k,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[i,k,a,c]*Vovov[j,c,k,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,j,a,c]*Vovov[i,c,k,b]
        P_OoVv[i,j,a,b] += -2.0*T1[l,b]*T2[i,k,a,c]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[k,i,a,c]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[i,k,d,b]*Vovvv[k,a,c,d]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[k,i,a,d]*Vovvv[k,b,d,c]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[i,k,a,d]*Vovvv[k,b,c,d]
        P_OoVv[i,j,a,b] += T1[j,c]*T2[l,k,a,b]*Vooov[l,k,i,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[i,k,a,c]*Vooov[k,l,j,c]
        P_OoVv[i,j,a,b] += -1.0*T1[k,a]*T2[i,j,d,c]*Vovvv[k,b,d,c]
        P_OoVv[i,j,a,b] += T1[k,a]*T2[i,l,c,b]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += 2.0*T1[j,c]*T2[i,k,a,d]*Vovvv[k,b,d,c]
        P_OoVv[i,j,a,b] += -1.0*T1[k,c]*T2[i,j,a,d]*Vovvv[k,b,d,c]
        P_OoVv[i,j,a,b] += 2.0*T1[k,c]*T2[i,j,a,d]*Vovvv[k,b,c,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T2[i,l,a,b]*Vooov[k,l,j,c]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T2[i,l,a,b]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += T2[j,k,c,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T1[j,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[j,d]*T2[i,l,a,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T1[l,a]*T2[i,j,d,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[l,a]*T2[i,j,d,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,b,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[i,c]*T1[k,a]*T2[j,l,b,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,d,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[l,b]*T2[k,j,a,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T2[i,k,d,c]*T2[l,j,a,b]*Voovv[k,l,c,d]

        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
end
#
#
############# END KERNEL FUNCTIONS ###################


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

"""
    HybridBranching(depth_cutoff::Int,
                    shallow_rule::VariableBranchingRule,
                    deep_rule::::VariableBranchingRule)

Branching strategy that switches between two variable branching strategies,
according to the depth of the node.

More specifically, if `node.depth <= depth_cutoff`, then `shallow_rule` is
applied. Otherwise, `deep_rule` is applied.
"""
mutable struct HybridBranching <: VariableBranchingRule
    depth_cutoff::Int
    shallow_rule::VariableBranchingRule
    deep_rule::VariableBranchingRule
end

HybridBranching(depth_cutoff::Int = 10) =
    HybridBranching(depth_cutoff, StrongBranching(), PseudocostBranching())

function find_branching_var(rule::HybridBranching, node::Node, pool::NodePool)::Variable
    if node.depth <= rule.depth_cutoff
        return find_branching_var(rule.shallow_rule, node, pool)
    else
        return find_branching_var(rule.deep_rule, node, pool)
    end
end


text:
# Here are the main exported functions and types. Note there are other exported functions and types in files such as "effective_waves/effective_waves_export" and "average_waves/average_waves_export.jl."

export  EffectiveWave, AverageWave # the two main types
export  MatchWave # a combination of the two types above

# for MatchWave
export  match_error, x_mesh_match

# for discrete method
export  x_mesh

# for material and particle properties
export  Specie, Medium, volume_fraction, Zn, t_vectors, Nn, p_speed, maximum_hankel_order

# for effective waves
export  wavenumbers, wavenumber, effective_waves, transmission_angle,
        effective_wavevectors, scattering_amplitudes_average, scale_amplitudes_effective

export  reflection_coefficient, reflection_coefficients

export  effective_medium, reflection_coefficient_halfspace

# List of shorthand for some materials
export  Brick, IronArmco, LeadAnnealed, RubberGum, FusedSilica, GlassPyrex,
        ClayRock, WaterDistilled, Glycerol, Hexadecane, Acetone, Benzene,
        Nitrobenzene, OliveOil, SodiumNitrate, AirDry,
        LimeStone, Clay, Calcite, SilicaQuartz

import Base.isequal, Base.(==), Base.zero
import SpecialFunctions: besselj, hankelh1
import Statistics: mean, std