text:
module OpticalBlochEquations

using StaticArrays
using StructArrays
using Parameters
using HalfIntegers
using Unitful
using LinearAlgebra
using WignerSymbols
using LoopVectorization
using PhysicalConstants.CODATA2018

const h = PlanckConstant.val
const ħ = h / 2π
const c = SpeedOfLightInVacuum.val
export h, ħ, c

macro with_unit(arg1, arg2)
    arg2 = @eval @u_str $arg2
    return convert(Float64, upreferred(eval(arg1) .* arg2).val)
end
export @with_unit

macro params(fields_tuple)
    fields = fields_tuple.args
    esc(
        quote
            NamedTuple{Tuple($fields)}(($fields_tuple))
        end)
end
export @params

const sph2cart = @SMatrix [
    -1/√2 +1/√2 0;
    +im/√2 +im/√2 0;
    0 0 1;
]
export sph2cart

const ϵ₊ = @SVector [-1/√2, -im/√2, 0]
const ϵ₋ = @SVector [ 1/√2, -im/√2, 0]
const ϵ₀ = @SVector [    0,      0, 1]
const ϵ = [ϵ₋, ϵ₀, ϵ₊]
export ϵ₊, ϵ₋, ϵ₀
export ϵ

const qs = @SVector [-1, 0, 1]
export qs

x̂ = @SVector [1,0,0]
ŷ = @SVector [0,1,0]
ẑ = @SVector [0,0,1]
#export x̂, ŷ, ẑ

# @with_kw struct Hamiltonian12
#     states::Vector{State}
#     H::StructArray{<:Complex} = StructArray(zeros(ComplexF64, length(states), length(states)))
#     Hsubs::Vector{StructArray{<:Complex}} = []
# end
# export Hamiltonian12

# function make_hamiltonian(states::Vector(State), lasers::Vector(Laser))
#     H = Hamiltonian(; states)
#     return H
# end

@with_kw struct Laser
    k::SVector{3, Float64}         # k-vector
    ϵ_re::SVector{3, Float64}      # real part of polarization
    ϵ_im::SVector{3, Float64}      # imaginary part of polarization
    ω::Float64                     # frequency
    s::Float64                     # saturation parameter
    kr::Float64                    # value of k ⋅ r, defaults to 0
    f_re_q::SVector{3, Float64}    # real part of f = exp(i(kr - ωt))
    f_im_q::SVector{3, Float64}    # imag part of f = exp(i(kr - ωt))
    re::Float64
    im::Float64
    Laser(k, ϵ, ω, s) = new(k, real.(ϵ), imag.(ϵ), ω, s, 0.0, SVector(0, 0, 0), SVector(0, 0, 0), 1.0, 0.0)
end
export Laser

@with_kw struct Field
    k::SVector{3, Float64}      # k-vector
    e::SVector{3, ComplexF64}   # polarization
    ω::Float64                  # frequency
    s::Float64                  # saturation parameter
    fre::SVector{3, Float64}
    fim::SVector{3, Float64}
    Hqm::Array{Float64, 2}
    Hq0::Array{Float64, 2}
    Hqp::Array{Float64, 2}
end

@with_kw struct State
    F::HalfInt                          # angular quantum number
    m::HalfInt                          # projection of angular quantum number
    ω::Float64                          # frequency
    μ::Float64                          # magnetic moment
    Γ::Union{Nothing, Float64}=nothing  # linewidth (defaults to `nothing` for ground states)
end
# export State

function define_field(k, e, ω, s)
    fre = zeros(Float64, 3)
    fim = zeros(Float64, 3)
    for q in eachindex(qs)
        dotted = e ⋅ ϵ[q]
        fre[q] = real(dotted)
        fim[q] = imag(dotted)
    end
    Hqm = zeros(Float64, (3, 3))
    Hq0 = zeros(Float64, (3, 3))
    Hqp = zeros(Float64, (3, 3))
    return Field(k, e, ω, s, fre, fim, Hqm, Hq0, Hqp)
end
export define_field

@with_kw struct Manifold
    F::HalfInt                          # angular quantum number
    states::Vector{State}
end
manifold(; F, ω, μ, Γ=nothing) = Manifold(F, [State(F=F, m=m, ω=ω, μ=μ, Γ=Γ) for m in -F:F])
export manifold

# Structure for quantum jumps from state `s` to state `s′` with rate `r`
struct Jump
    s ::Int64
    s′::Int64
    r ::Float64
end

# Round `val` to the nearest multiple of `prec`
round_to_mult(val, prec) = (inv_prec = 1 / prec; round.(val * inv_prec) / inv_prec)

function round_freq(ω, Γ, freq_res)
    ω_min = freq_res * Γ
    return round_to_mult(ω, ω_min)
end
export round_freq

function round_vel(v, λ, Γ, freq_res)
    v_min = freq_res * Γ * λ / 2π
    return round_to_mult(v, v_min)
end
export round_vel

function schrödinger(states, lasers, d, f, p)

    n_states = length(states)
    n_lasers = length(lasers)

    lasers = StructArray(lasers)

    r = @SVector [0.,0.,0.]
    v = @SVector [0.,0.,0.]

    type_complex = ComplexF64
    type_real    = Float64

    H = StructArray( zeros(type_complex, n_states, n_states) )

    # Define the optical Hamiltonian; it has dimensions (k, q)
    for i ∈ 1:n_states, j ∈ i:n_states
        for k in eachindex(lasers)
            lasers[k].Hqm[i,j] += d[1,i,j] * (lasers[k].e ⋅ ϵ[1])
            lasers[k].Hq0[i,j] += d[2,i,j] * (lasers[k].e ⋅ ϵ[2])
            lasers[k].Hqp[i,j] += d[3,i,j] * (lasers[k].e ⋅ ϵ[3])
            lasers[k].Hqm[j,i] += d[1,i,j] * (lasers[k].e ⋅ ϵ[1])
            lasers[k].Hq0[j,i] += d[2,i,j] * (lasers[k].e ⋅ ϵ[2])
            lasers[k].Hqp[j,i] += d[3,i,j] * (lasers[k].e ⋅ ϵ[3])
        end
    end

    ψ = zeros(type_complex, n_states)
    dψ = deepcopy(ψ)
    ψ[1] = 1
    ψ_soa = StructArray(ψ)
    dψ_soa = StructArray(dψ)

    p = @params (H, ψ_soa, dψ_soa, states, lasers, r, v, f, p)

    return (dψ, ψ, p)
end
export schrödinger

mutable struct Particle
    r0::SVector{3, Float64}
    r::SVector{3, Float64}
    v::SVector{3, Float64}
end

"""
    freq_res::Float: all frequencies are rounded by this value (in units of Γ)
"""
function obe(states, lasers, d, ρ; freq_res=1e-2)

    Γs = [s.Γ for s in states]
    Γ = maximum(filter(x -> x != nothing, Γs))

    λs = [c / s.ω for s in states]
    λs_nonzero = [c / s.ω for s in states] .!= Inf
    λ = 2π * minimum(λs[λs_nonzero])

    n_states = length(states)
    n_lasers = length(lasers)

    states = StructArray(states)
    lasers = StructArray(lasers)

    for i in eachindex(lasers)
        lasers.ω[i] = round_freq(lasers.ω[i], Γ, freq_res)
        lasers.ω[i] /= Γ
    end
    for i in eachindex(states)
        states.ω[i] = round_freq(states.ω[i], Γ, freq_res)
        states.ω[i] /= Γ
    end

    type_complex = ComplexF64
    type_real = Float64

    H       = StructArray( zeros(type_complex, n_states, n_states) )
    H_adj   = StructArray( zeros(type_complex, n_states, n_states) )
    # HJ      = zeros(type_real, n_states, n_states)
    HJ      = zeros(Float64, n_states)
    Hₒ      = [zeros(type_real, n_states, n_states) for l in lasers, q in qs]

    # Define the optical Hamiltonian; it has dimensions (k, q)
    # for k in eachindex(lasers)
    #     for s in eachindex(states), s′ in s:n_states
    #         q = Int64(states[s′].m - states[s].m)
    #         if abs(q) <= 1
    #             laser_ϵ = lasers[k].ϵ_re + im * lasers[k].ϵ_im
    #             print(laser_ϵ)
    #             print(ϵ[q+2])
    #             if !iszero(laser_ϵ ⋅ ϵ[q+2])
    #                 # Hₒ[k, q+2][s, s′] = d[s, s′, q+2]
    #                 # Hₒ[k, q+2][s′, s] = d[s, s′, q+2]
    #                 lasers[k].Hq[q+2][s, s′] = d[s, s′, q+2]
    #                 lasers[k].Hq[q+2][s′, s] = d[s, s′, q+2]
    #             end
    #         end
    #     end
    # end

    # Define the magnetic Hamiltonian; it has dimensions (q)
    # μ = (s, s′, q) -> - s.g * (-1)^(s.F - s′.m) * sqrt(s.F * (s.F + 1) * (2s.F + 1))
    #     * wigner3j(s.F, 1, s′.F, -s.m, q, s′.m)
    # for s in eachindex(states), s′ in eachindex(states), q in qs
    #     if states[s].F == states[s].F # States only mix if they belong to the same F state
    #         Hₘ[q+2][s, s′] = im * (-1)^q * μ(states[s], states[s′], -q)
    #     end
    # end

    # Construct an array containing all jump operators, as defined by `d`
    Js = Array{Jump}(undef, 0)
    for s′ in eachindex(states), s in s′:n_states, q in qs
        dme = d[s, s′, q+2]
        # println(dme)
        if dme != 0 & (states[s′].ω < states[s].ω) # only energy-allowed jumps are generated
            J = Jump(s, s′, dme)
            push!(Js, J)
        end
    end

    for J in Js
        # Adds the term - (iħ / 2) ∑ᵢ Jᵢ† Jᵢ
        # We assume jumps take the form Jᵢ = sqrt(Γ)|g⟩⟨e| such that Jᵢ† Jᵢ = Γ|e⟩⟨e|
        HJ[J.s] -= 0.5 * J.r^2
    end

    ω = [s.ω for s in states]
    eiωt  = StructArray(zeros(type_complex, n_states))

    dρ = deepcopy(ρ)
    ρ_soa = StructArray(ρ)
    dρ_soa = StructArray(dρ)

    # B(r) = [0,0,0]

    conj_mat = ones(Float64, n_states, n_states)
    for i in 1:n_states, j in 1:n_states
        if j < i
            conj_mat[i,j] = -1
        end
    end

    # Allocate some temporary arrays
    tmp1 = StructArray(zeros(ComplexF64, n_states, n_states))
    tmp2 = StructArray(zeros(ComplexF64, n_states, n_states))

     # Compute indices to indicate nonzero values in d
    d_nnz_m = Int64[]
    d_nnz_0 = Int64[]
    d_nnz_p = Int64[]
    d_m = @view d[:,:,1]
    d_0 = @view d[:,:,2]
    d_p = @view d[:,:,3]
    for i in 1:n_states^2
        if d_m[i] != 0
            push!(d_nnz_m, i)
        end
        if d_0[i] != 0
            push!(d_nnz_0, i)
        end
        if d_p[i] != 0
            push!(d_nnz_p, i)
        end
    end
    d_nnz = [d_nnz_m, d_nnz_0, d_nnz_p]

    particle = Particle([0.,0,0], [0.,0,0], [0.,0,0])

    t0 = [0.0, 0.0]

    p = @params (H, HJ, ρ, dρ, ρ_soa, dρ_soa, Js, ω, eiωt, states, lasers, particle, Γ, conj_mat, Hₒ, tmp1, tmp2, d, d_nnz, λ, t0)

    return (dρ, ρ, p)
end
export obe

function update_H_schrödinger!(τ, v, lasers, H)

    @turbo for l in eachindex(lasers)
        s = lasers.s[l]
        Ho_laser_m = lasers.Hqm[l]
        Ho_laser_0 = lasers.Hq0[l]
        Ho_laser_p = lasers.Hqp[l]
        for i in eachindex(H)
            H.re[i] += s * (Ho_laser_m[i] + Ho_laser_0[i] + Ho_laser_p[i])
        end
    end

    return nothing
end

function derivative_force(p, ρ, τ)

    @unpack ρ_soa, lasers, Γ, d, d_nnz = p

    r = p.particle.v .* τ
    update_lasers!(r, lasers, τ)

    F = SVector(0, 0, 0)

    for q in 1:3

        ampl = SVector(0, 0, 0)
        @inbounds for i in 1:length(lasers)
            s = lasers.s[i]
            k = lasers.k[i]
            ω = lasers.ω[i]
            x = h * Γ * s / (4π * √2)
            ampl += k * x * (im * lasers.f_re_q[i][q] - lasers.f_im_q[i][q])
        end

        d_q = @view d[:,:,q]
        d_nnz_q = d_nnz[q]
        @inbounds for i in d_nnz_q
            F += ampl * d_q[i] * ρ_soa[i] + conj(ampl * d_q[i] * ρ_soa[i])
            # F += ampl * d_q[i] * ρ[i] + conj(ampl * d_q[i] * ρ[i])
        end
    end

    return real(F[1])
end
export derivative_force

function force(p, ρ, τ)

    @unpack lasers, Γ, d, d_nnz = p

    r = p.particle.v .* τ
    update_lasers!(r, lasers, τ)

    base_to_soa!(ρ, p.ρ_soa)
    update_eiωt!(p.eiωt, p.ω, τ)
    Heisenberg!(p.ρ_soa, p.eiωt)

    F = SVector(0, 0, 0)

    for q in 1:3

        ampl = SVector(0, 0, 0)
        @inbounds for i in 1:length(lasers)
            s = lasers.s[i]
            k = lasers.k[i]
            ω = lasers.ω[i]
            # With SI units, should be x = -h * Γ * sqrt(s) / (2 * √2 * p.λ)
            x = sqrt(s) / (2 * √2)
            # x = h * Γ * sqrt(s) / (2 * √2 * p.λ)
            ampl += k * x * (im * lasers.f_re_q[i][q] - lasers.f_im_q[i][q])
        end
        # Note h * Γ / λ = 2π ħ * Γ / Λ = ħ * k * Γ, so this has units units of ħ k Γ

        d_q = @view d[:,:,q]
        d_nnz_q = d_nnz[q]
        # print(d_nnz_q)
        @inbounds for i in d_nnz_q
            # print(ampl * d_q[i] * p.ρ_soa[i] + conj(ampl * d_q[i] * p.ρ_soa[i]))
            F += ampl * d_q[i] * p.ρ_soa[i] + conj(ampl * d_q[i] * p.ρ_soa[i])
            # F += ampl * d_q[i] * ρ[i] + conj(ampl * d_q[i] * ρ[i])
        end
    end

    return real.(F[3])
end
export force

function update_lasers!(r, lasers, t)
    for i in 1:length(lasers)
        k = lasers.k[i]
        lasers.kr[i] = k ⋅ r
    end
    @turbo for i in 1:length(lasers)
        lasers.im[i], lasers.re[i] = sincos(lasers.kr[i] - lasers.ω[i] * t)
    end
    for i in 1:length(lasers)
        ϵ_re = lasers.ϵ_re[i]
        ϵ_im = lasers.ϵ_im[i]
        re = lasers.re[i]
        im = lasers.im[i]
        lasers.f_re_q[i] = ϵ_re * re - ϵ_im * im
        lasers.f_im_q[i] = ϵ_re * im + ϵ_im * re
    end
    return nothing
end

function update_H!(τ, r, Γ, lasers, H, conj_mat, d, d_nnz)

    update_lasers!(r, lasers, τ)

    @turbo for i in eachindex(H)
        H.re[i] = 0
        H.im[i] = 0
    end

    @inbounds for i in 1:length(lasers)
        s = lasers.s[i]
        x = sqrt(s) / (2 * √2)

        @inbounds for q in 1:3
            freq = lasers.f_re_q[i][q]
            fimq = lasers.f_im_q[i][q]
            if (freq > 1e-10) || (freq < -1e-10) || (fimq > 1e-10) || (fimq < -1e-10)
                @turbo for i in 1:size(H,1), j in 1:size(H,1)
                    H.re[i,j] += x * freq * d[i,j,q]
                    H.im[i,j] += x * fimq * d[i,j,q] * conj_mat[i,j]
                end
                # TODO: The approach below which iterates only over nonzero dipole moments is a bit faster
                # d_nnz_q = d_nnz[q]
                # d_q = @view d[:,:,q]
                # @inbounds @simd for i in d_nnz_q
                #     H.re[i] += x * freq * d_q[i]
                #     H.im[i] += x * fimq * d_q[i] * conj_mat[i]
                # end
            end
        end
    end

    # @inbounds @simd for q in qs
    #     β = (μB / (Γ * ħ)) * (B(r) ⋅ ϵ[q+2])
    #     axpy!(β, Hₘ[q+2], H)
    # end

    return nothing
end

function soa_to_base!(ρ::Array{<:Complex}, ρ_soa::StructArray{<:Complex})
    @inbounds for i in eachindex(ρ, ρ_soa)
        ρ[i] = ρ_soa.re[i] + im * ρ_soa.im[i]
    end
    return nothing
end

function base_to_soa!(ρ::Array{<:Complex}, ρ_soa::StructArray{<:Complex})
    @inbounds for i in eachindex(ρ, ρ_soa)
        ρ_soa.re[i] = real(ρ[i])
        ρ_soa.im[i] = imag(ρ[i])
    end
    return nothing
end

function update_eiωt!(eiωt::StructArray{<:Complex}, ω::Array{<:Real}, τ::Real)
    @turbo for i ∈ 1:size(ω, 1)
        eiωt.im[i], eiωt.re[i] = sincos( ω[i] * τ )
    end
    return nothing
end

function Heisenberg!(ρ::StructArray{<:Complex}, eiωt::StructArray{<:Complex}, im_factor=1)
    @inbounds for j ∈ 1:size(ρ, 2)
        jre = eiωt.re[j]
        jim = eiωt.im[j]
        for i ∈ 1:size(ρ, 1)
            ire = eiωt.re[i]
            iim = eiωt.im[i]
            cisim = im_factor * (iim * jre - ire * jim)
            cisre = ire * jre + iim * jim
            ρre_i = ρ.re[i,j]
            ρim_i = ρ.im[i,j]
            ρ.re[i,j] = ρre_i * cisre - ρim_i * cisim
            ρ.im[i,j] = ρre_i * cisim + ρim_i * cisre
        end
    end
    return nothing
end
export Heisenberg!

function im_commutator!(C, A, B, tmp1, tmp2, A_diag)
    @inbounds for i ∈ 1:size(A,1), j ∈ 1:size(B,2)
        Cre = 0.0
        Cim = 0.0
        for k ∈ 1:size(A,2)
            Aik_re = A.re[i,k]
            Aik_im = A.im[i,k]
            Bkj_re = B.re[k,j]
            Bkj_im = B.im[k,j]
            # Multiply by -1
            Cre -= Aik_re * Bkj_re - Aik_im * Bkj_im
            Cim -= Aik_re * Bkj_im + Aik_im * Bkj_re
        end
        C.re[i,j] = Cre
        C.im[i,j] = Cim
    end

    adjoint!(tmp1, C)
    C_add_A!(C, tmp1, -1)

    mul_by_im!(C)

    mul_diagonal!(tmp1, B, A_diag)
    adjoint!(tmp2, tmp1)

    C_add_AplusB!(C, tmp1, tmp2, 1, 1)

end
export im_commutator!

function ψ!(dψ, ψ, p, τ)

    # @avx p.H .= 0
    # @avx p.ψ_soa .= ψ
    base_to_soa!(ψ, p.ψ_soa)
    # Update the Hamiltonian according to the new time τ
    p.f(p.H, τ, p.p)
    # update_H_schrödinger!(τ, p.v, p.lasers, p.H)
    jgemvavx!(p.dψ_soa, p.H, p.ψ_soa)
    soa_to_base!(dψ, p.dψ_soa)

    return nothing
end
export ψ!

function ρ!(dρ, ρ, p, τ)

    p.particle.r = p.particle.r0 + p.particle.v .* τ

    base_to_soa!(ρ, p.ρ_soa)
    #p.ρ_soa .= ρ

    # Update the Hamiltonian according to the new time τ
    update_H!(τ, p.particle.r, p.Γ, p.lasers, p.H, p.conj_mat, p.d, p.d_nnz)

    # Apply a transformation to go to the Heisenberg picture
    update_eiωt!(p.eiωt, p.ω, τ)
    Heisenberg!(p.ρ_soa, p.eiωt)

    # Compute coherent evolution terms
    # im_commutator!(p.dρ_soa, p.H, p.ρ_soa, p.A12, p.B12, p.T1, p.T2, p.HJ, p.tmp1, p.tmp2)
    im_commutator!(p.dρ_soa, p.H, p.ρ_soa, p.tmp1, p.tmp2, p.HJ)

    # Add the terms ∑ᵢ Jᵢ ρ Jᵢ†
    # We assume jumps take the form Jᵢ = sqrt(Γ)|g⟩⟨e| such that JᵢρJᵢ† = Γ^2|g⟩⟨g|ρₑₑ
    @inbounds for i in eachindex(p.Js)
        J = p.Js[i]
        p.dρ_soa.re[J.s′, J.s′] += J.r^2 * p.ρ_soa.re[J.s, J.s]
    end

    # The left-hand side also needs to be transformed into the Heisenberg picture
    # To do this, we require the transpose of the `ω` matrix
    # Heisenberg!(p.dρ_soa, p.ω_trans, τ)
    Heisenberg!(p.dρ_soa, p.eiωt, -1)
    soa_to_base!(dρ, p.dρ_soa)

    return nothing
end
export ρ!

function mat_to_vec!(ρ, ρ_vec)
    @turbo for i in eachindex(ρ)
        ρ_vec[i] = ρ[i]
    end
    return nothing
end
export mat_to_vec!

function mat_to_vec_minus1!(ρ, ρ_vec)
    @turbo for i in 1:(length(ρ)-1)
        ρ_vec[i] = ρ[i]
    end
    return nothing
end
export mat_to_vec_minus1!

function ρ_and_force!(du, u, p, τ)

    p.particle.r = p.particle.v .* τ

    mat_to_vec_minus1!(u, p.ρ)
    base_to_soa!(p.ρ, p.ρ_soa)
    #p.ρ_soa .= ρ

    # Update the Hamiltonian according to the new time τ
    update_H!(τ, p.particle.r, p.Γ, p.lasers, p.H, p.conj_mat, p.d, p.d_nnz)

    # Apply a transformation to go to the Heisenberg picture
    update_eiωt!(p.eiωt, p.ω, τ)
    Heisenberg!(p.ρ_soa, p.eiωt)

    # Compute coherent evolution terms
    # im_commutator!(p.dρ_soa, p.H, p.ρ_soa, p.A12, p.B12, p.T1, p.T2, p.HJ, p.tmp1, p.tmp2)
    im_commutator!(p.dρ_soa, p.H, p.ρ_soa, p.tmp1, p.tmp2, p.HJ)

    # Add the terms ∑ᵢ Jᵢ ρ Jᵢ†
    # We assume jumps take the form Jᵢ = sqrt(Γ)|g⟩⟨e| such that JᵢρJᵢ† = Γ^2|g⟩⟨g|ρₑₑ
    @inbounds for i in eachindex(p.Js)
        J = p.Js[i]
        p.dρ_soa.re[J.s′, J.s′] += J.r^2 * p.ρ_soa.re[J.s, J.s]
    end

    # The left-hand side also needs to be transformed into the Heisenberg picture
    # To do this, we require the transpose of the `ω` matrix
    # Heisenberg!(p.dρ_soa, p.ω_trans, τ)
    Heisenberg!(p.dρ_soa, p.eiωt, -1)
    soa_to_base!(p.dρ, p.dρ_soa)

    mat_to_vec!(p.dρ, du)
    du[end] = derivative_force(p.ρ, p, τ)
    # u[end] = force(p.ρ, p, τ)

    return nothing
end
export ρ_and_force!

function C_copy_AplusB!(C::Array{<:Real}, A::Array{<:Real}, B::Array{<:Real}, α=1, β=1)
    @turbo for i in eachindex(A, B, C)
        C[i] = α * A[i] + β * B[i]
    end
end

function C_add_AplusB!(C::Array{<:Real}, A::Array{<:Real}, B::Array{<:Real}, α=1, β=1)
    @turbo for i in eachindex(A, B, C)
        C[i] += α * A[i] + β * B[i]
    end
end

function C_add_AplusB!(C::StructArray{<:Complex}, A::StructArray{<:Complex}, B::StructArray{<:Complex},
    α=1, β=1)
    @turbo for i in eachindex(A, B, C)
        C.re[i] += α * A.re[i] + β * B.re[i]
        C.im[i] += α * A.im[i] + β * B.im[i]
    end
end

function C_add_A!(C::StructArray{<:Complex}, A::StructArray{<:Complex}, factor=1)
    @turbo for i in eachindex(C, A)
        C.re[i] += factor * A.re[i]
        C.im[i] += factor * A.im[i]
    end
end

function C_copy_T₁T₂!(C::StructArray{<:Complex}, T1::Array{<:Real}, T2::Array{<:Real})
    @turbo for i ∈ eachindex(C)
        C.re[i] = T2[i] - T1[i]
        C.im[i] = T2[i] + T1[i]
    end
end

function mul_by_im!(C::StructArray{<:Complex})
    @turbo for i ∈ eachindex(C)
        a = C.re[i]
        C.re[i] = -C.im[i]
        C.im[i] = a
    end
end

function mul_diagonal!(C::StructArray{<:Complex}, A::StructArray{<:Complex}, B::Array{<:Real})
    """
    Computes A × B, where B is diagonal and real.
    """
    @turbo for j in axes(A,2)
        d = B[j]
        for i in axes(A,1)
            C.re[i,j] = d * A.re[i,j]
            C.im[i,j] = d * A.im[i,j]
        end
    end
end

function update_T₁T₂!(T1::Array{<:Real}, T2::Array{<:Real}, A::StructArray{<:Complex}, B::StructArray{<:Complex})
    @turbo for i ∈ 1:size(A,1), j ∈ 1:size(B,2)
        C1 = zero(eltype(A))
        C2 = zero(eltype(A))
        for k ∈ 1:size(A,2)
            C1 += A.re[i,k] * B.re[k,j]
            C2 += A.im[i,k] * B.im[k,j]
        end
        T1[i,j] = C1
        T2[i,j] = C2
    end
end

end


text:
export @code_yao, @code_qasm

function code_qasm(spec::RoutineSpec; optimize = false, gate = false, passes = Symbol[])
    ci, _ = code_yao(Semantic.main, spec; optimize, passes)
    target = gate ? TargetQASMGate() : TargetQASMTopLevel()
    return codegen(target, ci)
end

function code_qasm_m(exs...)
    options = []
    call = nothing

    for ex in exs
        if ex isa Expr && ex.head === :(=)
            option = ex.args[1]
            option in [:optimize, :gate, :passes] || error("invalid option: $option")
            push!(options, Expr(:kw, option, ex.args[2]))
        elseif ex isa Expr && ex.head === :call
            call = ex
        end
    end

    isnothing(call) && error("expect a routine call")

    return quote
        $code_qasm($call; $(options...))
    end
end

"""
    @code_qasm [optimize=false, gate=false, passes=:zx] <generic routine call>

Return the corresponding QASM code of given routine call. One can see the optimized
qasm by setting `optimize` to `true`. See [@code_yao](@ref) for detailed configuration on YaoCompiler passes
via `passes` keyword.

By default, the QASM code generation assumes the given routine evaluates
in toplevel. One can set `gate` to `true` to generate given routine as a
QASM gate routine.
"""
macro code_qasm(exs...)
    esc(code_qasm_m(exs...))
end

function code_yao(f, spec::RoutineSpec, args...; optimize::Bool = false, passes = Symbol[])
    if passes isa Symbol
        passes = [passes]
    end

    if optimize
        passes = isempty(passes) ? default_passes() : filter(x -> x !== :julia, passes)
        run_optimizer = true
    elseif :julia in passes
        run_optimizer = true
        passes = filter(x -> x !== :julia, passes)
    else
        if !isempty(passes)
            passes = filter(x -> x !== :julia, passes)
            run_optimizer = true
        else
            run_optimizer = false
        end
    end

    args_t = Base.typesof(spec, args...)
    atypes = Base.typesof(f, spec, args...)
    matches = methods(f, args_t)
    length(matches) == 1 || error("call is ambiguous")
    method = first(matches)
    mi = Core.Compiler.specialize_method(method, atypes, Core.svec())

    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    world = Core.Compiler.get_world_counter()
    interp = YaoCompiler.YaoInterpreter(; passes = passes)

    # NOTE: we need to run optimizer manually on current frame
    # if cache is false
    frame = Core.Compiler.InferenceState(result, false, interp) #=cached=#
    Core.Compiler.typeinf(interp, frame)
    if run_optimizer
        opt = OptimizationState(frame, OptimizationParams(interp), interp)
        Core.Compiler.optimize(interp, opt, OptimizationParams(interp), result.result)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, widenconst(result.result))
end

function code_yao_m(exs...)
    options = []
    call = nothing
    head = :main
    args = []
    for ex in exs
        if ex isa Expr && ex.head === :(=)
            ex.args[1] in [:optimize, :passes] || error("invalid option: $ex")
            push!(options, Expr(:kw, ex.args...))
        elseif ex isa Expr && ex.head === :call
            if ex.args[1] in (:gate, :ctrl)
                head = ex.args[1]
                call = ex.args[2]
                args = ex.args[3:end]
            else
                call = ex
            end
        end
    end

    return Expr(
        :call,
        code_yao,
        Expr(:parameters, options...),
        GlobalRef(Semantic, head),
        call,
        args...,
    )
end

"""
    @code_yao [ optimize=false， passes=[:julia, :zx] ] <generic routine call>

Return the quantum CodeInfo generated by YaoCompiler compiler for a given
generic routine call. One can also retreive the `gate`, `ctrl` call with
location mapping of a given routine specialization.

One can trun on optimization by setting `optimize=true`, this will call
all the optimization passes to simplify your program.
One can also configure passes to run via keyword `passes`.

## Valid Optimization Passes

- `:julia`: Julia's default optimization passes for classical programs
- `:zx`: YaoCompiler's ZX calculus based passes
"""
macro code_yao(exs...)
    return esc(code_yao_m(exs...))
end

ext: jl


text: @testset "Containers" begin include("containers.jl") end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################

# using pkgs
using Distributions
using HDF5
using Ipopt
using Polynomials
using PowerModels, PowerModelsDistribution
using PowerModelsDistributionStateEstimation
#using SCS #removed while SDP tests are not active
using Test

# pkg const
const _DST = Distributions
const _PMD = PowerModelsDistribution
const _PMDSE = PowerModelsDistributionStateEstimation
const _Poly = Polynomials

#network and feeder from ENWL for tests
ntw, fdr = 4, 2

season     = "summer"
time_step  = 144
elm        = ["load", "pv"]
pfs        = [0.95, 0.90]
rm_transfo = true
rd_lines   = true

# set solvers
ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer,"max_cpu_time"=>300.0,
                                                              "tol"=>1e-9,
                                                              "print_level"=>0)

# scs_solver = optimizer_with_attributes(SCS.Optimizer, "max_iters"=>20000, "eps"=>1e-5,
#                                                             "alpha"=>0.4, "verbose"=>0) #deactivated while SDP tests not active

@testset "PowerModelsDistributionStateEstimation" begin

    include("bad_data.jl")
    include("distributions.jl")
    include("estimation_criteria.jl")
    include("mixed_measurements.jl")
    include("non_exact_forms.jl")
    include("power_flow.jl")
    include("pseudo_measurements.jl")
    include("single_conductor_branches.jl")
    include("utils_and_start_val.jl")
    include("with_errors.jl")

end
ambiguities = Test.detect_ambiguities(PowerModelsDistributionStateEstimation);
if !isempty(ambiguities)
    println("ambiguities detected: $ambiguities")
end


text:
using Documenter, GtkObservables, TestImages
testimage("lighthouse")    # ensure all artifacts get downloaded before running tests

makedocs(sitename = "GtkObservables",
         format   = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         pages    = ["index.md", "controls.md", "drawing.md", "zoom_pan.md", "reference.md"]
         )

deploydocs(repo         = "github.com/JuliaGizmos/GtkObservables.jl.git",
           push_preview = true,
           )


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ ac3cc660-08e1-11eb-3db6-0d5949a89423
begin
	using Plots
	using Printf
	using CSV
	using DataFrames
end

# ╔═╡ 59909ca0-08e3-11eb-3348-41080f6a2bd8
md"""

## **_Hiemenz Flow_**

#### **Description:**

This notebook computes the Hiemenz Flow (two-dimensional stagnation point flow)

		hiemenz(ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

		hiemenz(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

#### **Hiemenz Equation**
Stagnation point flow problem towards an infinite flat plate can be represented using the ordinary differential equation (ODE) below:

$$f''' + f f'' - (f')² + 1 = 0$$

Boundary Conditions:

$$f(0) = 0$$
$$f'(0) = 0$$
$$f'(∞) = 1$$

#### **Solution Method:**

Let's assume:

$$p = f(η)$$
$$h = f'(η)$$

Then, we obtain two ODEs equations and boundary conditions
ODEs:

$$h'' + p h' - h² + 1 = 0$$
$$p' - h = 0$$

BCs:

$$p(0) = 0$$
$$h(0) = 0$$
$$h(∞) = 1$$

The equations given above are solved using Thomas algorithm.

Feel free to ask questions!

*Furkan Oz*
*Kursat Kara*

[foz@okstate.edu](foz@okstate.edu)
[kursat.kara@okstate.edu](kursat.kara@okstate.edu)

"""

# ╔═╡ d47c93d0-08e1-11eb-0e92-21ca5523e6c6
"""
Computes the Hiemenz Flow (Two-dimensional stagnation point flow)

    hiemenz(ηmax, N, itermax, ϵProfile, ϵBC)

If the arguments are missing, it will use the default values.

    hiemenz(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

Furkan Oz,
[foz@okstate.edu](foz@okstate.edu),

"""
function hiemenz(ηmax=10, N=30, itermax=40, ϵProfile=1e-6, ϵBC=1e-6)

	Δη = ηmax/N
    Δη²= Δη^2;

    #itermax = 40
    #ϵProfile = 1e-6
    #ϵBC = 1e-6

    iter = 0
    errorProfile = 1.
    errorBC = 1.;

    # Initialization of the arrays
    A = zeros(N+1)
    B = zeros(N+1)
    C = zeros(N+1)
    D = zeros(N+1)

    G = zeros(N+1)
    H = zeros(N+1)

    p = zeros(N+1)
    h = zeros(N+1)
    η = zeros(N+1)

    η = [(i-1)*Δη for i=1:N+1];

    #BCs
    h[1]   = 0.0   # h(η=0) = 0
    h[N+1] = 1.0   # h(η=∞) = 1
    p[1]   = 0.0   # p(η=0) = 0

    G[1] = 0.0
    H[1] = 0.0;

    #Solution initialization

    # assume a linear profile for h(η)
    h = [(i-1)/N  for i=1:N+1]

    for i = 2:N+1
        p[i] = p[i-1] + (h[i] + h[i-1])*Δη/2
    end

    println("iter     error            convergence of h(η→∞)")
    println("-----------------------------------------------")

    #while ϵProfile<=errorProfile && ϵBC<=errorBC && iter<itermax
    while ϵProfile<=errorProfile && iter<itermax
        A = [ 1/Δη² + p[i]/(2*Δη) for i=1:N+1]
        B = [-2/Δη² - h[i] for i=1:N+1]
        C = [ 1/Δη² - p[i]/(2*Δη) for i=1:N+1]
        D = [1 for i=1:N+1]

        for i=2:N
            G[i] = - ( C[i]*G[i-1] + D[i] )/(B[i] + C[i] * H[i-1])
            H[i] = -                 A[i]  /(B[i] + C[i] * H[i-1])
        end

        hold = copy(h)

        for i=N:-1:2
            h[i] = G[i] + H[i] * h[i+1]
        end

        errorProfile = maximum(abs.(hold-h))

        errorBC = abs(h[N+1]-h[N])

        for i = 2:N
            p[i] = p[i-1] + (h[i] + h[i-1])*Δη/2
        end

        iter += 1

        #println("iter: $iter,  errorProfile: $errorProfile, errorBC: $errorBC")
        @printf("%4.4d %16.6e %16.6e \n", iter, errorProfile, errorBC)
    end

     if errorProfile<=ϵProfile
        println("")
        println("Solution converged!")
        println("The maximum change between consecutive profiles is less than the error criteria ϵProfile=$ϵProfile.")
     end

     if errorBC<=ϵBC
        println("")
        println("Solution for the boundary condition converged!")
        println("The difference between h(N) and h(N+1) is less than the error criteria ϵBC=$ϵBC.")
     end
    return η,h
end

# ╔═╡ df424c60-08e1-11eb-319f-2fa7a0f0fb08
hiemenz();

# ╔═╡ 07798e4e-08e2-11eb-3066-0120057d15b2
ηtest, htest = hiemenz(10,40);

# ╔═╡ f0a39bd0-08e1-11eb-3d16-ffb6b9e631e0
begin
	data = CSV.read("hiemenz-white-dataset.csv",DataFrame)
	ηref = convert(Array,data[1:2:end,1])
	uref = convert(Array,data[1:2:end,2])
	plot_ref = plot([htest,uref],[ηtest,ηref],
        title = "Hiemenz Similarity Profile",
        label = ["Thomas Algorithm" "White(2006)"],
        legend = :topleft,
        xlabel = "h = f '",
        ylabel = "\\eta",
        linewidth = [2.0 0.000],
        linecolor = :auto,
        markershape = :auto,
		markeralpha = [1.0 0.6],
        markercolor = :auto,
	)
	#png(plot_ref, "hiemenzprofile")
end

# ╔═╡ 638b8090-08e2-11eb-2c26-7b7f6caa7d14
md"""
#### The effect of ηmax:
N is selected as 30 and the error limit is 1e-6.
"""

# ╔═╡ 7bbba7b0-08e9-11eb-3d8c-75a72baa41fe
begin
	η5, h5 = hiemenz(5);
	η10, h10 = hiemenz(10);
	η20, h20 = hiemenz(20);
end

# ╔═╡ d8a1b2d0-08e9-11eb-1c9a-91c4d297b41e
plot([h5, h10, h20], [η5, η10, η20],
        title = "Hiemenz Flow - Effect of Viscous Layer Height",
        label = ["\\etamax = 5" "\\etamax = 10" "\\etamax = 20"],
        legend = :topleft,
        xlabel = "h",
        ylabel = "\\eta",
        linewidth = 1,
        linecolor = :black,
        markershape = [:x :utriangle :dtriangle],
        markersize = [5 5 5],
        markeralpha = [0.9 0.6 0.6],
        markercolor = [:black :red :green]
    )

# ╔═╡ 2a95ec50-08ea-11eb-06a4-050c654d9c1f
md"""
#### The effect of N:
ηmax is selected as 10 and the error limit is 1e-6.
"""

# ╔═╡ 3649fc80-08ea-11eb-28b5-f58a62c010da
begin
	ηN10, hN10 = hiemenz(10, 10);
	ηN20, hN20 = hiemenz(10, 20);
	ηN40, hN40 = hiemenz(10, 40);
	ηN80, hN80 = hiemenz(10, 80);
end

# ╔═╡ 5241b720-08ea-11eb-3fad-5f8164fe40f6
begin
	gr()
plot([hN10, hN20, hN40, hN80], [ηN10, ηN20, ηN40, ηN80],
        seriestype = :scatter,
        title = "Hiemenz Flow - Effect of Grid Resolution",
        label = ["N = 10" "N = 20" "N = 40" "N = 80"],
        legend = :topleft,
        xlabel = "h",
        ylabel = "\\eta",
        markershape = [:circle :utriangle :dtriangle :square],
        markersize = [10 8 6 3],
        markeralpha = [0.6 0.6 0.6 0.6],
        markercolor = [:purple :red :green :blue]
    )
end

# ╔═╡ Cell order:
# ╟─59909ca0-08e3-11eb-3348-41080f6a2bd8
# ╠═ac3cc660-08e1-11eb-3db6-0d5949a89423
# ╠═d47c93d0-08e1-11eb-0e92-21ca5523e6c6
# ╠═df424c60-08e1-11eb-319f-2fa7a0f0fb08
# ╠═07798e4e-08e2-11eb-3066-0120057d15b2
# ╠═f0a39bd0-08e1-11eb-3d16-ffb6b9e631e0
# ╟─638b8090-08e2-11eb-2c26-7b7f6caa7d14
# ╠═7bbba7b0-08e9-11eb-3d8c-75a72baa41fe
# ╠═d8a1b2d0-08e9-11eb-1c9a-91c4d297b41e
# ╟─2a95ec50-08ea-11eb-06a4-050c654d9c1f
# ╠═3649fc80-08ea-11eb-28b5-f58a62c010da
# ╠═5241b720-08ea-11eb-3fad-5f8164fe40f6

max_stars_repo_name: frkanz/A-CFD-Tutorial-in-Julia

max_issues_repo_name: frkanz/A-CFD-Tutorial-in-Julia

max_forks_repo_name: frkanz/A-CFD-Tutorial-in-Julia


text:
using DelimitedFiles

Ls = [4,6,8,10,12,14,16]
nlayers = [10, 50, 100, 200, 500, 1000, 2000, 5000]
repeats = [i for i in 1:4]

function read()
    Es = zeros(length(Ls), length(nlayers), length(repeats))
    i = 0
    for (k, L) in enumerate(Ls)
        for (i, nlayer) in enumerate(nlayers)
            for (j, r) in enumerate(repeats)
                E = 0
                try
                    E =
                    readdlm("/scratch/mbeach/new_new_tickettoride/L-$L/layers-$nlayer-r-$r.txt")[end]
                catch
                end
                # println(E)
                # println(k, "  ", i, "  ",  j)
                Es[k, i, j] = E
            end
        end
    end
    return Es
end

Eex = [-0.4040063509461097,
       -0.4155961889813218,
       -0.42186657483598655,
        -0.4258035207282875,
        -0.4285075527367124]
# Eex = [Eex[i] for i in 1:5, for j in 1:10]
Es = read()
Emins = dropdims(minimum(Es, dims=3), dims=3)


text:
export DotProduct, dot

"""
Description:

    out = in1'*in2

    in1: d-dimensional vector
    in2: d-dimensional vector
    out: scalar

           in2
           |
      in1  V   out
    ----->[⋅]----->

    f(out, in1, in2) =  δ(out - in1'*in2)

Interfaces:

    1 i[:out], 2 i[:in1], 3 i[:in2]

Construction:

    DotProduct(out, in1, in2, id=:my_node)
"""
mutable struct DotProduct <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function DotProduct(out, in1, in2; id=generateId(DotProduct))
        @ensureVariables(out, in1, in2)
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:in2] = self.interfaces[3] = associate!(Interface(self), in2)

        return self
    end
end

slug(::Type{DotProduct}) = "dot"

function dot(a::Variable, in1::Variable)
    out = Variable()
    DotProduct(out, in1, a)
    return out
end

function dot(a::Variable, in1)
    @ensureVariables(in1)
    dot(a, in1)
end

function dot(a, in1::Variable)
    @ensureVariables(a)
    dot(a, in1)
end


text:
module asteroids

export Point, tau, chart, maxima, sweep, destroy

const tau = 2 * atan(0, -1)

struct Point
    x::Int64
    y::Int64
end

function chart(map)
    coords = Set{Point}()
    for i = 1:length(map)
        for j = 1:length(map[i])
            if map[i][j] == '#'
                push!(coords, Point(j, i))
            end
        end
    end
    return coords
end

function angle(origin, target)
    r = atan(target.y - origin.y, target.x - origin.x)
    return r < atan(1, 0) ? r + tau : r
end

function visible(coords, origin)
    delete!(coords, origin)
    return length(unique([
        angle(origin, target) for target in coords
    ]))
end

function maxima(coords)
    acc = Dict()
    for ast in coords
        acc[ast] = visible(copy(coords), ast)
    end
    max = maximum(values(acc))
    return (first(keys(filter(ast -> ast.second == max, acc))), max)
end

function distance(origin, target)
    return abs(target.x - origin.x) + abs(target.y - origin.y)
end

function closest(origin, left, right)
    if left === nothing && right === nothing
        return nothing
    elseif left === nothing
        return right
    elseif right === nothing
        return left
    else
        return distance(origin, left) > distance(origin, right) ?
            right :
            left
    end
end

function sweep(coords, origin, start, finish)
    acc = Dict()
    for ast in coords
        dx = ast.x - origin.x
        dy = ast.y - origin.y
        r = atan(dy, dx)
        # rotate anything not in the -x,+y quadrant by tau so they sort from
        # 'north' (atan(-11, 0)) clockwise from highest to lowest
        nr = r < atan(-1, 0) ? r + tau : r
        current = get(acc, nr, nothing)
        acc[nr] = closest(origin, current, ast)
    end
    directions = keys(filter(entry -> entry.first >= finish && entry.first < start, acc))
    result = [acc[direction] for direction in sort(collect(directions))]
    return result
end

function destroy(coords, origin, start, finish)
    inrange = sweep(coords, origin, start, finish)
    [delete!(coords, ast) for ast in inrange]
end

end


text:
include("base.jl")
include("blas.jl")
include("linalg.jl")
include("statistics.jl")
include("nnlib.jl")
include("sparse.jl")
include("mapreduce.jl")
include("bennett.jl")


text:
# ------------ Type alias and APIs ------------ #

const Type = MLIR.API.MlirType

get_context(t::Type) = MLIR.API.mlirTypeGetContext(t)
is_null(t::Type) = MLIR.API.mlirTypeIsNull(t)
Base.:(==)(t1::Type, t2::Type) = MLIR.API.mlirTypeEqual(t1, t2)
dump(t::Type) = MLIR.API.mlirTypeDump(t)
parse_type(ctx::Context, t::String) = MLIR.API.mlirTypeParseGet(ctx, StringRef(t))

# Constructor.
Type(ctx::Context, t::String) = parse_type(ctx, t)

@doc(
"""
const Type = MLIR.API.MlirType
""", Type)

# ------------ Builtins ------------ #

# Integers.
is_int(t::Type) = MLIR.API.mlirTypeIsAInteger(t)
get_int_type(ctx::Context, bitwidth::Int64) = MLIR.API.mlirIntegerTypeGet(ctx, bitwidth)
get_sint_type(ctx::Context, bitwidth::Int64) = MLIR.API.mlirIntegerSignedGet(ctx, bitwidth)
get_uint_type(ctx::Context, bitwidth::Int64) = MLIR.API.mlirIntegerUnsignedGet(ctx, bitwidth)
get_bitwidth(t::Type) = MLIR.API.mlirIntegerTypeGetWidth(t)

# Floating point.
get_f32_type(ctx::Context) = MLIR.API.mlirF32TypeGet(ctx)
get_f64_type(ctx::Context) = MLIR.API.mlirF64TypeGet(ctx)


text:
"""
    image(cmd0::String="", arg1=nothing; kwargs...)

Place images or EPS files on maps.

Full option list at [`psimage`]($(GMTdoc)image.html)

Parameters
----------

- $(GMT.opt_B)
- **D** | **pos** | **position** :: [Type => Str]

    Sets reference point on the map for the image using one of four coordinate systems.
    ($(GMTdoc)image.html#d)
- **F** | **box** :: [Type => Str | []]

    Without further options, draws a rectangular border around the image using MAP_FRAME_PEN.
    ($(GMTdoc)image.html#f)
- **I** | **invert_1bit** :: [Type => Str | Number]

    Invert 1-bit image before plotting.
    ($(GMTdoc)image.html#i)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- **M** | **monochrome** :: [Type => Bool]

    Convert color image to monochrome grayshades using the (television) YIQ-transformation.
    ($(GMTdoc)image.html#m)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_savefig)
"""
function image(cmd0::String="", arg1=nothing; first=true, kwargs...)

	gmt_proggy = (IamModern[1]) ? "image "  : "psimage "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd = parse_BJR(d, "", "", O, " -JX" * split(def_fig_size, '/')[1] * "/0")[1]
	cmd = parse_common_opts(d, cmd, [:F :UVXY :JZ :c :p :t :params], first)[1]
	cmd = parse_these_opts(cmd, d,  [[:G :bit_color], [:I :invert_1bit], [:M :monochrome]])
	cmd = parse_type_anchor(d, cmd, [:D :pos :position],
	                        (map=("g", arg2str, 1), outside=("J", arg2str, 1), inside=("j", arg2str, 1), norm=("n", arg2str, 1), paper=("x", arg2str, 1), anchor=("", arg2str, 2), dpi="+r", width=("+w", arg2str), justify="+j", replicate=("+n", arg2str), offset=("+o", arg2str)), 'j')

	cmd, _, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted

	return finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1)
end

# ---------------------------------------------------------------------------------------------------
image!(cmd0::String="", arg1=nothing; kw...) = image(cmd0, arg1; first=false, kw...)

# ---------------------------------------------------------------------------------------------------
const psimage  = image			# Alias
const psimage! = image!			# Alias


text:
## Pooling API
#
#  We provide the following generic methods, for 3d, 4d, and 5d tensors, calculating 1d,
#  2d and 3d pooling, based on the rank of the input tensors, in both mutating and
#  non-mutating auto-allocating variants:
#   - Pooling:
#     - maxpool(x, pdims)
#     - maxpool!(y, x, pdims)
#     - meanpool(x, pdims)
#     - meanpool!(y, x, pdims)
#   - Pooling input backprop
#     - ∇maxpool(dy, y, x, pdims)
#     - ∇maxpool!(dx, dy, y, x, pdims)
#     - ∇meanpool(dy, y, x, pdims)
#     - ∇meanpool!(dx, dy, y, x pdims)
#
#   All methods require a `PoolDims` object to define the dimensions and optional
#   elements of the convolution (stride, dilation, etc...), which is easily constructable
#   through something like `PoolDims(x, w)`.

# First, we will define mappings from the generic API names to our accelerated backend
# implementations.  At the moment this is only the direct implementation, however this
# exists here so that other packages (NNPACK, MAGMA, etc...) can override this easily.
for (front_name, backend) in (
        # This maps from public, front-facing name, to internal backend name
        :maxpool  => :direct,
        :meanpool => :direct,
    )

    # We only define 3d pooling primitives, we reshape lower down to get 1d and 2d pooling
    @eval begin
        function $(Symbol("$(front_name)!"))(
                y::AbstractArray{T,5}, x::AbstractArray{T,5},
                pdims::PoolDims; kwargs...) where {T}
            $(Symbol("$(front_name)_$(backend)!"))(y, x, pdims; kwargs...)
        end
    end
end

# Do the same for backprops
for (front_name, backend) in (
        :∇maxpool  => :direct,
        :∇meanpool => :direct,
    )
    @eval begin
        function $(Symbol("$(front_name)!"))(
                        dx::AbstractArray{T,5}, dy::AbstractArray{T,5},
                        y::AbstractArray{T,5}, x::AbstractArray{T,5},
                        pdims::PoolDims; kwargs...) where {T}
            $(Symbol("$(front_name)_$(backend)!"))(dx, dy, y, x, pdims; kwargs...)
        end
    end
end

# Our strategy for pooling is to reshape to an array with three spatial dimensions, which
# makes things MUCH EASIER for us on the backend side, and is in general pretty fast,
# since we can specialize on sizes.
for front_name in (:maxpool, :meanpool)
    for backend in (Symbol(), :_direct)
        for N in (3, 4)
            @eval begin
                function $(Symbol("$(front_name)$(backend)!"))(
                                y::AbstractArray{T,$N}, x::AbstractArray{T,$N},
                                pdims::PoolDims; kwargs...) where {T}
                    $(Symbol("$(front_name)$(backend)!"))(
                        insert_singleton_spatial_dimension(y, $(5 - N)),
                        insert_singleton_spatial_dimension(x, $(5 - N)),
                        insert_singleton_spatial_dimension(pdims, $(5 - N));
                        kwargs...
                    )

                    # We explicitly return `y` here, because the backend call
                    # itself may return a reshaped view, which we don't want.
                    return y
                end

                # backprops too
                function $(Symbol("∇$(front_name)$(backend)!"))(
                                dx::AbstractArray{T,$N}, dy::AbstractArray{T,$N},
                                y::AbstractArray{T,$N}, x::AbstractArray{T,$N},
                                pdims::PoolDims; kwargs...) where {T}
                    $(Symbol("∇$(front_name)$(backend)!"))(
                        insert_singleton_spatial_dimension(dx, $(5 - N)),
                        insert_singleton_spatial_dimension(dy, $(5 - N)),
                        insert_singleton_spatial_dimension(y, $(5 - N)),
                        insert_singleton_spatial_dimension(x, $(5 - N)),
                        insert_singleton_spatial_dimension(pdims, $(5 - N));
                        kwargs...
                    )

                    # We explicitly return `dx` here, because the backend call
                    # itself may return a reshaped view, which we don't want.
                    return dx
                end
            end
        end
    end
end

# Finally, let's generate auto-allocating versions of all our functions, for all backends:
for backend in (Symbol(), :_direct, :_nnpack)
    # First make auto-allocating versions of the basic pooling calls:
    for name in (:maxpool, :meanpool)
        @eval begin
            function $(Symbol("$(name)$(backend)"))(
                            x::AbstractArray{xT,N},
                            pdims::PoolDims; kwargs...) where {xT, N}
                y = similar(x, output_size(pdims)..., channels_out(pdims), size(x, N))
                fill!(y, xT(0))
                return $(Symbol("$(name)$(backend)!"))(y, x, pdims; kwargs...)
            end

            # Backprops too
            function $(Symbol("∇$(name)$(backend)"))(
                            dy::AbstractArray{T,N}, y::AbstractArray{T,N},
                            x::AbstractArray{T,N}, pdims::PoolDims;
                            kwargs...) where {T, N}
                dx = similar(x, input_size(pdims)..., channels_in(pdims), size(dy, N))
                fill!(dx, T(0))
                return $(Symbol("∇$(name)$(backend)!"))(dx, dy, y, x, pdims; kwargs...)
            end
        end
    end
end

## Use NNPACK if it is available and operation is supported.
## The corresponding gradient is not available in NNPACK
## Commented out due to #210
# if is_nnpack_available()
#     function maxpool(x::Array{Float32, 4}, pdims::PoolDims{2, K, S, P, (1, 1)}; kwargs...) where {T, K, S, P}
#         func = nnpack_supported_operation(pdims) ? maxpool_nnpack : maxpool_direct
#         return func(x, pdims; kwargs...)
#     end
# end

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

"""
    maxpool(x, k::NTuple; pad=0, stride=k)

Perform max pool operation with window size `k` on input tensor `x`.
"""
function maxpool(x, k::NTuple{N, Integer}; pad=0, stride=k) where N
    pad = expand(Val(N), pad)
    stride = expand(Val(N), stride)
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return maxpool(x, pdims)
end

"""
    meanpool(x, k::NTuple; pad=0, stride=k)

Perform mean pool operation with window size `k` on input tensor `x`.
"""
function meanpool(x, k::NTuple{N, Integer}; pad=0, stride=k) where N
    pad = expand(Val(N), pad)
    stride = expand(Val(N), stride)
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return meanpool(x, pdims)
end

for pool in [:maxpool, :meanpool]
    ∇pool = Symbol(:∇, pool)
    pullback = Symbol(pool, :_pullback)
    @eval function rrule(::typeof($pool), x, pdims::PoolDims; kw...)
        Ω = $pool(x, pdims; kw...)
        $pullback(Δ) = (NoTangent(), $∇pool(unthunk(Δ), Ω, x, pdims; kw...), NoTangent())
        return Ω, $pullback
    end
end


text:
module TestCG

using IterativeSolvers
using LinearMaps
using Test
using LinearAlgebra
using SparseArrays
using Random

import LinearAlgebra.ldiv!

include("laplace_matrix.jl")

struct JacobiPrec
    diagonal
end

ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal

@testset "Conjugate Gradients" begin

Random.seed!(1234321)

@testset "Small full system" begin
    n = 10

    @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T, n, n)
        A = A' * A + I
        b = rand(T, n)
        reltol = √eps(real(T))

        x,ch = cg(A, b; reltol=reltol, maxiter=2n, log=true)
        @test isa(ch, ConvergenceHistory)
        @test norm(A*x - b) / norm(b) ≤ reltol
        @test ch.isconverged

        # If you start from the exact solution, you should converge immediately
        x,ch = cg!(A \ b, A, b; abstol=2n*eps(real(T)), reltol=zero(real(T)), log=true)
        @test niters(ch) ≤ 1
        @test nprods(ch) ≤ 2

        # Test with cholfact as preconditioner should converge immediately
        F = cholesky(A, Val(false))
        x,ch = cg(A, b; Pl=F, log=true)
        @test niters(ch) ≤ 2
        @test nprods(ch) ≤ 2

        # All-zeros rhs should give all-zeros lhs
        x0 = cg(A, zeros(T, n))
        @test x0 == zeros(T, n)
    end
end

@testset "Sparse Laplacian" begin
    A = laplace_matrix(Float64, 10, 2)
    P = JacobiPrec(diag(A))

    rhs = randn(size(A, 2))
    rmul!(rhs, inv(norm(rhs)))
    abstol = 1e-5
    reltol = 1e-5

    @testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, Float32), Ti in (Int64, Int32)
        xCG = cg(A, rhs; reltol=reltol, maxiter=100)
        xJAC = cg(A, rhs; Pl=P, reltol=reltol, maxiter=100)
        @test norm(A * xCG - rhs) ≤ reltol
        @test norm(A * xJAC - rhs) ≤ reltol
    end

    Af = LinearMap(A)
    @testset "Function" begin
        xCG = cg(Af, rhs; reltol=reltol, maxiter=100)
        xJAC = cg(Af, rhs; Pl=P, reltol=reltol, maxiter=100)
        @test norm(A * xCG - rhs) ≤ reltol
        @test norm(A * xJAC - rhs) ≤ reltol
    end

    @testset "Function with specified starting guess" begin
        x0 = randn(size(rhs))
        xCG, hCG = cg!(copy(x0), Af, rhs; abstol=abstol, reltol=0.0, maxiter=100, log=true)
        xJAC, hJAC = cg!(copy(x0), Af, rhs; Pl=P, abstol=abstol, reltol=0.0, maxiter=100, log=true)
        @test norm(A * xCG - rhs) ≤ reltol
        @test norm(A * xJAC - rhs) ≤ reltol
        @test niters(hJAC) == niters(hCG)
    end
end

@testset "CG with a view" begin
    A = rand(10, 10)
    A = A + A' + 100I
    x = view(rand(10, 2), :, 1)
    b = rand(10)
    x, hist = cg!(x, A, b, log = true)
    @test hist.isconverged
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = cg!(x, A, b, log=true)
        @test 2 ≤ niters(ch) ≤ n

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = cg!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end

end

end # module


text:
function costcLshift(U_rad,
    S_U::Vector{Complex{T}},
    updatedfunc, updateβfunc, updateλfunc, updateκfunc,
    p::Vector{T},
    Es,
    E_BLS::Matrix{Complex{T}}, κ_BLS::Vector{T}, b_BLS,
    f)::T where T <: Real

    #println("p = ", p)
    ### update main optim variables.
    updatedfunc(p)
    updateβfunc(p)
    updateλfunc(p)

    updateκfunc(p)
    parseκ!(Es, κ_BLS)

    # ## l-2 costfunc.
    # cost = zero(T)
    # for m = 1:length(S_U)
    #
    #     f_u = f(U_rad[m])
    #
    #     cost += abs2( f_u - S_U[m] )
    #     #cost += (abs2(f_u) - abs2(S_U[m]))^2
    #     #cost += abs2( f_buffer[m] - DFT_s[m] ) + (abs2(f_buffer[m]) - abs2(DFT_s[m]))^2
    # end
    B = reinterpret(T, E_BLS)
    tmp = B*κ_BLS - b_BLS
    cost = dot(tmp, tmp)

    return cost
end

# WIP.
function setupcostalign0(Es::Vector{NMRSpectraSimulator.καFIDModelType{T, SST}},
    Bs::Vector{NMRSpectraSimulator.FIDModelType{T, SST}},
    As::Vector{NMRSpectraSimulator.SHType{T}},
    κs_β_orderings::Vector{Vector{Vector{Int}}},
    κs_β_DOFs::Vector{Vector{Int}},
    fs::T,
    SW::T,
    LS_inds,
    U0_rad,
    y0::Vector{Complex{T}},
    Δsys_cs::Vector{Vector{T}};
    w = ones(T, length(Es))) where {T <: Real, SST}

    # model.
    f = uu->NMRSpectraSimulator.evalitpproxymixture(uu, Bs, Es; w = w)

    ##### update functions.
    N_d = sum( getNd(Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

    st_ind_d = 1
    fin_ind_d = st_ind_d + N_d - 1
    updatedfunc = pp->updatemixtured!(Bs,
        p, st_ind, fs, SW, Δsys_cs)

    st_ind_β = fin_ind_d + 1
    fin_ind_β = st_ind_β + N_β - 1
    updateβfunc = pp->updateβ!(Bs, pp, st_ind_β)

    N_vars_set = [N_d; N_β]

    ### LS κ.
    U_rad_LS = U0_rad[LS_inds]

    N_κ, N_κ_singlets = countκ(Es)
    N_κ_vars = N_κ + N_κ_singlets
    E_BLS, κ_BLS, b_BLS = setupupdateLS(length(U_rad_LS), N_κ_vars, y0[LS_inds])

    κ_lb = ones(N_κ_vars) .* κ_lb_default
    κ_ub = ones(N_κ_vars) .* κ_ub_default

    updateκfunc = xx->updateκ!(E_BLS, b_BLS, κ_BLS,
    U_rad_LS, Es, w, κ_lb, κ_ub)

    #### extract parameters from p.
    getshiftfunc = pp->pp[st_ind_d:fin_ind_d]
    getβfunc = pp->pp[st_ind_β:fin_ind_β]
    getλfunc = pp->pp[st_ind_λ:fin_ind_λ]

    return f, updatedfunc, updateβfunc, updateλfunc, updateκfunc,
    E_BLS, κ_BLS, b_BLS, getshiftfunc, getβfunc, getλfunc, N_vars_set
end


text:
using BandedMatrices, LinearAlgebra, ArrayLayouts, Test
import BandedMatrices: BandedColumns, BandedRows

@testset "Triangular" begin
    @testset "Upper" begin
        for (t, ud) in ((UpperTriangular, 'N'), (UnitUpperTriangular, 'U'))
            A = t(brand(10,10,1,2))
            @test isbanded(A)
            @test MemoryLayout(typeof(A)) == TriangularLayout{'U',ud,BandedColumns{DenseColumnMajor}}()
            @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (0,2)
            @test bandwidth(A,2) == 2
            @test BandedMatrix(A) == A
            @test A[2,1] == 0
            if( ud == 'U' )
                @test A[2,2] == 1
            end
            x=rand(10)
            @test A*x ≈ BandedMatrix(A)*x ≈ Matrix(A)*x
            @test transpose(A)*x ≈ transpose(BandedMatrix(A))*x ≈ transpose(Matrix(A))*x
            @test all(A*x .=== lmul!(A, copy(x)) .=== lmul(A,x) .===
                        BandedMatrices.tbmv!('U', 'N', ud, 10, 2, parent(A).data, copy(x)))
            @test_throws DimensionMismatch BandedMatrices.tbmv('U', 'N', ud, 10, 2, parent(A).data, rand(9))

            @test A\x ≈ Matrix(A) \ x ≈ BandedMatrix(A) \ x
            @test all((A \ x) .=== ldiv!(A, copy(x)) .=== copyto!(similar(x), Ldiv(A, x)) .=== (similar(x) .= Ldiv(A, x)) .===
                        BandedMatrices.tbsv!('U', 'N', ud, 10, 2, parent(A).data, copy(x)))
            @test_throws DimensionMismatch BandedMatrices.tbsv('U', 'N', ud, 10, 2, parent(A).data, rand(9))
        end
    end
    @testset "Lower" begin
        for (t, ud) in ((LowerTriangular, 'N'), (UnitLowerTriangular, 'U'))
            A = t(brand(10,10,1,2))
            @test isbanded(A)
            @test MemoryLayout(typeof(A)) == TriangularLayout{'L', ud,BandedColumns{DenseColumnMajor}}()
            @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (1,0)
            @test bandwidth(A,2) == 0
            @test BandedMatrix(A) == A
            @test A[1,2] == 0
            if( ud == 'U' )
                @test A[2,2] == 1
            end

            x=rand(10)
            @test A*x ≈ Matrix(A)*x
            @test transpose(A)*x ≈ transpose(BandedMatrix(A))*x ≈ transpose(Matrix(A))*x
            @test all(A*x .=== lmul!(A, copy(x)) .=== lmul(A,x) .===
                        BandedMatrices.tbmv!('L', 'N', ud, 10, 1, view(parent(A).data, 3:4,:), copy(x)))
            @test_throws DimensionMismatch BandedMatrices.tbmv('L', 'N', ud, 10, 1, view(parent(A).data, 3:4,:), rand(9))

            @test A\x ≈ Matrix(A) \ x ≈ BandedMatrix(A) \ x
            @test all((A \ x) .=== ldiv!(A, copy(x)) .=== copyto!(similar(x), Ldiv(A, x)) .=== (similar(x) .= Ldiv(A, x)) .===
                        BandedMatrices.tbsv!('L', 'N', ud, 10, 1, view(parent(A).data, 3:4,:), copy(x)))
            @test_throws DimensionMismatch BandedMatrices.tbsv('L', 'N', ud, 10, 1, view(parent(A).data, 3:4,:), rand(9))
        end
    end

    @testset "subarray" begin
        A = brand(6,6,2,1)
        L = LowerTriangular(A)
        U = UpperTriangular(A)
        V = view(L,2:4,1:3)
        @test MemoryLayout(V) isa BandedColumns{DenseColumnMajor}
        @test L[2:4,1:3] == BandedMatrix(V) == LowerTriangular(Matrix(A))[2:4,1:3]
        @test L[2:4,1:3] isa BandedMatrix
        @test bandwidths(L[2:4,1:3]) == bandwidths(V) == (1,1)
        @test U[2:4,1:3] == UpperTriangular(Matrix(A))[2:4,1:3]
        @test U[2:4,1:3] isa BandedMatrix
        @test bandwidths(U[2:4,1:3]) == (-1,2)
    end

    @testset "row major" begin
        A = brand(5,5,2,1)
        b = randn(5)
        U = UpperTriangular(A')
        @test MemoryLayout(U) isa TriangularLayout{'U','N',BandedRows{DenseColumnMajor}}
        @test ArrayLayouts.lmul!(U,copy(b)) ≈ U*b
    end
end


text:
@testset "AdditiveWhiteGaussianNoise" begin
    @info "Test: AdditiveWhiteGaussianNoise"

    @testset "API" begin
        @test_throws MethodError AdditiveWhiteGaussianNoise()
        @test_throws ArgumentError AdditiveWhiteGaussianNoise(0.0, -0.1)
        @test AdditiveWhiteGaussianNoise(0.1) == AdditiveWhiteGaussianNoise(0.0, 0.1)

        n = AdditiveWhiteGaussianNoise(0.1)
        A = rand(Gray, 255, 255)
        B1 = similar(A)
        B2 = copy(A)
        out_B1 = apply_noise!(B1, A, n)
        out_B2 = apply_noise!(B2, n)
        # issue #7
        @test out_B2 == B2
        @test out_B1 == B1

        B3 = apply_noise(A, n)

        B1 = similar(A)
        B2 = copy(A)
        apply_noise!(B1, A, n; rng = MersenneTwister(0))
        apply_noise!(B2, n; rng = MersenneTwister(0))
        B3 = apply_noise(A, n; rng = MersenneTwister(0))
        B4 = apply_noise(Float64, A, n; rng = MersenneTwister(0))
        @test B1 == B2 == B3 == B4
    end

    @testset "types" begin
        n = AdditiveWhiteGaussianNoise(0.1)

        # Gray
        type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
        A = [1.0 1.0 1.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
        for T in type_list
            a = T.(A)
            apply_noise(a, n)

            b1 = similar(A, floattype(eltype(A)))
            b2 = copy(A)
            apply_noise!(b1, A, n; rng = MersenneTwister(0))
            apply_noise!(b2, n; rng = MersenneTwister(0))
            b3 = apply_noise(A, n; rng = MersenneTwister(0))
            b4 = apply_noise(floattype(eltype(A)), A, n; rng = MersenneTwister(0))
            @test b1 == b2 == b3 == b4
        end

        # RGB
        type_list = generate_test_types([Float32, N0f8], [RGB])
        A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
            RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
            RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
        for T in type_list
            a = T.(A)
            apply_noise(a, n)

            b1 = similar(A, floattype(eltype(A)))
            b2 = copy(A)
            apply_noise!(b1, A, n; rng = MersenneTwister(0))
            apply_noise!(b2, n; rng = MersenneTwister(0))
            b3 = apply_noise(A, n; rng = MersenneTwister(0))
            b4 = apply_noise(floattype(eltype(A)), A, n; rng = MersenneTwister(0))
            @test b1 == b2 == b3 == b4
        end
    end
end


text:
module CoreNLPClient

using HTTP

export CoreNLP, NLPAnnotation, getNLPAnnotations

# ===========================================================================================
# function parseStringProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the string properties
#                    in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:"property_value".
# output:
#   (property_name, property_value, number_of_lines_parsed)
function parseStringProperty(input::SubString{String})::Tuple{String,String,Int}
    result = ("","",0)
    stringPropertyMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*\"(.|[:])*\"", input)
    if stringPropertyMatch !== nothing
        splits = split(stringPropertyMatch.match, ":")
        key = split(splits[1], "\"")[2]
        value = splits[2]
        for i = 3:length(splits)
            value = value * ":" * splits[i]
        end
        value = split(value, "\"")[2]
        result = (key,value,1)
    end
    result
end

# ===========================================================================================
# function parseFloatProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the float properties
#                    in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:property_value, where the
#          property_value is a float number.
# output:
#   (property_name, property_value, number_of_lines_parsed)
# note:
#   The parsing should use the function before parseIntProperty, otherwise the result will be
#   wrong.
function parseFloatProperty(input::SubString{String})::Tuple{String,Float64,Int}
    result = ("",0,0)
    floatPropertyMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*(\d)*\.(\d)+", input)
    if floatPropertyMatch !== nothing
        key = split(split(floatPropertyMatch.match, ":")[1], "\"")[2]
        value = parse(Float64,split(floatPropertyMatch.match, ":")[2])
        result = (key,value,1)
    end
    result
end

# ===========================================================================================
# function parseIntIntervalProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the int interval
#                    properties in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:property_value, where the
#          property_value is an int interval in form of [i1, i2].
# output:
#   (property_name, property_value, number_of_lines_parsed)
function parseIntIntervalProperty(input::SubString{String})::Tuple{String,Vector{Int},Int}
    result = ("",[0,0],0)
    IntIntervalPropertyMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*\[[ \t]*(-)?(\d)+[ \t]*,[ \t]*(-)?(\d)+[ \t]*\]", input)
    if IntIntervalPropertyMatch !== nothing
        splits = split(IntIntervalPropertyMatch.match, ":")
        key = split(splits[1], "\"")[2]
        values = Int[]
        matches = eachmatch(r"(-)?(\d)+", splits[2])
        for m in matches
            push!(values, parse(Int, m.match))
        end
        @assert length(values) == 2
        result = (key,values,1)
    end
    result
end

# ===========================================================================================
# function parseIntProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the int properties
#                    in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:property_value, where the
#          property_value is an int.
# output:
#   (property_name, property_value, number_of_lines_parsed)
function parseIntProperty(input::SubString{String})::Tuple{String,Int,Int}
    result = ("",0,0)
    intPropertyMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*(\d)+", input)
    if intPropertyMatch !== nothing
        key = split(split(intPropertyMatch.match, ":")[1], "\"")[2]
        value = parse(Int,split(intPropertyMatch.match, ":")[2])
        result = (key,value,1)
    end
    result
end

# ===========================================================================================
# function parseVectorProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the vector properties
#                    in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:property_value, where the
#          property_value is in form of [{item 1}, {item 2}, ..., {item k}].
# output:
#   (property_name, property_value, number_of_lines_parsed)
function parseVectorProperty(input::Vector{SubString{String}})::Tuple{String,Vector{Any},Int}
    result = []
    vectorPropertyMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*\[", input[1])
    if vectorPropertyMatch !== nothing
        myKey = split(split(vectorPropertyMatch.match, ":")[1], "\"")[2]
        n = length(input)
        i = 2
        while i <= n
            key,value,numLinesOfValue = parseDict(input[i:n])
            if numLinesOfValue > 0
                @assert key == ""
                push!(result,value)
                i += numLinesOfValue
            else
                closeMatch = match(r"^[ \t]*\]", input[i])
                i += 1
                if closeMatch !== nothing
                    break
                end
            end
        end

        (myKey,result,i-1)
    else
        ("",[],0)
    end
end

# ===========================================================================================
# function parseDictProperty
# brief description: This is a utility function for the main function CoreNLP(). In CoreNLP(),
#                    annotations of the input text are retrieved from the Stanford CoreNLP
#                    server, and this utility function is used for parsing the dict properties
#                    in the annotations.
# input:
#   input: A substring in the annotations in format of property_name:property_value, where the
#          property_value is in form of {prop1:value1, prop2:value2, ..., propK:valueK}.
# output:
#   (property_name, property_value, number_of_lines_parsed)
function parseDict(input::Vector{SubString{String}})::Tuple{String,Dict{String,Any},Int}
    myKey = ""
    result = Dict{String,Any}()

    openMatch = match(r"^[ \t]*{", input[1])
    if openMatch === nothing
        openMatch = match(r"^[ \t]*\"(.)+\"[ \t]*:[ \t]*{", input[1])
        if openMatch !== nothing
            myKey = split(split(openMatch.match, ":")[1], "\"")[2]
        end
    end
    if openMatch === nothing
        return ("",Dict{String,Any}(), 0)
    end
    #println("Dict $(input[1]) starts:")

    n = length(input)
    i = 2
    closeMatch = nothing
    while i <= n
        #println(input[i])
        key,value,numLinesOfValue = parseStringProperty(input[i])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        key,value,numLinesOfValue = parseFloatProperty(input[i])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        key,value,numLinesOfValue = parseIntProperty(input[i])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        key,value,numLinesOfValue = parseIntIntervalProperty(input[i])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        key,value,numLinesOfValue = parseVectorProperty(input[i:n])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        key,value,numLinesOfValue = parseDict(input[i:n])
        if numLinesOfValue > 0
            result[key] = value
            i += numLinesOfValue
            continue
        end
        closeMatch = match(r"^[ \t]*}", input[i])
        i += 1
        if closeMatch !== nothing
            break
        end
    end
    #println("Dict $(input[1]) ends.")
    @assert closeMatch !== nothing
    (myKey,result,i-1)
end

# ===========================================================================================
# function CoreNLP
# brief description: This is the main function of this package. It is used for querying the
#                    Stanford CoreNLP server for the annotations of an input text.
# input:
#   serverURL: The URL of the CoreNLP server.
#   input:  The input text.
# output:
#   A dictionary of the result retrieved from the CoreNLP server.
# note:
#   A CoreNLP server is needed for using the package. Stanford CoreNLP can be downloaded from:
#   https://stanfordnlp.github.io/CoreNLP/
#   The server must be started whenever the package is used. The detail of starting the server
#   could be found on webpage:
#   https://stanfordnlp.github.io/CoreNLP/corenlp-server.html
function CoreNLP(serverURL::String, input::String)::Dict{String,Any}
    url = serverURL
    if !startswith(url, "http://")
        url = "http://" * url
    end
    splits = split(url, ":")
    @assert length(splits) <= 3
    if length(splits) < 3
        url = url * ":9000"
    end

    res=HTTP.post(url, [], input)
    body = "{"*split(String(res), "\n{")[end]
    annotations = split(body, "\n")
    key,value,numLinesOfValue = parseDict(annotations)
    @assert key == "" && numLinesOfValue > 0
    value
end

# ===========================================================================================
# function CoreNLP
# brief description: This is the shorthand function for the full-argument function above.
function CoreNLP(input::String)::Dict{String,Any}
    CoreNLP("http://localhost:9000", input)
end

# ===========================================================================================
# struct TNLPAnnotation
# brief description: This struct is used as a node type for an annotation graph extracted
#                    using the Stanford CoreNLP toolkit.
# fields:
#   token: A text input to CoreNLP is partitioned into words and punctuations. A token is
#          either a word or a punctuation.
#   POS: Part of Speech. Please visit https://sites.google.com/site/partofspeechhelp/ for
#        more details of POS.
#   governors: CoreNLP analyzes dependencies among tokens. The dependencies are represented
#              with a directed graph. The nodes of this graph are tokens. Every directed edge
#              is pointed from its governor to its dependent. This field is a vector of the
#              indices of the governors. If an index is 0, the governor is ROOT of the
#              dependency tree. Otherwise, the index must be a positive integer to indicate
#              the governor's location in a sentence.
#   dependents: The nodes connected by this node through the dependencies as explained in the
#               description of governors. This field is a vector of the indices of the
#               dependents in a sentence.
#   depLabels: One label for each dependent to describe the type of dependency.
struct NLPAnnotation
    token::String
    POS::String
    governors::Vector{Int}
    dependents::Vector{Int}
    depLabels::Vector{String}
end

# ===========================================================================================
# function getNLPAnnotations
# brief description: Annotate text using Stanford CoreNLP toolkit.
# input:
#   text: The text to annotate.
# output:
#   An annotated document with the input text. The document is represented by a vector of
#   sentences. Each sentence is represented by a vector of TNLPAnnotation.
# note:
#   (1) By default, we use localhost:9000 as the URL of CoreNLP. If a customized URL needs to
#       be set, we could use the "CORENLP_URL" environment variable to specify the URL.
function getNLPAnnotations(text::String)::Vector{Vector{NLPAnnotation}}
    # ---------------------------------------------------------------------------------------
    # step 1: Fetch CoreNLP server URL from environment variables.
    coreNLPUrl = ""
    if "CORENLP_URL" in keys(ENV)
        coreNLPUrl = ENV[CORENLP_URL]
    end

    # ---------------------------------------------------------------------------------------
    # step 2: Get an annotated result from CoreNLP.
    nlpResult = if coreNLPUrl == ""
        CoreNLP(text)
    else
        CoreNLP(coreNLPUrl, text)
    end

    # ---------------------------------------------------------------------------------------
    # step 3: Iterate through the annotated result to fill the output sentence by sentence.
    output = Vector{NLPAnnotation}[]
    for sentence in nlpResult["sentences"]
        # (3.1) create a vector for the annotations of a sentence
        annotationOfSentence = NLPAnnotation[]

        # (3.2) get tokens and dependencies from the CoreNLP annotations of this sentence
        tokens = sentence["tokens"]
        deps = sentence["enhancedPlusPlusDependencies"]

        # (3.3) add the tokens to the vector of annotations
        for token in tokens
            annotationOfToken = NLPAnnotation(token["word"], token["pos"], Int[], Int[], String[])
            push!(annotationOfSentence, annotationOfToken)
        end

        # (3.4) set the dependencies in the vector of annotations
        for dep in deps
            governor = dep["governor"]
            dependent = dep["dependent"]
            label = dep["dep"]
            depNode = annotationOfSentence[dependent]
            push!(depNode.governors, governor)
            if governor > 0
                govNode = annotationOfSentence[governor]
                push!(govNode.dependents, dependent)
                push!(govNode.depLabels, label)
            end
        end

        # (3.5) insert this vector of annotations into the output
        push!(output, annotationOfSentence)
    end

    # ---------------------------------------------------------------------------------------
    # step 4: return the output
    output
end

end # module


text:
#using FixedSizeArrays
using CompScienceMeshes
using Test

p = point(0.0, 1.0)
q = point(2.0, 1.0)
r = point(2.0, 2.0)
a = point(1.0, 0.0)
b = point(2.0, 0.0)
c = point(1.0, 2.0)

@test CompScienceMeshes.leftof(p, a,b) == true
@test CompScienceMeshes.intersectlines(p,r, a,c) == [1.0, 1.5]

α = CompScienceMeshes.sutherlandhodgman2d([p,q,r], [a,b,c])
@test length(α) == 4
@test α[1] == [1.2, 1.6]
@test α[2] == [1.0, 1.5]
@test α[3] == [1.0, 1.0]
@test α[4] == [1.5, 1.0]

β = CompScienceMeshes.sutherlandhodgman([p,q,r], [a,b,c])
@test length(β) == length(α)
for i in eachindex(β) @test β[i] == α[i] end

P = simplex(p, q, r)
Q = simplex(a, b, c)
R = intersection(P,Q)

@test length(R) == 2
@test R[1].vertices == [α[1], α[2], α[3]]
@test R[2].vertices == [α[1], α[3], α[4]]

t = point(1.0, 2.0, 3.0)
p = point(0.0, 1.0, 0.0) + t
q = point(2.0, 1.0, 0.0) + t
r = point(2.0, 2.0, 0.0) + t
a = point(1.0, 0.0, 0.0) + t
b = point(2.0, 0.0, 0.0) + t
c = point(1.0, 2.0, 0.0) + t

P = simplex(p,q,r)
Q = simplex(a,b,c)
S = intersection(P,Q)

# Test whether translation and intersection commute
@test length(S) == 2
for i in 1:2
    for j in 1:3
        for k in 1:2
            @test S[i].vertices[j][k] ≈ R[i].vertices[j][k]+t[k]
        end
        @test S[i].vertices[j][3] ≈ t[3]
    end
end


text:
"""
	$(SIGNATURES)

Abstract type for admissions rules.
Not a ModelObject. This is combined with an `AbstractAdmProbFct` that potentially contains calibrated parameters.
"""
abstract type AbstractAdmissionsRule{I1, F1 <: Real} end

"""
	$(SIGNATURES)

Abstract type for switches from which admissions rules are constructed.
"""
abstract type AbstractAdmissionsSwitches{I1, F1 <: Real} end

Lazy.@forward AbstractAdmissionsRule.switches (
    StructLH.describe, open_admission, min_coll_set_prob, n_colleges
    );

"""
	$(SIGNATURES)

Initialize an admission rule from its switches.
"""
function make_admissions(switches :: AbstractAdmissionsSwitches) end

# """
# 	$(SIGNATURES)

# Variable used to use in admissions rule. Distinct from the ranking variable used in sequential admissions.
# """
# percentile_var(switches :: AbstractAdmissionsSwitches) =
#     switches.pctVar;

n_colleges(a :: AbstractAdmissionsSwitches) = a.nColleges;

# By default, college sets are of the form `1:n`.
n_college_sets(a :: AbstractAdmissionsRule) = n_colleges(a);
college_set(a :: AbstractAdmissionsRule, iSet :: Integer) = 1 : iSet;
open_admission(switches :: AbstractAdmissionsSwitches) = false;

# Iterate over college sets
function Base.iterate(a :: AbstractAdmissionsRule, j)
    if j > n_college_sets(a)
        return nothing
    else
        return college_set(a, j), j+1
    end
end

Base.iterate(a :: AbstractAdmissionsRule) = Base.iterate(a, 1);

# min_coll_set_prob(a :: AbstractAdmissionsRule{I1, F1}) where {I1, F1} =
#     min_coll_set_prob(a.switches);
min_coll_set_prob(a :: AbstractAdmissionsSwitches{I1, F1}) where {I1, F1} =
    a.minCollSetProb;

# Probability of being admitted at each college
function admission_probs(a :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    pctV :: AbstractVector) where {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1}

    nSets = n_college_sets(a);
    nc = n_colleges(a);
    J = length(pctV);
    prob_jcM = zeros(F1, J, nc);
    for iSet = 1 : nSets
        probV = prob_coll_set(a, admProbFct, iSet, pctV);
        iCollV = college_set(a, iSet);
        prob_jcM[:, iCollV] .+= probV;
    end
    @assert check_float_array(prob_jcM, 0.0, 1.0001);
    bracket_array!(prob_jcM, 0.0, 1.0);
    return prob_jcM
end

# Prob of each college set. Vector input. Each type defines this for scalar input.
# Inefficient +++
function prob_coll_set(a :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    iSet :: Integer,
    scoreV :: AbstractVector{F2}) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    J = length(scoreV);
    probV = Vector{F1}(undef, J);
    for j = 1 : J
        probV[j] = prob_coll_set(a, admProbFct, iSet, scoreV[j]);
    end
    return probV
end

# Prob college sets for several types.
function prob_coll_sets(a :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    scoreV :: AbstractVector{F2}) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    nSets = n_college_sets(a);
    n = length(scoreV);
    prob_jsM = Matrix{F1}(undef, n, nSets);
    for j = 1 : n
        prob_jsM[j, :] = prob_coll_sets(a, admProbFct, hsGpaPctV[j]);
    end

    @assert all_at_least(prob_jsM, a.switches.minCollSetProb)
    @assert all_at_most(prob_jsM, 1.0)
    @check all(isapprox.(sum(prob_jsM, dims = 2), 1.0, atol = 1e-6))
    return prob_jsM
end

# ----------


text:
# Sane behavior when run from the REPL
source_dir = typeof(Base.source_dir()) == Void ? joinpath(Pkg.dir("NeuralNet"), "src") : Base.source_dir()

"""
align(algorithm::Compat.String, gap_penalty::Float64, gap_cont_penalty::Float64,
            seq1::Compat.String, seq2::Compat.String, smat::Compat.String, outputfile::Compat.String)
Align SEQ1 to SEQ2 using the scoring matrix SMAT and paramters GAP_PENALTY and GAP_CONT_PENALTY
for the ALGORITHM defined by user.
"""
function simple_encoder(N::Int, alpha::Float64, outputfile::Compat.String)

    include(joinpath(source_dir, "algorithms", "nn.jl"));
    # run
    Base.invokelatest(encode_8x3x8, N, alpha, outputfile);

end

function train_model(N::Int, alpha::Float64, hl_size::Int, write_output::Bool,
            balance::Int, C::Int, pos_seqs::Compat.String, total_seq::Compat.String, outputfile::Compat.String)

    include(joinpath(source_dir, "algorithms", "nn.jl"));
    include(joinpath(source_dir, "utils", "parse.jl"));

    tdata, labels = Base.invokelatest(parse_training, pos_seqs, total_seq, balance);

    Base.invokelatest(train_nn, tdata, labels, hl_size, N, alpha,
                        outputfile, write_output, C);

end

function predict(N::Int, alpha::Float64, hl_size::Int, write_output::Bool,
            balance::Int, C::Int, pos_seqs::Compat.String, total_seq::Compat.String,
            predict_seq::Compat.String, outputfile::Compat.String)

    include(joinpath(source_dir, "algorithms", "nn.jl"));
    include(joinpath(source_dir, "utils", "parse.jl"));

    tdata, labels = Base.invokelatest(parse_training, pos_seqs, total_seq, balance);
    test_data = Base.invokelatest(parse_testing, predict_seq)

    Base.invokelatest(nn_predict_on_data, tdata, labels, test_data, hl_size, N, alpha,
                        outputfile, write_output, C);

end


text:
# We roughly follow DLMF notation
# Q[x,n] = k[n] * x^(n-1)
# Note that we have
# Q[x,n] = (1/γ[n] - β[n-1]/γ[n]) * x * Q[x,n-1] - γ[n-1]/γ[n] * Q[x,n]
#
#

function lanczos!(Ns, X::AbstractMatrix{T}, W::AbstractMatrix{T}, γ::AbstractVector{T}, β::AbstractVector{T}, R::AbstractMatrix{T}) where T
    for n = Ns
        v = view(R,:,n);
        p1 = view(R,:,n-1);
        muladd!(one(T), X, p1, zero(T), v); # TODO: `mul!(v, X, p1)`
        β[n-1] = -dot(v,W,p1)
        axpy!(β[n-1],p1,v);
        if n > 2
            p0 = view(R,:,n-2)
            axpy!(-γ[n-1],p0,v)
        end
        γ[n] = sqrt(dot(v,W,v));
        lmul!(inv(γ[n]), v)
    end
    γ,β,R
end

const PaddedVector{T} = CachedVector{T,Vector{T},Zeros{T,1,Tuple{OneToInf{Int}}}}
const PaddedMatrix{T} = CachedMatrix{T,Matrix{T},Zeros{T,2,NTuple{2,OneToInf{Int}}}}

mutable struct LanczosData{T}
    X
    W
    γ::PaddedVector{T}
    β::PaddedVector{T}
    R::UpperTriangular{T,PaddedMatrix{T}}
    ncols::Int

    function LanczosData{T}(X, W, γ, β, R) where {T}
        R[1,1] = 1;
        p0 = view(R,:,1);
        γ[1] = sqrt(dot(p0,W,p0))
        lmul!(inv(γ[1]), p0)
        new{T}(X, W, γ, β, R, 1)
    end
end

LanczosData(X, W, γ::AbstractVector{T}, β, R) where {T} = LanczosData{T}(X, W, γ, β, R)
LanczosData(X::AbstractMatrix{T}, W::AbstractMatrix{T}) where T = LanczosData(X, W, zeros(T,∞), zeros(T,∞), UpperTriangular(zeros(T,∞,∞)))

function LanczosData(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
    x = axes(P,1)
    wP = weighted(P)
    X = jacobimatrix(P)
    W = wP \ (w .* P)
    LanczosData(X, W)
end

function resizedata!(L::LanczosData, N)
    N ≤ L.ncols && return L
    resizedata!(L.R, N, N)
    resizedata!(L.γ, N)
    resizedata!(L.β, N)
    resizedata!(L.W, N, N)
    lanczos!(L.ncols+1:N, L.X, L.W, L.γ, L.β, L.R)
    L.ncols = N
    L
end

struct LanczosConversion{T} <: LayoutMatrix{T}
    data::LanczosData{T}
end

size(::LanczosConversion) = (∞,∞)
copy(R::LanczosConversion) = R

Base.permutedims(R::LanczosConversion{<:Number}) = transpose(R)

bandwidths(::LanczosConversion) = (0,∞)
colsupport(L::LanczosConversion, j) = 1:maximum(j)
rowsupport(L::LanczosConversion, j) = minimum(j):∞

function _lanczosconversion_getindex(R, k, j)
    resizedata!(R.data, max(maximum(k), maximum(j)))
    R.data.R[k,j]
end

getindex(R::LanczosConversion, k::Integer, j::Integer) = _lanczosconversion_getindex(R, k, j)
getindex(R::LanczosConversion, k::AbstractUnitRange, j::AbstractUnitRange) = _lanczosconversion_getindex(R, k, j)

inv(R::LanczosConversion) = ApplyArray(inv, R)

Base.BroadcastStyle(::Type{<:LanczosConversion}) = LazyArrays.LazyArrayStyle{2}()

struct LanczosConversionLayout <: AbstractLazyLayout end

LazyArrays.simplifiable(::Mul{LanczosConversionLayout,<:PaddedLayout}) = Val(true)
function copy(M::Mul{LanczosConversionLayout,<:PaddedLayout})
    resizedata!(M.A.data, maximum(colsupport(M.B)))
    M.A.data.R * M.B
end

function copy(M::Ldiv{LanczosConversionLayout,<:PaddedLayout})
    resizedata!(M.A.data, maximum(colsupport(M.B)))
    M.A.data.R \ M.B
end

function getindex(L::Ldiv{LanczosConversionLayout,<:AbstractBandedLayout}, ::Colon, j::Integer)
    m = maximum(colrange(L.B,j))
    [L.A[1:m,1:m] \ L.B[1:m,j]; Zeros{eltype(L)}(∞)]
end

MemoryLayout(::Type{<:LanczosConversion}) = LanczosConversionLayout()
triangulardata(R::LanczosConversion) = R
sublayout(::LanczosConversionLayout, ::Type{<:Tuple{KR,Integer}}) where KR =
    sublayout(PaddedLayout{UnknownLayout}(), Tuple{KR})

function sub_paddeddata(::LanczosConversionLayout, S::SubArray{<:Any,1,<:AbstractMatrix,<:Tuple{Any,Integer}})
    P = parent(S)
    (kr,j) = parentindices(S)
    resizedata!(P.data, j)
    paddeddata(view(UpperTriangular(P.data.R.data.data), kr, j))
end

# struct LanczosJacobiMatrix{T} <: AbstractBandedMatrix{T}
#     data::LanczosData{T}
# end

struct LanczosJacobiBand{T} <: LazyVector{T}
    data::LanczosData{T}
    diag::Symbol
end

size(P::LanczosJacobiBand) = (∞,)
resizedata!(A::LanczosJacobiBand, n) = resizedata!(A.data, n)

# γ[n+1]*Q[x,n+1] + β[n]*Q[x,n] + γ[n-1]*Q[x,n-1] = x  * Q[x,n]
function _lanczos_getindex(C::LanczosJacobiBand, I)
    resizedata!(C, maximum(I)+1)
    if C.diag == :du
        C.data.γ.data[I .+ 1]
    else # :d
        -C.data.β.data[I]
    end
end

getindex(A::LanczosJacobiBand, I::Integer) = _lanczos_getindex(A, I)
getindex(A::LanczosJacobiBand, I::AbstractVector) = _lanczos_getindex(A, I)
getindex(K::LanczosJacobiBand, k::AbstractInfUnitRange{<:Integer}) = view(K, k)
getindex(K::SubArray{<:Any,1,<:LanczosJacobiBand}, k::AbstractInfUnitRange{<:Integer}) = view(K, k)

copy(A::LanczosJacobiBand) = A # immutable

struct LanczosPolynomial{T,Weight,Basis} <: OrthogonalPolynomial{T}
    w::Weight # Weight of orthogonality
    P::Basis # Basis we use to represent the OPs
    data::LanczosData{T}
end

function LanczosPolynomial(w_in::AbstractQuasiVector, P::AbstractQuasiMatrix)
    Q = normalize(P)
    wQ = weighted(Q)
    w = wQ * (wQ \ w_in) # expand weight in basis
    LanczosPolynomial(w, Q, LanczosData(w, Q))
end

LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, orthonormalpolynomial(singularities(w)))

==(A::LanczosPolynomial, B::LanczosPolynomial) = A.w == B.w
==(::LanczosPolynomial, ::OrthogonalPolynomial) = false # TODO: fix
==(::OrthogonalPolynomial, ::LanczosPolynomial) = false # TODO: fix
==(::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, ::LanczosPolynomial) = false # TODO: fix
==(::LanczosPolynomial, ::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = false # TODO: fix

normalize(Q::LanczosPolynomial) = Q
normalize(Q::AbstractQuasiMatrix) = Normalized(Q)

OrthogonalPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w)
orthogonalpolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w)
orthogonalpolynomial(w::SubQuasiArray) = orthogonalpolynomial(parent(w))[parentindices(w)[1],:]
orthonormalpolynomial(w::AbstractQuasiVector) = normalize(orthogonalpolynomial(w))

orthogonalityweight(Q::LanczosPolynomial) = Q.w

axes(Q::LanczosPolynomial) = (axes(Q.w,1),OneToInf())

_p0(Q::LanczosPolynomial) = inv(Q.data.γ[1])*_p0(Q.P)

jacobimatrix(Q::LanczosPolynomial) = SymTridiagonal(LanczosJacobiBand(Q.data, :d), LanczosJacobiBand(Q.data, :du))

Base.summary(io::IO, Q::LanczosPolynomial{T}) where T = print(io, "LanczosPolynomial{$T} with weight with singularities $(singularities(Q.w))")
Base.show(io::IO, Q::LanczosPolynomial{T}) where T = summary(io, Q)

Base.array_summary(io::IO, C::SymTridiagonal{T,<:LanczosJacobiBand}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
    print(io, Base.dims2string(length.(inds)), " SymTridiagonal{$T} Jacobi operator from Lanczos")

Base.array_summary(io::IO, C::LanczosConversion{T}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
    print(io, Base.dims2string(length.(inds)), " LanczosConversion{$T}")

Base.array_summary(io::IO, C::LanczosJacobiBand{T}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
    print(io, Base.dims2string(length.(inds)), " LanczosJacobiBand{$T}")

# Sometimes we want to expand out, sometimes we don't

QuasiArrays.ApplyQuasiArray(Q::LanczosPolynomial) = ApplyQuasiArray(*, arguments(ApplyLayout{typeof(*)}(), Q)...)

function \(A::LanczosPolynomial{T}, B::LanczosPolynomial{V}) where {T,V}
    A == B && return Eye{promote_type(T,V)}(∞)
    inv(LanczosConversion(A.data)) * (A.P \ B.P)  * LanczosConversion(B.data)
end
\(A::OrthogonalPolynomial, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
\(A::Normalized, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)
\(Q::LanczosPolynomial, A::OrthogonalPolynomial) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
\(Q::LanczosPolynomial, A::Normalized) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
\(Q::LanczosPolynomial, A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = inv(LanczosConversion(Q.data)) * (Q.P \ A)
\(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, Q::LanczosPolynomial) = (A \ Q.P) * LanczosConversion(Q.data)

ArrayLayouts.mul(Q::LanczosPolynomial, C::AbstractArray) = ApplyQuasiArray(*, Q, C)
function ldiv(Qn::SubQuasiArray{<:Any,2,<:LanczosPolynomial,<:Tuple{<:Inclusion,<:Any}}, C::AbstractQuasiArray)
    _,jr = parentindices(Qn)
    Q = parent(Qn)
    LanczosConversion(Q.data)[jr,jr] \ (Q.P[:,jr] \ C)
end

struct LanczosLayout <: AbstractOPLayout end

MemoryLayout(::Type{<:LanczosPolynomial}) = LanczosLayout()
arguments(::ApplyLayout{typeof(*)}, Q::LanczosPolynomial) = Q.P, LanczosConversion(Q.data)
copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:BroadcastLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{ApplyLayout{typeof(*)},BroadcastLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},BroadcastLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{LanczosLayout,<:ExpansionLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))
LazyArrays._mul_arguments(Q::LanczosPolynomial) = arguments(ApplyLayout{typeof(*)}(), Q)
LazyArrays._mul_arguments(Q::QuasiAdjoint{<:Any,<:LanczosPolynomial}) = arguments(ApplyLayout{typeof(*)}(), Q)

broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P.P, Q.P)
broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q) = broadcastbasis(+, P.P, Q)
broadcastbasis(::typeof(+), P, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P, Q.P)


text:
using DataFrames

f =  readtable("hightemp.txt", separator = '\t', header = false)
sort!(f, cols = :x3, rev = true)
println(f)


text:
# Gets the robot config for the specified robot
# This contains params like camera matrices etc.

"""
    $(SIGNATURES)

Plugin callback for fetching and caching a robot config for use by other plugins.
"""
function getRobotConfiguration(vis::MeshCat.Visualizer,
                        params::Dict{Symbol, Any},
                        rose_fgl::Tuple{<:AbstractString, <:AbstractString} )
    if !haskey(params, :robotConfig)
		@info "Getting robot config for robot $(getGraffConfig().robotId)..."
		params[:robotConfig] = getRobotConfig()
	end
    return nothing
end

"""
    $(SIGNATURES)

Plugin callback for fetching robot config in the case of a local graph.
User must specify params[:robotConfig] beforehand in this case.
"""
function getRobotConfiguration(vis::MeshCat.Visualizer,
                        params::Dict{Symbol, Any},
                        rose_fgl::FactorGraph )
    #
	if !haskey(params, :robotConfig)
		error("For local graphs, please manually specify your robot config as a parameter params[:robotConfig]...")
	end
    return nothing
end