text:
"""
    Fermi.Geometry

Module handling molecule and atoms data

# exports:

    Atom      Object storing information about an atom
    Molecule  Object storing information about a molecule (group of atoms)
"""
module Geometry

export Molecule
export Atom

using Fermi
using Fermi.Output
using Fermi.PhysicalConstants: atomic_number, bohr_to_angstrom, angstrom_to_bohr
using LinearAlgebra
using Lints
using Formatting

struct InvalidMolecule <: Exception
    msg::String
end
Base.showerror(io::IO, e::InvalidMolecule) = print(io, "InvalidMolecule: ", e.msg)

"""
    Fermi.Atom

Object storing information about an atom.

# Fields:

    AtomicSymbol    Atomic symbol
    Z               Atomic number
    xyz             xyz array in Angstrom
"""
struct Atom
    AtomicSymbol::String
    Z::Int
    xyz::Tuple{Float64,Float64,Float64}
end

"""
    Fermi.Molecule

Object storing information about a molecule (group of atoms).

# Fields:

    atoms        Array with Fermi.Atom objects
    charge       Charge of the molecule
    multiplicity Multiplicity ``(2Ms + 1)``
    Nα           Number of Alpha electrons
    Nβ           Number of Beta electrons
    Vnuc         Nuclear repulsion energy
"""
struct Molecule
    atoms::Tuple
    charge::Int
    multiplicity::Int
    Nα::Int
    Nβ::Int
    Vnuc::Float64
end

"""
    Fermi.Molecule()

Uses data stored in Fermi.CurrentOptions to return a Molecule object.

"""
function Molecule()
    molstring = Fermi.CurrentOptions["molstring"]
    unit = Fermi.CurrentOptions["unit"]
    Molecule(molstring, unit)
end

"""
    Fermi.Molecule(molstring::String, unit::String="angstrom")

Uses data stored in Fermi.CurrentOptions for charge and multiplicity to build
a Molecule object given XYZ string.
"""
function Molecule(molstring::String, unit::String="angstrom")
    charge = Fermi.CurrentOptions["charge"]
    multiplicity = Fermi.CurrentOptions["multiplicity"]
    Molecule(molstring, charge, multiplicity, unit)
end

"""
    Fermi.Molecule(molstring::String, charge::Int, multiplicity::Int, unit::String="angstrom")

Produces a Molecule object from the XYZ string.
"""
function Molecule(molstring::String, charge::Int, multiplicity::Int, unit::String="angstrom")
    atoms = get_atoms(molstring, unit=unit)
    Molecule(atoms, charge, multiplicity)
end

"""
    Fermi.Molecule(molstring::String, charge::Int, multiplicity::Int, unit::String="angstrom")

Produces a Molecule object from an array of Atom objects.
"""
function Molecule(atoms::Array{Atom,1}, charge::Int, multiplicity::Int)

    # Compute Nuclear repulsion
    Vnuc = 0.0
    for i in eachindex(atoms)
        for j in 1:(i-1)
            @inbounds Vnuc += nuclear_repulsion(atoms[i], atoms[j])
        end
    end

    # Compute number of electrons

    nelec = 0
    for i in eachindex(atoms)
        @inbounds nelec += atoms[i].Z
    end

    nelec -= charge

    if nelec ≤ 0
        throw(Fermi.InvalidMolecule("Invalid charge ($charge) for given molecule"))
    end

    αexcess = multiplicity-1

    if isodd(nelec) != isodd(αexcess)
        throw(Fermi.InvalidMolecule("Incompatible charge $(charge) and multiplicity $(multiplicity)"))
    end

    Nβ = (nelec - αexcess)/2
    Nα = nelec - Nβ

    out =  Molecule(Tuple(atoms), charge, multiplicity, Nα, Nβ, Vnuc)

    return out
end

function to_lints_molecule(M::Molecule)
    natoms = length(M.atoms)
    zs = zeros(Int64,natoms)
    pos = fill(Float64[],natoms)
    for i=1:natoms
        zs[i] = M.atoms[i].Z
        pos[i] = collect(M.atoms[i].xyz)
    end
    Lints.Molecule(zs,pos)
end

function print_out(M::Molecule)
    @output "   • Molecule:\n\n"
    output(get_xyz(M))
    @output "\n"
    @output "\nCharge: {}   " out.charge
    @output "Multiplicity: {}   \n" out.multiplicity
    @output "Nuclear repulsion: {:15.10f}\n\n" out.Vnuc
end
"""
    Fermi.Geometry.nuclear_repulsion(A::Atom, B::Atom)

Returns the repulsion energy between the two given atoms.
"""
function nuclear_repulsion(A::Atom, B::Atom)
    return (A.Z*B.Z)/(√((A.xyz.-B.xyz)⋅(A.xyz.-B.xyz))*angstrom_to_bohr)
end

"""
    Fermi.get_atoms(molstring::String; unit::String="angstrom")

From a XYZ string, produces an array of Atom objects.
"""
function get_atoms(molstring::String; unit::String="angstrom")

    if unit == "bohr"
        conv = bohr_to_angstrom
    elseif unit == "angstrom"
        conv = 1.0
    else
        error("Invalid unit given to Fermi.Molecule.get_atoms: $unit")
    end

    atoms = Atom[]
    #exp = r"(\w{1,2})\s+([-+]??\d+\.\d+\s+[-+]??\d+\.\d+\s+[-+]??\d+\.\d+)"
    for line in split(strip(molstring), "\n")
        #m = match(exp, line)
        m = split(line)
        if length(m) != 4
            throw(InvalidMolecule("4 columns expected on XYZ string. Got $(length(m))"))
        end

        if m != nothing
            # Convert SubString to String
            AtomicSymbol = String(m[1])
            Z = atomic_number(AtomicSymbol)

            # Convert Substring to String and then String to Float64 Array
            xyz = [0.0, 0.0, 0.0]
            try
                xyz .= parse.(Float64, m[2:4]).*conv
            catch ArgumentError
                throw(InvalidMolecule("Failed to process XYZ string"))
            end

            push!(atoms, Atom(AtomicSymbol, Z, Tuple(xyz)))
        end
    end
    return atoms
end

"""
    Fermi.get_xyz(M::Molecule)

Returns a XYZ string in angstrom for the given Molecule.
"""
function get_xyz(M::Molecule)

    molstring = ""

    for i in eachindex(M.atoms)

        A = M.atoms[i]
        molstring *= format("{}   {: 15.12f}   {: 15.12f}   {: 15.12f}\n", A.AtomicSymbol, A.xyz...)
    end
    return molstring
end

end #Module


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.001 # almost isothermal when gamma reaches 1
equations = CompressibleEulerEquations2D(gamma)

# This is a hand made colliding flow setup without reference. Features Mach=70 inflow from both
# sides, with relative low temperature, such that pressure keeps relatively small
# Computed with gamma close to 1, to simulate isothermal gas
function initial_condition_colliding_flow_astro(x, t, equations::CompressibleEulerEquations2D)
  # change discontinuity to tanh
  # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF)
  # domain size is [-64,+64]^2
  @unpack gamma = equations
  # the quantities are chosen such, that they are as close as possible to the astro examples
  # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...)
  rho = 0.0247
  c = 0.2
  p = c^2 / gamma * rho
  vel = 13.907432274789372
  slope = 1.0
  v1 = -vel*tanh(slope * x[1])
  # add small initial disturbance to the field, but only close to the interface
  if abs(x[1]) < 10
    v1 = v1 * (1 + 0.01 * sin(pi * x[2]))
  end
  v2 = 0.0
  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_colliding_flow_astro

boundary_conditions = (
                       x_neg=BoundaryConditionDirichlet(initial_condition_colliding_flow_astro),
                       x_pos=BoundaryConditionDirichlet(initial_condition_colliding_flow_astro),
                       y_neg=boundary_condition_periodic,
                       y_pos=boundary_condition_periodic,
                      )

surface_flux = flux_lax_friedrichs # HLLC needs more shock capturing (alpha_max)
volume_flux  = flux_ranocha # works with Chandrashekar flux as well
polydeg = 3
basis = LobattoLegendreBasis(polydeg)

# shock capturing necessary for this tough example, however alpha_max = 0.5 is fine
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.0001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-64.0, -64.0)
coordinates_max = ( 64.0,  64.0)

# only refinment in a patch. Needs x=-17/+17 to trigger refinment due to coarse base mesh
refinement_patches = (
  (type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)),
  (type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)),
  (type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)),
  (type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)),
  #(type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)), # very high resolution, takes about 1000s on 2 cores
)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                refinement_patches=refinement_patches,
                periodicity=(false,true),
                n_cells_max=100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 25.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution)

# positivity limiter necessary for this tough example
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
                                                     variables=(Trixi.density, pressure))

###############################################################################
# run the simulation
# use adaptive time stepping based on error estimates, time step roughly dt = 5e-3
sol = solve(ode, SSPRK43(stage_limiter!),
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

ext: jl


text:
using SpecialFunctions:loggamma
using LinearAlgebra:dot

struct PolynomialLimbDark{T,VT <: AbstractVector{T},MT <: AbstractMatrix{T},AT <: AbstractArray{T}} <: AbstractLimbDark
    n_max::Int
    u_n::VT
    g_n::VT
    Mn_coeff::AT
    Nn_coeff::MT
    norm::T
    Mn::VT
    Nn::VT
end

@doc raw"""
    PolynomialLimbDark(u::AbstractVector)

Polynomial limb darkening using analytical integrals. The length of the `u` vector is equivalent to the order of polynomial used; e.g., `[0.2, 0.3]` corresponds to quadratic limb darkening.

# Mathematical form
```math
I(\mu) \propto 1 - u_1(1-\mu) - u_2(1-\mu)^2 - \dots - u_N(1-\mu)^N
```
which is equivalent to the series
```math
I(\mu) \propto -\sum_{i=0}^N{u_i(1-\mu)^i}
```
with the definition $u_0 \equiv -1$.

# Examples
```jldoctest poly
u = [0.4, 0.26] # quadratic and below is 100% analytical
ld = PolynomialLimbDark(u)
ld(0.1, 0.01)

# output
0.9998787880717668
```
```jldoctest poly
u2 = vcat(u, ones(12) ./ 12)
ld2 = PolynomialLimbDark(u2)
ld2(0.1, 0.01)

# output
0.9998740059086433
```

# References

> [Agol, Luger, Foreman-Mackey (2020)](https://ui.adsabs.harvard.edu/abs/2020AJ....159..123A)
>
>   "Analytic Planetary Transit Light Curves and Derivatives for Stars with Polynomial Limb Darkening"

> [Luger et al. (2019)](https://ui.adsabs.harvard.edu/abs/2019AJ....157...64L)
>
>   "starry: Analytic Occultation Light Curves"
"""
function PolynomialLimbDark(u::AbstractVector{S}; maxiter=100) where S
    T = float(S)
    # add constant u_0 term
    n_max = length(u)
    u_n = pushfirst!(copy(u), -one(S))
    # get Green's basis coefficients
    g_n = compute_gn(u_n)
    # compute series expansion for M_n and N_n
    Mn_coeff = compute_Mn_coeffs(T, n_max; maxiter=maxiter)
    Nn_coeff = compute_Nn_coeffs(T, n_max; maxiter=maxiter)

    # calculate flux normalization factor, which only depends on first two terms
    if length(g_n) > 1
        norm = inv(π * (g_n[1] + 2 / 3 * g_n[2]))
    else
        norm = inv(π * g_n[1])
    end

    # pre-allocate temp arrays
    Mn = similar(g_n)
    Nn = similar(g_n)

    return PolynomialLimbDark(n_max, u_n, g_n, Mn_coeff, Nn_coeff, norm, Mn, Nn)
end

function compute(ld::PolynomialLimbDark, b::S, r) where S
    T = float(S)

    ## check for trivial cases
    if b ≥ 1 + r || iszero(r)
        # completely unobscured
        return one(T)
    elseif r ≥ 1 + b
        # completely obscured
        return zero(T)
    end

    r2 = r^2
    b2 = b^2

    if iszero(b)
        onemr2 = 1 - r2
        sqrt1mr2 = sqrt(onemr2)

        # annular ellipse
        if ld.n_max == 0
            flux = ld.g_n[begin] * onemr2
        else
            flux = ld.g_n[begin] * onemr2 + 2 / 3 * ld.g_n[begin + 1] * sqrt1mr2^3
        end
        fac = 2 * r2 * onemr2
        for i in 2:ld.n_max
            flux -= ld.g_n[begin + i] * fac
            fac *= sqrt1mr2
        end
        return flux * π * ld.norm
    end

    # take a moment to calculate values used repeatedly in the
    # numerical routines ahead
    onembmr2 = (r - b + 1) * (1 - r + b)
    onembmr2inv = inv(onembmr2)
    sqonembmr2 = sqrt(onembmr2)
    br = b * r
    fourbr = 4 * br
    fourbrinv = inv(fourbr)
    sqbr = sqrt(br)
    sqbrinv = inv(sqbr)
    onembpr2 = (1 - r - b) * (1 + b + r)
    sqarea = sqarea_triangle(one(T), r, b)
    k2 = max(zero(T), onembmr2 * fourbrinv)
    k = sqrt(k2)
    omemr2 = (1 - r) * (1 + r)
    onemr2mb2 = omemr2 - b2
    onemr2pb2 = omemr2 + b2
    if k2 ≥ 1
        if k2 > 2
            kc2 = 1 - inv(k2)
        else
            kc2 = onembpr2 * onembmr2inv
        end
    else
        if k2 > 0.5
            kc2 = (r - 1 + b) * (b + r + 1) * fourbrinv
        else
            kc2 = 1 - k2
        end
    end
    kc = sqrt(abs(kc2))

    ## Compute uniform term
    s0, kap0, kite_area2, kck = compute_uniform(b, r; b2, r2, sqarea, fourbrinv)

    flux = ld.g_n[begin] * s0

    if ld.n_max == 0
        return flux * ld.norm
    end

    ## compute linear term
    s1, Eofk, Em1mKdm = compute_linear(b, r; k2, kc, kc2, r2, b2, br, fourbr, sqbr, onembmr2, onembpr2, onembmr2inv)

    flux += ld.g_n[begin + 1] * s1

    if ld.n_max == 1
        return flux * ld.norm
    end

    ## calculate quadratic term
    s2 = compute_quadratic(b, r; s0, r2, b2, kap0, kite_area2, k2)
    flux += ld.g_n[begin + 2] * s2
    if ld.n_max == 2
        return flux * ld.norm
    end

    ## higher order (N > 3) terms require solving Mn and Nn integrals
    if k2 < 0.5 && ld.n_max > 3
        downwardM!(ld.Mn, ld.Mn_coeff; n_max=ld.n_max, sqonembmr2, sqbr, onemr2mb2, k, k2, sqarea, kap0, Eofk, Em1mKdm, kite_area2)
    else
        upwardM!(ld.Mn; sqbr, n_max=ld.n_max, sqonembmr2, onemr2mb2, sqarea, k2, kap0, Eofk, Em1mKdm, kite_area2)
    end

    # compute remaining terms
    for n in 3:ld.n_max
        pofgn_M = 2 * r2 * ld.Mn[begin + n] - n / (n + 2) *
                       (onemr2mb2 * ld.Mn[begin + n] + sqarea * ld.Mn[begin + n - 2])
        flux -= ld.g_n[begin + n] * pofgn_M
    end

    return flux * ld.norm
end

###
###   Type methods
###

function Base.show(io::IO, ld::PolynomialLimbDark)
    print(io, "PolynomialLimbDark(", ld.u_n, ")")
end
function Base.show(io::IO, ::MIME"text/plain", ld::PolynomialLimbDark)
    print(io, "PolynomialLimbDark\n u_n: ", ld.u_n)
end
function Base.show(io::IO, ::MIME"text/html", ld::PolynomialLimbDark)
    print(io, "PolynomialLimbDark\n un: ", ld.u_n)
end

###
###   Initialization helpers
###

"""
Transform the `u_n` coefficients to `g_n`, which are coefficients in Green's basis from REF.
"""
compute_gn(u::AbstractVector) = compute_gn!(zero(u), u)

function compute_gn!(g_n::AbstractVector{T}, u_n) where T
    n = length(g_n) - 1
    ## First: calculate the a_n terms (eqn. 10)
    a_n = zero(g_n)
    # constant term
    a_n[begin] = one(T)
    for i in 1:n
        bcoeff = one(T) # binomial coefficient
        for j in 0:i
            a_n[begin + j] -= (-1)^j * bcoeff * u_n[begin + i]
            bcoeff *= (i - j) / (j + 1)
        end
    end

    ## Second: calculate the g_n terms (eqn. 15)
    for j in n:-1:2
        if j ≥ n - 1
            g_n[begin + j] = a_n[begin + j] / (j + 2)
        else
            g_n[begin + j] = a_n[begin + j] / (j + 2) + g_n[begin + j + 2]
        end
    end
    if n ≥ 3
        g_n[begin + 1] = a_n[begin + 1] + 3 * g_n[begin + 3]
    elseif n ≥ 1
        g_n[begin + 1] = a_n[begin + 1]
    end
    if n ≥ 2
        g_n[begin] = a_n[begin] + 2 * g_n[begin + 2]
    else
        g_n[begin] = a_n[begin]
    end

    return g_n
end

###
###   Computation helpers
###

function sort3(p0, p1, p2)
    if p0 < p1
        p0, p1 = p1, p0
    end
    if p1 < p2
        p1, p2 = p2, p1
    end
    if p0 < p1
        p0, p1 = p1, p0
    end
    return p0, p1, p2
end

function sqarea_triangle(p0, p1, p2)
    p0, p1, p2 = sort3(p0, p1, p2)
    sqarea = (p0 + (p1 + p2)) * (p2 - (p0 - p1)) *
             (p2 + (p0 - p1)) * (p0 + (p1 - p2))
    return sqarea
end

function compute_uniform(b::T, r; sqarea, r2=r^2, b2=b^2, fourbrinv=inv(4 * b * r)) where T
    if b ≤ 1 - r
        flux = π * (1 - r2)
        kap0 = convert(T, π)
        kck = zero(T)
        kite_area2 = zero(T)
    else
        kite_area2 = sqrt(sqarea)
        r2m1 = (r - 1) * (r + 1)
        kap0 = atan(kite_area2, r2m1 + b2)
        Πmkap1 = atan(kite_area2, r2m1 - b2)
        kck = kite_area2 * fourbrinv
        flux = Πmkap1 - r2 * kap0 + 0.5 * kite_area2
    end

    return flux, kap0, kite_area2, kck
end

function compute_linear(b::T, r; k2, kc, kc2, onembmr2, onembpr2, onembmr2inv, r2=r^2, b2=b^2, br=b * r, fourbr=4 * br, sqbr=sqrt(br)) where T
    if iszero(b) # case 10
        Λ1 = -2 * π * sqrt(1 - r2)^3
        Eofk = 0.5 * π
        Em1mKdm = 0.25 * π
    elseif b == r
        if r == 0.5 # case 6
            Λ1 = π - 4 / 3
            Eofk = one(T)
            Em1mKdm = one(T)
        elseif r < 0.5 # case 5
            m = 4 * r2
            Eofk = cel(m, one(T), one(T), 1 - m)
            Em1mKdm = cel(m, one(T), one(T), zero(T))
            Λ1 = π + 2 / 3 * ((2 * m - 3) * Eofk - m * Em1mKdm) +
                 (b - r) * 4 * r * (Eofk - 2 * Em1mKdm)
        else # case 7
            m = 4 * r2
            minv = inv(m)
            Eofk = cel(minv, one(T), one(T), 1 - minv)
            Em1mKdm = cel(minv, one(T), one(T), zero(T))
            Λ1 = π + 1 / 3 * ((2 * m - 3) * Em1mKdm - m * Eofk) / r -
                 (b - r) * 2 * (2 * Eofk - Em1mKdm)
        end
    else
        if b + r > 1 # case 2, case 8
            Πofk, Eofk, Em1mKdm = cel(k2, kc, (b - r)^2 * kc2, zero(T), one(T), one(T), 3 * kc2 * (b - r) * (b + r), kc2, zero(T))
            Λ1 = onembmr2 * (Πofk + (-3 + 6 * r2 + 2 * br) * Em1mKdm - fourbr * Eofk) / (3 * sqbr)
        elseif b + r < 1 # case 3 case 9
            bmrdbpr = (b - r) / (b + r)
            μ = 3 * bmrdbpr * onembmr2inv
            p = bmrdbpr^2 * onembpr2 * onembmr2inv
            Πofk, Eofk, Em1mKdm = cel(inv(k2), kc, p, 1 + μ, one(T), one(T), p + μ, kc2, zero(T))
            Λ1 = 2 * sqrt(onembmr2) * (onembpr2 * Πofk - (4 - 7 * r2 - b2) * Eofk) / 3
        else # case
            Λ1 = 2 * acos(1 - 2 * r) - 2 * π * (r > 0.5) - (4 / 3 * (3 + 2 * r - 8 * r2) + 8 * (r + b - 1) * r) * sqrt(r * (1 - r))
            Eofk = one(T)
            Em1mKdm = one(T)
        end
    end

    flux = ((1 - T(r > b)) * 2π - Λ1) / 3

    return flux, Eofk, Em1mKdm
end

function compute_quadratic(b, r; s0, r2, b2, kap0, kite_area2, k2)
    r2pb2 = r2 + b2
    η2 = r2 * (r2pb2 + b2)
    if k2 > 1
        four_pi_eta = 2 * π * (η2 - 1)
    else
        Πmkap1 = atan(kite_area2, (r - 1) * (r + 1) - b2)
        four_pi_eta = 2 * (-Πmkap1 + η2 * kap0 - 0.25 * kite_area2 * (1 + 5 * r2 + b2))
    end
    return 2 * s0 + four_pi_eta
end


text:
include("ArrayFire/examples/ML/ANNDemo.jl")
using AF, ANNDemo

af = ArrayFire{CPU}()

setDevice(af, 0)
setSeed(af, 42)

runDemo(af)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
ENV["MPLBACKEND"] = "qt4agg"   # To fix a bug when connecting to Anaconda Python
using PyPlot, CSV

function make_scatterplot(lmm_results)
    fit_srn  = lmm_results[lmm_results[!,:rnntype].=="SRN", :]
    fit_gru  = lmm_results[lmm_results[!,:rnntype].=="GRU", :]
    fit_lstm = lmm_results[lmm_results[!,:rnntype].=="LSTM",:]

    scatter(-fit_srn[!,:avsurp], fit_srn[!,:chi2] .*sign.(fit_srn[!,:b_surp]),  s=10, c=:blue,  marker="x", label="SRN")
    scatter(-fit_gru[!,:avsurp], fit_gru[!,:chi2] .*sign.(fit_gru[!,:b_surp]),  s=10, c=:red,   marker="o", label="GRU")
    scatter(-fit_lstm[!,:avsurp],fit_lstm[!,:chi2].*sign.(fit_lstm[!,:b_surp]), s=10, c=:green, marker="^", label="LSTM")

    grid("on")
end

function make_gamplot(gam_curve)
    gam_srn  = gam_curve[gam_curve[!,:rnntype].=="SRN", :]
    gam_gru  = gam_curve[gam_curve[!,:rnntype].=="GRU", :]
    gam_lstm = gam_curve[gam_curve[!,:rnntype].=="LSTM", :]

    plot(gam_srn[!,:x],  gam_srn[!,:fit],  "b-",  label="SRN")
    plot(gam_gru[!,:x],  gam_gru[!,:fit],  "r:",  label="GRU")
    plot(gam_lstm[!,:x], gam_lstm[!,:fit], "g--", label="LSTM")

    fill_between(gam_srn[!,:x],  gam_srn[!,:fit] -2gam_srn[!,:se],  gam_srn[!,:fit] +2gam_srn[!,:se],  color=:blue,  alpha=0.3, linewidth=0)
    fill_between(gam_gru[!,:x],  gam_gru[!,:fit] -2gam_gru[!,:se],  gam_gru[!,:fit] +2gam_gru[!,:se],  color=:red,   alpha=0.3, linewidth=0)
    fill_between(gam_lstm[!,:x], gam_lstm[!,:fit] -2gam_lstm[!,:se], gam_lstm[!,:fit]+2gam_lstm[!,:se], color=:green, alpha=0.3, linewidth=0)

    grid("on")
end

#-------------------------------------------------------------------
lmm_SPR = CSV.read("lmm_SPR.csv", delim='\t')
lmm_ET  = CSV.read("lmm_ET.csv",  delim='\t')
lmm_EEG = CSV.read("lmm_EEG.csv", delim='\t')
gam_SPR = CSV.read("gam_SPR.csv", delim='\t')
gam_ET  = CSV.read("gam_ET.csv",  delim='\t')
gam_EEG = CSV.read("gam_EEG.csv", delim='\t')

lmm_EEG[!,:chi2] = -lmm_EEG[!,:chi2]
gam_EEG[!,:fit]  = -gam_EEG[!,:fit]

figure()

subplot(231)
make_scatterplot(lmm_SPR)
ylabel("Goodness-of-fit")
title("Self-paced reading time")

subplot(232)
make_scatterplot(lmm_ET)
title("Gaze duration")
legend(loc="upper left")

subplot(233)
make_scatterplot(lmm_EEG)
title("N400 size")
#legend(loc="upper right")

subplot(234)
make_gamplot(gam_SPR)
xlabel("Weighted average log-prob")
ylabel("Goodness-of-fit")

subplot(235)
make_gamplot(gam_ET)
xlabel("Weighted average log-prob")
legend(loc="upper left")

subplot(236)
make_gamplot(gam_EEG)
xlabel("Weighted average log-prob")
#legend(loc="upper right")

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
using OdinSon
using Distributions

# Graphics[Table[Circle[RandomReal[10, 2]], {100}]]
f1 = Canvas([Circle(rand(Uniform(0, 10), 2), 1, style=Style(fill=nothing)) for i = 1:100])
render(f1) # the render in the Mathematica notebook is implicit
f2 =Canvas(mapslices(p->Circle(p, 1, style=Style(fill=nothing)), rand(Uniform(0, 10), 100, 2), 2))
render(f2)
# So very verbose in comparison. Thow if instead we have like Mathematica implicit defaults
# which is very close, but I don't like the defaults as they don't feel natural to me,
# save for maybe having the style not fill by default
#=
```
Graphics[Table[Circle[RandomReal[10, 2]], {100}]]
Canvas([Circle(rand(Uniform(0, 10), 2)) for i = 1:100])
```
=#

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
########## functions ##########

function start_fishes(data)
    fishes = zeros(BigInt, 9)
    for day = 1:9
        fishes[day] = sum(data .== (day - 1))
    end
    return fishes
end

function live_another_day!(fishes)
    n_spawns = fishes[1]
    fishes[1:8] = fishes[2:9]
    fishes[9] = n_spawns
    fishes[7] += n_spawns
    return nothing
end

########## Script ##########

function main()

    # open file and store content into an input string
    input = open("inputs//day6.txt", "r") do io
        read(io, String)
    end

    # split input into a data vector of integers
    data = parse.(Int, String.(split(input, ",")))

    # initialize fishes table
    fishes = start_fishes(data)

    # part 1
    for _ = 1:80
        live_another_day!(fishes)
    end
    println("part 1: ", sum(fishes))

    # part 2
    for _ = 81:256
        live_another_day!(fishes)
    end
    println("part 1: ", sum(fishes))

end

main()


text:
@testset "transducer.jl" begin
    Random.seed!(1234)

    @testset "Basic sampling" begin
        N = 1_000
        local chain
        Logging.with_logger(TerminalLogger()) do
            xf = AbstractMCMC.Sample(MyModel(), MySampler();
                                     sleepy = true, logger = true)
            chain = withprogress(1:N; interval=1e-3) |> xf |> collect
        end

        # test output type and size
        @test chain isa Vector{<:MySample}
        @test length(chain) == N

        # test some statistical properties
        tail_chain = @view chain[2:end]
        @test mean(x.a for x in tail_chain) ≈ 0.5 atol=6e-2
        @test var(x.a for x in tail_chain) ≈ 1 / 12 atol=5e-3
        @test mean(x.b for x in tail_chain) ≈ 0.0 atol=5e-2
        @test var(x.b for x in tail_chain) ≈ 1 atol=6e-2
    end

    @testset "drop" begin
        xf = AbstractMCMC.Sample(MyModel(), MySampler())
        chain = 1:10 |> xf |> Drop(1) |> collect
        @test chain isa Vector{MySample{Float64,Float64}}
        @test length(chain) == 9
    end

    # Reproduce iterator example
    @testset "iterator example" begin
        # filter missing values and split transitions
        xf = opcompose(
            AbstractMCMC.Sample(MyModel(), MySampler()),
            OfType(MySample{Float64,Float64}),
            Map(x -> (x.a, x.b)),
        )
        as, bs = foldl(xf, 1:999; init = (Float64[], Float64[])) do (as, bs), (a, b)
            push!(as, a)
            push!(bs, b)
            as, bs
        end

        @test length(as) == length(bs) == 998

        @test mean(as) ≈ 0.5 atol=2e-2
        @test var(as) ≈ 1 / 12 atol=5e-3
        @test mean(bs) ≈ 0.0 atol=5e-2
        @test var(bs) ≈ 1 atol=5e-2
    end
end


text:
lines = readlines(ARGS[1])

function process_lines(lines)
    increased = 0
    prev = nothing
    for i in lines
        int_val = parse(Int, i)
        if isnothing(prev)
            prev = int_val
        else
            if int_val > prev
                increased += 1
            end
            prev = int_val
        end
    end
    return increased
end

function process_lines_2(lines)
    increased = 0
    prev = []
    for i in lines
        int_val = parse(Int, i)
        if length(prev) < 3
            append!(prev, int_val)
        else
            append!(prev, int_val)
            if sum(prev[2:4]) > sum(prev[1:3])
                increased += 1
            end
            popfirst!(prev)
        end
    end
    return increased
end

println(process_lines(lines))
println(process_lines_2(lines))


text:
export AlgAss, AlgAssElem, AlgGrp, AlgGrpElem, AlgMat, AlgMatElem

abstract type AbsAlgAss{T} <: Ring end

abstract type AbsAlgAssElem{T} <: RingElem end

################################################################################
#
#  AlgAss / AlgAssElem
#
################################################################################

# Associative algebras defined by structure constants (multiplication table)
mutable struct AlgAss{T} <: AbsAlgAss{T}
  base_ring::Ring
  mult_table::Array{T, 3} # e_i*e_j = sum_k mult_table[i, j, k]*e_k
  one::Vector{T}
  has_one::Bool
  iszero::Bool
  iscommutative::Int       # 0: don't know
                           # 1: known to be commutative
                           # 2: known to be not commutative
  basis#::Vector{AlgAssElem{T, AlgAss{T}}
  gens#::Vector{AlgAssElem{T, AlgAss{T}} # "Small" number of algebra generators
  trace_basis_elem::Vector{T}
  issimple::Int
  issemisimple::Int

  decomposition#::Vector{Tuple{AlgAss{T}, mor(AlgAss{T}, AlgAss{T})}
  center#Tuple{AlgAss{T}, mor(AlgAss{T}, AlgAss{T})
  maps_to_numberfields
  maximal_order

  # Constructor with default values
  function AlgAss{T}(R::Ring) where {T}
    A = new{T}()
    A.base_ring = R
    A.iscommutative = 0
    A.issimple = 0
    A.issemisimple = 0
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}, one::Vector{T}) where {T}
    A = AlgAss{T}(R)
    A.mult_table = mult_table
    A.one = one
    A.has_one = true
    A.iszero = (size(mult_table, 1) == 0)
    return A
  end

  function AlgAss{T}(R::Ring, mult_table::Array{T, 3}) where {T}
    A = AlgAss{T}(R)
    A.mult_table = mult_table
    A.iszero = (size(mult_table, 1) == 0)
    return A
  end
end

mutable struct AlgAssElem{T, S} <: AbsAlgAssElem{T}
  parent::S
  coeffs::Array{T, 1}

  function AlgAssElem{T, S}(A::S) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  function AlgAssElem{T, S}(A::AlgAss{T}) where {T, S}
    z = new{T, AlgAss{T}}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  # This does not make a copy of coeffs
  function AlgAssElem{T, S}(A::S, coeffs::Array{T, 1}) where {T, S}
    z = new{T, AlgAss{T}}()
    z.parent = A
    z.coeffs = coeffs
    return z
  end
end

################################################################################
#
#  AlgGrp / AlgGrpElem
#
################################################################################

# Group rings
mutable struct AlgGrp{T, S, R} <: AbsAlgAss{T}
  base_ring::Ring
  group::S
  group_to_base::Dict{R, Int}
  one::Vector{T}
  basis#::Vector{AlgAssElem{T, AlgAss{T}}
  gens#::Vector{AlgAssElem{T, AlgAss{T}} # "Small" number of algebra generators
  mult_table::Array{Int, 2} # b_i * b_j = b_(mult_table[i, j])
  iscommutative::Int
  trace_basis_elem::Array{T, 1}
  issimple::Int
  issemisimple::Int

  decomposition
  center
  maps_to_numberfields
  maximal_order

  function AlgGrp(K::Ring, G::GrpAbFinGen)
    A = AlgGrp(K, G, op = +)
    A.iscommutative = true
    return A
  end

  function AlgGrp(K::Ring, G; op = *)
    A = new{elem_type(K), typeof(G), elem_type(G)}()
    A.iscommutative = 0
    A.issimple = 0
    A.issemisimple = 0
    A.base_ring = K
    A.group = G
    d = Int(order(G))
    A.group_to_base = Dict{elem_type(G), Int}()
    A.mult_table = zeros(Int, d, d)

    for (i, g) in enumerate(G)
      A.group_to_base[deepcopy(g)] = i
    end

    v = Vector{elem_type(K)}(undef, d)
    for i in 1:d
      v[i] = zero(K)
    end
    v[1] = one(K)

    A.one = v

    for (i, g) in enumerate(G)
      for (j, h) in enumerate(G)
        l = op(g, h)
        A.mult_table[i, j] = A.group_to_base[l]
      end
    end

    @assert all(A.mult_table[1, i] == i for i in 1:dim(A))

    return A
  end
end

mutable struct AlgGrpElem{T, S} <: AbsAlgAssElem{T}
  parent::S
  coeffs::Array{T, 1}

  function AlgGrpElem{T, S}(A::S) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = Array{T, 1}(undef, size(A.mult_table, 1))
    for i = 1:length(z.coeffs)
      z.coeffs[i] = A.base_ring()
    end
    return z
  end

  function AlgGrpElem{T, S}(A::S, g::U) where {T, S, U}
    return A[A.group_to_base[g]]
  end

  # This does not make a copy of coeffs
  function AlgGrpElem{T, S}(A::S, coeffs::Array{T, 1}) where {T, S}
    z = new{T, S}()
    z.parent = A
    z.coeffs = coeffs
    return z
  end
end

################################################################################
#
#  AbsAlgAssIdl
#
################################################################################

# S is the type of the algebra, T = elem_type(S) and U is the type of matrices
# over the base ring of the algebra
mutable struct AbsAlgAssIdl{S, T, U}
  algebra::S
  basis::Vector{T}
  basis_matrix::U

  isleft::Int                      # 0 Not known
                                   # 1 Known to be a left ideal
                                   # 2 Known not to be a left ideal
  isright::Int                     # as for isleft

  iszero::Int

  function AbsAlgAssIdl{S, T, U}(A::S) where {S, T, U}
    I = new{S, T, U}()
    I.algebra = A
    I.isleft = 0
    I.isright = 0
    I.iszero = 0
    return I
  end

  function AbsAlgAssIdl{S, U}(A::S, M::U) where {S, U}
    I = new{S, elem_type(S), U}()
    I.algebra = A
    I.basis_matrix = M
    I.isleft = 0
    I.isright = 0
    I.iszero = 0
    return I
  end
end

################################################################################
#
#  AlgMat / AlgMatElem
#
################################################################################

# T == elem_type(base_ring), S == dense_matrix_type(coefficient_ring)
mutable struct AlgMat{T, S} <: AbsAlgAss{T}
  base_ring::Ring
  coefficient_ring::Ring
  one::S
  basis
  basis_matrix # matrix over the base_ring
  dim::Int
  degree::Int
  issimple::Int
  issemisimple::Int
  iscommutative::Int
  decomposition
  maximal_order
  mult_table::Array{T, 3} # e_i*e_j = sum_k mult_table[i, j, k]*e_k
  canonical_basis::Int # whether A[(j - 1)*n + i] == E_ij, where E_ij = (e_kl)_kl with e_kl = 1 if i =k and j = l and e_kl = 0 otherwise.
  center#Tuple{AlgAss{T}, mor(AlgAss{T}, AlgAss{T})
  trace_basis_elem::Vector{T}

  function AlgMat{T, S}(R::Ring) where {T, S}
    A = new{T, S}()
    A.base_ring = R
    A.coefficient_ring = R
    A.issimple = 0
    A.issemisimple = 0
    A.iscommutative = 0
    A.canonical_basis = 0
    return A
  end

  function AlgMat{T, S}(R1::Ring, R2::Ring) where {T, S}
    A = new{T, S}()
    A.base_ring = R1
    A.coefficient_ring = R2
    A.issimple = 0
    A.issemisimple = 0
    A.iscommutative = 0
    A.canonical_basis = 0
    return A
  end
end

mutable struct AlgMatElem{T, S, Mat} <: AbsAlgAssElem{T}
  parent::S
  matrix::Mat # over the coefficient ring of the parent
  coeffs::Vector{T} # over the base ring of the parent
  has_coeffs::Bool

  function AlgMatElem{T, S, Mat}(A::S) where {T, S, Mat}
    z = new{T, S, Mat}()
    z.parent = A
    z.matrix = zero_matrix(base_ring(A), degree(A), degree(A))
    z.has_coeffs = false
    return z
  end

  function AlgMatElem{T, S, Mat}(A::S, M::Mat) where {T, S, Mat}
    z = new{T, S, Mat}()
    z.parent = A
    z.matrix = M
    z.has_coeffs = false
    return z
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RegularSampling(n1, n2, ..., np)

Generate samples regularly using `n1` points along the first
parametric dimension, `n2` points along the second parametric
dimension, ..., `np` poitns along the last parametric dimension.

## Example

Sample sphere regularly with 360 longitudes and 180 latitudes:

```julia
sample(Sphere((0,0,0), 1), RegularSampling(360, 180))
```
"""
struct RegularSampling{N} <: ContinuousSamplingMethod
  sizes::NTuple{N,Int}
end

RegularSampling(sizes::Vararg{Int,N}) where {N} = RegularSampling(sizes)

function sample(segment::Segment, method::RegularSampling)
  n = first(method.sizes)
  a, b = extrema(segment)

  δ = (b - a) / (n - 1)

  (a + (i-1)*δ for i in 1:n)
end

function sample(box::Box, method::RegularSampling)
  sz = _adjust_sizes(method.sizes, paramdim(box))
  l, u = extrema(box)

  # origin and spacing
  or, sp = l, (u - l) ./ (sz .- 1)

  ivec(or + (ind.I .- 1) .* sp for ind in CartesianIndices(sz))
end

function sample(sphere::Sphere{2,T}, method::RegularSampling) where {T}
  sz = _adjust_sizes(method.sizes, paramdim(sphere))
  c, r = center(sphere), radius(sphere)

  V = T <: AbstractFloat ? T : Float64
  θmin, θmax = V(0), V(2π)
  δθ = (θmax - θmin) / sz[1]
  θrange = range(θmin, stop=θmax-δθ, length=sz[1])

  r⃗(θ) = Vec{2,V}(r*cos(θ), r*sin(θ))

  ivec(c + r⃗(θ) for θ in θrange)
end

# spherical coordinates in ISO 80000-2:2019 convention
function sample(sphere::Sphere{3,T}, method::RegularSampling) where {T}
  sz = _adjust_sizes(method.sizes, paramdim(sphere))
  c, r = center(sphere), radius(sphere)

  V = T <: AbstractFloat ? T : Float64
  θmin, θmax = V(0), V(π)
  φmin, φmax = V(0), V(2π)
  δθ = (θmax - θmin) / (sz[1] + 1)
  δφ = (φmax - φmin) / (sz[2]    )
  θrange = range(θmin+δθ, stop=θmax-δθ, length=sz[1])
  φrange = range(φmin, stop=φmax-δφ, length=sz[2])

  r⃗(θ, φ) = Vec{3,V}(r*sin(θ)*cos(φ), r*sin(θ)*sin(φ), r*cos(θ))

  ivec(c + r⃗(θ, φ) for θ in θrange, φ in φrange)
end

function sample(ball::Ball{Dim,T}, method::RegularSampling) where {Dim,T}
  sz = _adjust_sizes(method.sizes, paramdim(ball))
  c, r = center(ball), radius(ball)

  V = T <: AbstractFloat ? T : Float64
  smin, smax = V(0), V(1)
  δs = (smax - smin) / (last(sz) - 1)
  srange = range(smin+δs, stop=smax, length=last(sz))

  # reuse samples on the boundary
  points = sample(Sphere(c, r), RegularSampling(sz[1:Dim-1]))

  scale(p, s) = Point(s * coordinates(p))

  ivec(scale(p, s) for p in points, s in srange)
end

# helper function to adjust sizes to a given length
function _adjust_sizes(sizes, len)
  N, S = len, length(sizes)
  S != N ? ntuple(i -> i ≤ S ? sizes[i] : last(sizes), N) : sizes
end


text:
function get_json(url::String; retries=3)
    r = HTTP.get(url)
    if r.status != 200
        if retries > 0
            return get_json(url; retries=retries-1)
        end
        error("Unable to GET $(url)")
    end
    return JSON3.read(r.body)
end


text:
mutable struct EndsleyPatch{T<:IO} <: Patch
    io::T
    new_size::Int64
end
EndsleyPatch(io::IO) = EndsleyPatch(io, typemax(Int64))

format_magic(::Type{EndsleyPatch}) = "ENDSLEY/BSDIFF43"

function write_start(
    ::Type{EndsleyPatch},
    patch_io::IO,
    old_data::AbstractVector{UInt8},
    new_data::AbstractVector{UInt8},
)
    new_size = length(new_data)
    write_int(patch_io, new_size)
    stream = TranscodingStream(compressor(), patch_io)
    patch = EndsleyPatch(stream, new_size)
    finalizer(finalize_patch, patch)
    return patch
end

function read_start(::Type{EndsleyPatch}, patch_io::IO)
    new_size = read_int(patch_io)
    EndsleyPatch(TranscodingStream(decompressor(), patch_io), new_size)
end

function finalize_patch(patch::EndsleyPatch)
    if patch.io isa TranscodingStream
        # must be called to avoid leaking memory
        TranscodingStreams.changemode!(patch.io, :close)
    end
end

function write_finish(patch::EndsleyPatch)
    if patch.io isa TranscodingStream
        write(patch.io, TranscodingStreams.TOKEN_END)
    end
    flush(patch.io)
    finalize(patch)
end

function encode_control(
    patch::EndsleyPatch,
    diff_size::Int,
    copy_size::Int,
    skip_size::Int,
)
    write_int(patch.io, diff_size)
    write_int(patch.io, copy_size)
    write_int(patch.io, skip_size)
end

function decode_control(patch::EndsleyPatch)
    eof(patch.io) && return nothing
    diff_size = read_int(patch.io)
    copy_size = read_int(patch.io)
    skip_size = read_int(patch.io)
    return diff_size, copy_size, skip_size
end

function encode_diff(
    patch::EndsleyPatch,
    diff_size::Int,
    new::AbstractVector{UInt8}, new_pos::Int,
    old::AbstractVector{UInt8}, old_pos::Int,
)
    for i = 1:diff_size
        write(patch.io, new[new_pos + i] - old[old_pos + i])
    end
end

function decode_diff(
    patch::EndsleyPatch,
    diff_size::Int,
    new::IO,
    old::AbstractVector{UInt8},
    old_pos::Int,
)
    for i = 1:diff_size
        write(new, old[old_pos + i] + read(patch.io, UInt8))
    end
end

function encode_data(
    patch::EndsleyPatch,
    copy_size::Int,
    new::AbstractVector{UInt8}, pos::Int,
)
    for i = 1:copy_size
        write(patch.io, new[pos + i])
    end
end

function decode_data(
    patch::EndsleyPatch,
    copy_size::Int,
    new::IO,
)
    for i = 1:copy_size
        write(new, read(patch.io, UInt8))
    end
end


text:
function su2u1_grossneveu(;g2SPT=0,g2AFM=0)
    ph       = Rep[SU₂×U₁]( (1//2,0)=>1, (0,-1)=>1, (0,1)=>1 )
    bigonleg = Rep[SU₂×U₁]( (0,0)=>1, (1//2,-1)=>1, (1//2,1)=>1 )
    unit     = oneunit(ph)

    LK = TensorMap(ones, ComplexF64, unit*ph, bigonleg*ph)
    blocks(LK)[Irrep[SU₂](0)⊠Irrep[U₁](-1)]    =  [im*2/sqrt(2) 1]
    blocks(LK)[Irrep[SU₂](1//2)⊠Irrep[U₁](0)] =  [1. im -im]
    blocks(LK)[Irrep[SU₂](0)⊠Irrep[U₁](1)]    =  [im*2/sqrt(2) 1]

    RK = TensorMap(ones, ComplexF64, bigonleg*ph, unit*ph)
    blocks(RK)[Irrep[SU₂](0)⊠Irrep[U₁](-1)][:]    =  [2/sqrt(2) 1][:]
    blocks(RK)[Irrep[SU₂](1//2)⊠Irrep[U₁](0)][:] =  [1 -1 1][:]
    blocks(RK)[Irrep[SU₂](0)⊠Irrep[U₁](1)][:]    =  [2/sqrt(2) 1][:]

    Cplus = TensorMap(ones, ComplexF64, bigonleg*ph, bigonleg*ph)
    blocks(Cplus)[Irrep[SU₂](1//2)⊠Irrep[U₁](0)]  = [0 im*0.5 -im*0.5; -0.5 0 0; 0.5 0 0]
    blocks(Cplus)[Irrep[SU₂](0)⊠Irrep[U₁](-1)]     = [0 0.5*sqrt(2); im*0.5*sqrt(2) 0]
    blocks(Cplus)[Irrep[SU₂](1)⊠Irrep[U₁](-1)]     = zeros(1,1)
    blocks(Cplus)[Irrep[SU₂](0)⊠Irrep[U₁](1)]     = [0 0.5*sqrt(2); im*0.5*sqrt(2) 0]
    blocks(Cplus)[Irrep[SU₂](1)⊠Irrep[U₁](1)]     = zeros(1,1)
    blocks(Cplus)[Irrep[SU₂](1//2)⊠Irrep[U₁](-2)] = zeros(1,1)
    blocks(Cplus)[Irrep[SU₂](1//2)⊠Irrep[U₁](2)]  = zeros(1,1)

    Cmin = TensorMap(ones, ComplexF64, bigonleg*ph, bigonleg*ph)
    blocks(Cmin)[Irrep[SU₂](1//2)⊠Irrep[U₁](0)]  = conj([0 im*0.5 -im*0.5; -0.5 0 0; 0.5 0 0])
    blocks(Cmin)[Irrep[SU₂](0)⊠Irrep[U₁](-1)]     = conj([0 0.5*sqrt(2); im*0.5*sqrt(2) 0])
    blocks(Cmin)[Irrep[SU₂](1)⊠Irrep[U₁](-1)]     = zeros(1,1)
    blocks(Cmin)[Irrep[SU₂](0)⊠Irrep[U₁](1)]     = conj([0 0.5*sqrt(2); im*0.5*sqrt(2) 0])
    blocks(Cmin)[Irrep[SU₂](1)⊠Irrep[U₁](1)]     = zeros(1,1)
    blocks(Cmin)[Irrep[SU₂](1//2)⊠Irrep[U₁](-2)] = zeros(1,1)
    blocks(Cmin)[Irrep[SU₂](1//2)⊠Irrep[U₁](2)]  = zeros(1,1)

    f1 = isomorphism(fuse(unit, unit), unit*unit)
    f2 = isomorphism(bigonleg*bigonleg, fuse(bigonleg, bigonleg))
    f3 = isomorphism(fuse(bigonleg, bigonleg), bigonleg*bigonleg)
    f4 = isomorphism(unit*unit, fuse(unit, unit))

    @tensor Ldiffsq[-1 -2;-3 -4] := f1[-1,1,2]*LK[1,3,5,-4]*LK[2,-2,4,3]*f2[5,4,-3]
    @tensor Cdiffsq[-1 -2;-3 -4] := f3[-1,1,2]*Cmin[1,3,5,-4]*Cmin[2,-2,4,3]*f2[5,4,-3]
    @tensor Rdiffsq[-1 -2;-3 -4] := f3[-1,1,2]*RK[1,3,5,-4]*RK[2,-2,4,3]*f4[5,4,-3]

    #and now with the extra O(4) breaking part ie the O operator
    O_op = TensorMap(zeros, ComplexF64, unit*ph, unit*ph)
    blocks(O_op)[Irrep[SU₂](1//2)⊠Irrep[U₁](0)] = -zeros(1,1)
    blocks(O_op)[Irrep[SU₂](0)⊠Irrep[U₁](-1)]    =  -1*ones(1,1)
    blocks(O_op)[Irrep[SU₂](0)⊠Irrep[U₁](1)]    =  1*ones(1,1)

    MPOHamiltonian([LK, Cplus, RK]) +
    MPOHamiltonian([-0.25*g2SPT^2*Ldiffsq, Cdiffsq, Rdiffsq]) +
    MPOHamiltonian([-0.5*g2AFM^2*O_op*O_op]) +
    MPOHamiltonian([+0.5*g2AFM^2*O_op, O_op])
end

function su2u1_orderpars()
    ph       = ℂ[SU₂×U₁]( (1//2,0)=>1, (0,-1)=>1, (0,1)=>1 )
    onleg    = ℂ[SU₂×U₁]( (1//2,-1)=>1, (1//2,1)=>1 )
    unit     = oneunit(ph)

    LK = TensorMap(ones, ComplexF64, unit*ph, onleg*ph)
    blocks(LK)[SU₂(0)×U₁(-1)]   =  ones(1,1)*im*2/sqrt(2)
    blocks(LK)[SU₂(1//2)×U₁(0)] =  [im -im]
    blocks(LK)[SU₂(0)×U₁(1)]    =  ones(1,1)*im*2/sqrt(2)

    RK = TensorMap(ones, ComplexF64, onleg*ph, unit*ph)
    blocks(RK)[SU₂(0)×U₁(-1)][:]    =  ones(1,1)*2/sqrt(2)
    blocks(RK)[SU₂(1//2)×U₁(0)][:] =  [-1, 1][:]
    blocks(RK)[SU₂(0)×U₁(1)][:]    =  ones(1,1)*2/sqrt(2)

    #and now with the extra O(4) breaking part ie the O operator
    O_op = TensorMap(zeros, ComplexF64, unit*ph, unit*ph)
    blocks(O_op)[SU₂(1//2)×U₁(0)] = -zeros(1,1)
    blocks(O_op)[SU₂(0)×U₁(-1)]    =  -1*ones(1,1)
    blocks(O_op)[SU₂(0)×U₁(1)]    =  1*ones(1,1)

    return [LK, RK] , [O_op]
end


text:
function transition(interpreter::Interpreter,
                    domain::Domain, state::State, action::Term; options...)
    return execute(interpreter, domain, state, action; options...)
end

function transition!(interpreter::Interpreter,
                     domain::Domain, state::State, action::Term; options...)
    return execute!(interpreter, domain, state, action; options...)
end


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

module TestTridiagonal

using Test, LinearAlgebra, Random

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")

isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
using .Main.Quaternions

include("testutils.jl") # test_approx_eq_modphase

#Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences
function test_approx_eq_vecs(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) where {S<:Real,T<:Real}
    n = size(a, 1)
    @test n==size(b,1) && size(a,2)==size(b,2)
    error===nothing && (error=n^3*(eps(S)+eps(T)))
    for i=1:n
        ev1, ev2 = a[:,i], b[:,i]
        deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2)))
        if !isnan(deviation)
            @test deviation ≈ 0.0 atol=error
        end
    end
end

@testset for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)
    n = 12 #Size of matrix problem to test
    Random.seed!(123)
    if elty == Int
        Random.seed!(61516384)
        d = rand(1:100, n)
        dl = -rand(0:10, n-1)
        du = -rand(0:10, n-1)
        v = rand(1:100, n)
        B = rand(1:100, n, 2)
        a = rand(1:100, n-1)
        b = rand(1:100, n)
        c = rand(1:100, n-1)
    else
        d = convert(Vector{elty}, 1 .+ randn(n))
        dl = convert(Vector{elty}, randn(n - 1))
        du = convert(Vector{elty}, randn(n - 1))
        v = convert(Vector{elty}, randn(n))
        B = convert(Matrix{elty}, randn(n, 2))
        a = convert(Vector{elty}, randn(n - 1))
        b = convert(Vector{elty}, randn(n))
        c = convert(Vector{elty}, randn(n - 1))
        if elty <: Complex
            a += im*convert(Vector{elty}, randn(n - 1))
            b += im*convert(Vector{elty}, randn(n))
            c += im*convert(Vector{elty}, randn(n - 1))
        end
    end
    @test_throws DimensionMismatch SymTridiagonal(dl, fill(elty(1), n+1))
    @test_throws ArgumentError SymTridiagonal(rand(n, n))
    @test_throws ArgumentError Tridiagonal(dl, dl, dl)
    @test_throws ArgumentError convert(SymTridiagonal{elty}, Tridiagonal(dl, d, du))

    if elty != Int
        @testset "issue #1490" begin
            @test det(fill(elty(1),3,3)) ≈ zero(elty) atol=3*eps(real(one(elty)))
            @test det(SymTridiagonal(elty[],elty[])) == one(elty)
        end
    end

    @testset "constructor" begin
        for (x, y) in ((d, dl), (GenericArray(d), GenericArray(dl)))
            ST = (SymTridiagonal(x, y))::SymTridiagonal{elty, typeof(x)}
            @test ST == Matrix(ST)
            @test ST.dv === x
            @test ST.ev === y
            TT = (Tridiagonal(y, x, y))::Tridiagonal{elty, typeof(x)}
            @test TT == Matrix(TT)
            @test TT.dl === y
            @test TT.d  === x
            @test TT.du === y
        end
        ST = SymTridiagonal{elty}([1,2,3,4], [1,2,3])
        @test eltype(ST) == elty
        @test SymTridiagonal{elty, Vector{elty}}(ST) === ST
        @test SymTridiagonal{Int64, Vector{Int64}}(ST) isa SymTridiagonal{Int64, Vector{Int64}}
        TT = Tridiagonal{elty}([1,2,3], [1,2,3,4], [1,2,3])
        @test eltype(TT) == elty
        ST = SymTridiagonal{elty,Vector{elty}}(d, GenericArray(dl))
        @test isa(ST, SymTridiagonal{elty,Vector{elty}})
        TT = Tridiagonal{elty,Vector{elty}}(GenericArray(dl), d, GenericArray(dl))
        @test isa(TT, Tridiagonal{elty,Vector{elty}})
        @test_throws MethodError SymTridiagonal(d, GenericArray(dl))
        @test_throws MethodError SymTridiagonal(GenericArray(d), dl)
        @test_throws MethodError Tridiagonal(GenericArray(dl), d, GenericArray(dl))
        @test_throws MethodError Tridiagonal(dl, GenericArray(d), dl)
        @test_throws MethodError SymTridiagonal{elty}(d, GenericArray(dl))
        @test_throws MethodError Tridiagonal{elty}(GenericArray(dl), d,GenericArray(dl))
        STI = SymTridiagonal([1,2,3,4], [1,2,3])
        TTI = Tridiagonal([1,2,3], [1,2,3,4], [1,2,3])
        TTI2 = Tridiagonal([1,2,3], [1,2,3,4], [1,2,3], [1,2])
        @test SymTridiagonal(STI) === STI
        @test Tridiagonal(TTI)    === TTI
        @test Tridiagonal(TTI2)   === TTI2
        @test isa(SymTridiagonal{elty}(STI), SymTridiagonal{elty})
        @test isa(Tridiagonal{elty}(TTI), Tridiagonal{elty})
        TTI2y = Tridiagonal{elty}(TTI2)
        @test isa(TTI2y, Tridiagonal{elty})
        @test TTI2y.du2 == convert(Vector{elty}, [1,2])
    end
    @testset "interconversion of Tridiagonal and SymTridiagonal" begin
        @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl)
        @test SymTridiagonal(d, dl) == Tridiagonal(dl, d, dl)
        @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du)
        @test SymTridiagonal(d, dl) + Tridiagonal(dl, d, du) == Tridiagonal(dl + dl, d+d, dl+du)
        @test convert(SymTridiagonal,Tridiagonal(SymTridiagonal(d, dl))) == SymTridiagonal(d, dl)
        @test Array(convert(SymTridiagonal{ComplexF32},Tridiagonal(SymTridiagonal(d, dl)))) == convert(Matrix{ComplexF32}, SymTridiagonal(d, dl))
    end
    @testset "tril/triu" begin
        zerosd = fill!(similar(d), 0)
        zerosdl = fill!(similar(dl), 0)
        zerosdu = fill!(similar(du), 0)
        @test_throws ArgumentError tril!(SymTridiagonal(d, dl), -n - 2)
        @test_throws ArgumentError tril!(SymTridiagonal(d, dl), n)
        @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), -n - 2)
        @test_throws ArgumentError tril!(Tridiagonal(dl, d, du), n)
        @test tril(SymTridiagonal(d,dl))    == Tridiagonal(dl,d,zerosdl)
        @test tril(SymTridiagonal(d,dl),1)  == Tridiagonal(dl,d,dl)
        @test tril(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,zerosd,zerosdl)
        @test tril(SymTridiagonal(d,dl),-2) == Tridiagonal(zerosdl,zerosd,zerosdl)
        @test tril(Tridiagonal(dl,d,du))    == Tridiagonal(dl,d,zerosdu)
        @test tril(Tridiagonal(dl,d,du),1)  == Tridiagonal(dl,d,du)
        @test tril(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,zerosd,zerosdu)
        @test tril(Tridiagonal(dl,d,du),-2) == Tridiagonal(zerosdl,zerosd,zerosdu)

        @test_throws ArgumentError triu!(SymTridiagonal(d, dl), -n)
        @test_throws ArgumentError triu!(SymTridiagonal(d, dl), n + 2)
        @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), -n)
        @test_throws ArgumentError triu!(Tridiagonal(dl, d, du), n + 2)
        @test triu(SymTridiagonal(d,dl))    == Tridiagonal(zerosdl,d,dl)
        @test triu(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,d,dl)
        @test triu(SymTridiagonal(d,dl),1)  == Tridiagonal(zerosdl,zerosd,dl)
        @test triu(SymTridiagonal(d,dl),2)  == Tridiagonal(zerosdl,zerosd,zerosdl)
        @test triu(Tridiagonal(dl,d,du))    == Tridiagonal(zerosdl,d,du)
        @test triu(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,d,du)
        @test triu(Tridiagonal(dl,d,du),1)  == Tridiagonal(zerosdl,zerosd,du)
        @test triu(Tridiagonal(dl,d,du),2)  == Tridiagonal(zerosdl,zerosd,zerosdu)

        @test !istril(SymTridiagonal(d,dl))
        @test istril(SymTridiagonal(d,zerosdl))
        @test !istril(SymTridiagonal(d,dl),-2)
        @test !istriu(SymTridiagonal(d,dl))
        @test istriu(SymTridiagonal(d,zerosdl))
        @test !istriu(SymTridiagonal(d,dl),2)
        @test istriu(Tridiagonal(zerosdl,d,du))
        @test !istriu(Tridiagonal(dl,d,zerosdu))
        @test istriu(Tridiagonal(zerosdl,zerosd,du),1)
        @test !istriu(Tridiagonal(dl,d,zerosdu),2)
        @test istril(Tridiagonal(dl,d,zerosdu))
        @test !istril(Tridiagonal(zerosdl,d,du))
        @test istril(Tridiagonal(dl,zerosd,zerosdu),-1)
        @test !istril(Tridiagonal(dl,d,zerosdu),-2)

        @test isdiag(SymTridiagonal(d,zerosdl))
        @test !isdiag(SymTridiagonal(d,dl))
        @test isdiag(Tridiagonal(zerosdl,d,zerosdu))
        @test !isdiag(Tridiagonal(dl,d,zerosdu))
        @test !isdiag(Tridiagonal(zerosdl,d,du))
        @test !isdiag(Tridiagonal(dl,d,du))

        # Test methods that could fail due to dv and ev having the same length
        # see #41089

        badev = zero(d)
        badev[end] = 1
        S = SymTridiagonal(d, badev)

        @test istriu(S, -2)
        @test istriu(S, 0)
        @test !istriu(S, 2)

        @test isdiag(S)
    end

    @testset "iszero and isone" begin
        Tzero = Tridiagonal(zeros(elty, 9), zeros(elty, 10), zeros(elty, 9))
        Tone = Tridiagonal(zeros(elty, 9), ones(elty, 10), zeros(elty, 9))
        Tmix = Tridiagonal(zeros(elty, 9), zeros(elty, 10), zeros(elty, 9))
        Tmix[end, end] = one(elty)

        Szero = SymTridiagonal(zeros(elty, 10), zeros(elty, 9))
        Sone = SymTridiagonal(ones(elty, 10), zeros(elty, 9))
        Smix = SymTridiagonal(zeros(elty, 10), zeros(elty, 9))
        Smix[end, end] = one(elty)

        @test iszero(Tzero)
        @test !isone(Tzero)
        @test !iszero(Tone)
        @test isone(Tone)
        @test !iszero(Tmix)
        @test !isone(Tmix)

        @test iszero(Szero)
        @test !isone(Szero)
        @test !iszero(Sone)
        @test isone(Sone)
        @test !iszero(Smix)
        @test !isone(Smix)

        badev = zeros(elty, 3)
        badev[end] = 1

        @test isone(SymTridiagonal(ones(elty, 3), badev))
        @test iszero(SymTridiagonal(zeros(elty, 3), badev))
    end

    @testset for mat_type in (Tridiagonal, SymTridiagonal)
        A = mat_type == Tridiagonal ? mat_type(dl, d, du) : mat_type(d, dl)
        fA = map(elty <: Complex ? ComplexF64 : Float64, Array(A))
        @testset "similar, size, and copyto!" begin
            B = similar(A)
            @test size(B) == size(A)
            copyto!(B, A)
            @test B == A
            @test isa(similar(A), mat_type{elty})
            @test isa(similar(A, Int), mat_type{Int})
            @test isa(similar(A, (3, 2)), Matrix)
            @test isa(similar(A, Int, (3, 2)), Matrix{Int})
            @test size(A, 3) == 1
            @test size(A, 1) == n
            @test size(A) == (n, n)
            @test_throws ArgumentError size(A, 0)
        end
        @testset "getindex" begin
            @test_throws BoundsError A[n + 1, 1]
            @test_throws BoundsError A[1, n + 1]
            @test A[1, n] == convert(elty, 0.0)
            @test A[1, 1] == d[1]
        end
        @testset "setindex!" begin
            @test_throws BoundsError A[n + 1, 1] = 0 # test bounds check
            @test_throws BoundsError A[1, n + 1] = 0 # test bounds check
            @test_throws ArgumentError A[1, 3]   = 1 # test assignment off the main/sub/super diagonal
            if mat_type == Tridiagonal
                @test (A[3, 3] = A[3, 3]; A == fA) # test assignment on the main diagonal
                @test (A[3, 2] = A[3, 2]; A == fA) # test assignment on the subdiagonal
                @test (A[2, 3] = A[2, 3]; A == fA) # test assignment on the superdiagonal
                @test ((A[1, 3] = 0) == 0; A == fA) # test zero assignment off the main/sub/super diagonal
            else # mat_type is SymTridiagonal
                @test ((A[3, 3] = A[3, 3]) == A[3, 3]; A == fA) # test assignment on the main diagonal
                @test_throws ArgumentError A[3, 2] = 1 # test assignment on the subdiagonal
                @test_throws ArgumentError A[2, 3] = 1 # test assignment on the superdiagonal
            end
        end
        @testset "diag" begin
            @test (@inferred diag(A))::typeof(d) == d
            @test (@inferred diag(A, 0))::typeof(d) == d
            @test (@inferred diag(A, 1))::typeof(d) == (mat_type == Tridiagonal ? du : dl)
            @test (@inferred diag(A, -1))::typeof(d) == dl
            @test (@inferred diag(A, n-1))::typeof(d) == zeros(elty, 1)
            @test_throws ArgumentError diag(A, -n - 1)
            @test_throws ArgumentError diag(A, n + 1)
            GA = mat_type == Tridiagonal ? mat_type(GenericArray.((dl, d, du))...) : mat_type(GenericArray.((d, dl))...)
            @test (@inferred diag(GA))::typeof(GenericArray(d)) == GenericArray(d)
            @test (@inferred diag(GA, -1))::typeof(GenericArray(d)) == GenericArray(dl)
        end
        @testset "Idempotent tests" begin
            for func in (conj, transpose, adjoint)
                @test func(func(A)) == A
            end
        end
        @testset "permutedims(::[Sym]Tridiagonal)" begin
            @test permutedims(permutedims(A)) === A
            @test permutedims(A) == transpose.(transpose(A))
            @test permutedims(A, [1, 2]) === A
            @test permutedims(A, (2, 1)) == permutedims(A)
        end
        if elty != Int
            @testset "Simple unary functions" begin
                for func in (det, inv)
                    @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(real(one(elty))))
                end
            end
        end
        ds = mat_type == Tridiagonal ? (dl, d, du) : (d, dl)
        for f in (real, imag)
            @test f(A)::mat_type == mat_type(map(f, ds)...)
        end
        if elty <: Real
            for f in (round, trunc, floor, ceil)
                fds = [f.(d) for d in ds]
                @test f.(A)::mat_type == mat_type(fds...)
                @test f.(Int, A)::mat_type == f.(Int, fA)
            end
        end
        fds = [abs.(d) for d in ds]
        @test abs.(A)::mat_type == mat_type(fds...)
        @testset "Multiplication with strided matrix/vector" begin
            @test (x = fill(1.,n); A*x ≈ Array(A)*x)
            @test (X = fill(1.,n,2); A*X ≈ Array(A)*X)
        end
        @testset "Binary operations" begin
            B = mat_type == Tridiagonal ? mat_type(a, b, c) : mat_type(b, a)
            fB = map(elty <: Complex ? ComplexF64 : Float64, Array(B))
            for op in (+, -, *)
                @test Array(op(A, B)) ≈ op(fA, fB)
            end
            α = rand(elty)
            @test Array(α*A) ≈ α*Array(A)
            @test Array(A*α) ≈ Array(A)*α
            @test Array(A/α) ≈ Array(A)/α

            @testset "Matmul with Triangular types" begin
                @test A*LinearAlgebra.UnitUpperTriangular(Matrix(1.0I, n, n)) ≈ fA
                @test A*LinearAlgebra.UnitLowerTriangular(Matrix(1.0I, n, n)) ≈ fA
                @test A*UpperTriangular(Matrix(1.0I, n, n)) ≈ fA
                @test A*LowerTriangular(Matrix(1.0I, n, n)) ≈ fA
            end
            @testset "mul! errors" begin
                Cnn, Cnm, Cmn = Matrix{elty}.(undef, ((n,n), (n,n+1), (n+1,n)))
                @test_throws DimensionMismatch LinearAlgebra.mul!(Cnn,A,Cnm)
                @test_throws DimensionMismatch LinearAlgebra.mul!(Cnn,A,Cmn)
                @test_throws DimensionMismatch LinearAlgebra.mul!(Cnn,B,Cmn)
                @test_throws DimensionMismatch LinearAlgebra.mul!(Cmn,B,Cnn)
                @test_throws DimensionMismatch LinearAlgebra.mul!(Cnm,B,Cnn)
            end
        end
        @testset "Negation" begin
            mA = -A
            @test mA isa mat_type
            @test -mA == A
        end
        if mat_type == SymTridiagonal
            @testset "Tridiagonal/SymTridiagonal mixing ops" begin
                B = convert(Tridiagonal{elty}, A)
                @test B == A
                @test B + A == A + B
                @test B - A == A - B
            end
            if elty <: LinearAlgebra.BlasReal
                @testset "Eigensystems" begin
                    zero, infinity = convert(elty, 0), convert(elty, Inf)
                    @testset "stebz! and stein!" begin
                        w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a)
                        evecs = LAPACK.stein!(b, a, w)

                        (e, v) = eigen(SymTridiagonal(b, a))
                        @test e ≈ w
                        test_approx_eq_vecs(v, evecs)
                    end
                    @testset "stein! call using iblock and isplit" begin
                        w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a)
                        evecs = LAPACK.stein!(b, a, w, iblock, isplit)
                        test_approx_eq_vecs(v, evecs)
                    end
                    @testset "stegr! call with index range" begin
                        F = eigen(SymTridiagonal(b, a),1:2)
                        fF = eigen(Symmetric(Array(SymTridiagonal(b, a))),1:2)
                        test_approx_eq_modphase(F.vectors, fF.vectors)
                        @test F.values ≈ fF.values
                    end
                    @testset "stegr! call with value range" begin
                        F = eigen(SymTridiagonal(b, a),0.0,1.0)
                        fF = eigen(Symmetric(Array(SymTridiagonal(b, a))),0.0,1.0)
                        test_approx_eq_modphase(F.vectors, fF.vectors)
                        @test F.values ≈ fF.values
                    end
                    @testset "eigenvalues/eigenvectors of symmetric tridiagonal" begin
                        if elty === Float32 || elty === Float64
                            DT, VT = @inferred eigen(A)
                            @inferred eigen(A, 2:4)
                            @inferred eigen(A, 1.0, 2.0)
                            D, Vecs = eigen(fA)
                            @test DT ≈ D
                            @test abs.(VT'Vecs) ≈ Matrix(elty(1)I, n, n)
                            test_approx_eq_modphase(eigvecs(A), eigvecs(fA))
                            #call to LAPACK.stein here
                            test_approx_eq_modphase(eigvecs(A,eigvals(A)),eigvecs(A))
                        elseif elty != Int
                            # check that undef is determined accurately even if type inference
                            # bails out due to the number of try/catch blocks in this code.
                            @test_throws UndefVarError fA
                        end
                    end
                end
            end
            if elty <: Real
                Ts = SymTridiagonal(d, dl)
                Fs = Array(Ts)
                Tldlt = factorize(Ts)
                @testset "symmetric tridiagonal" begin
                    @test_throws DimensionMismatch Tldlt\rand(elty,n+1)
                    @test size(Tldlt) == size(Ts)
                    if elty <: AbstractFloat
                        @test LinearAlgebra.LDLt{elty,SymTridiagonal{elty,Vector{elty}}}(Tldlt) === Tldlt
                        @test LinearAlgebra.LDLt{elty}(Tldlt) === Tldlt
                        @test typeof(convert(LinearAlgebra.LDLt{Float32,Matrix{Float32}},Tldlt)) ==
                            LinearAlgebra.LDLt{Float32,Matrix{Float32}}
                        @test typeof(convert(LinearAlgebra.LDLt{Float32},Tldlt)) ==
                            LinearAlgebra.LDLt{Float32,SymTridiagonal{Float32,Vector{Float32}}}
                    end
                    for vv in (copy(v), view(v, 1:n))
                        invFsv = Fs\vv
                        x = Ts\vv
                        @test x ≈ invFsv
                        @test Array(Tldlt) ≈ Fs
                    end

                    @testset "similar" begin
                        @test isa(similar(Ts), SymTridiagonal{elty})
                        @test isa(similar(Ts, Int), SymTridiagonal{Int})
                        @test isa(similar(Ts, (3, 2)), Matrix)
                        @test isa(similar(Ts, Int, (3, 2)), Matrix{Int})
                    end

                    @test first(logabsdet(Tldlt)) ≈ first(logabsdet(Fs))
                    @test last(logabsdet(Tldlt))  ≈ last(logabsdet(Fs))
                    # just test that the det method exists. The numerical value of the
                    # determinant is unreliable
                    det(Tldlt)
                end
            end
        else # mat_type is Tridiagonal
            @testset "tridiagonal linear algebra" begin
                for (BB, vv) in ((copy(B), copy(v)), (view(B, 1:n, 1), view(v, 1:n)))
                    @test A*vv ≈ fA*vv
                    invFv = fA\vv
                    @test A\vv ≈ invFv
                    # @test Base.solve(T,v) ≈ invFv
                    # @test Base.solve(T, B) ≈ F\B
                    Tlu = factorize(A)
                    x = Tlu\vv
                    @test x ≈ invFv
                end
            end
        end
        @testset "generalized dot" begin
            x = fill(convert(elty, 1), n)
            y = fill(convert(elty, 1), n)
            @test dot(x, A, y) ≈ dot(A'x, y)
        end
    end
end

@testset "SymTridiagonal block matrix" begin
    M = [1 2; 2 4]
    n = 5
    A = SymTridiagonal(fill(M, n), fill(M, n-1))
    @test @inferred A[1,1] == Symmetric(M)
    @test @inferred A[1,2] == M
    @test @inferred A[2,1] == transpose(M)
    @test @inferred diag(A, 1) == fill(M, n-1)
    @test @inferred diag(A, 0) == fill(Symmetric(M), n)
    @test @inferred diag(A, -1) == fill(transpose(M), n-1)
    @test_throws ArgumentError diag(A, -2)
    @test_throws ArgumentError diag(A, 2)
    @test_throws ArgumentError diag(A, n+1)
    @test_throws ArgumentError diag(A, -n-1)
end

@testset "Issue 12068" begin
    @test SymTridiagonal([1, 2], [0])^3 == [1 0; 0 8]
end

@testset "convert for SymTridiagonal" begin
    STF32 = SymTridiagonal{Float32}(fill(1f0, 5), fill(1f0, 4))
    @test convert(SymTridiagonal{Float64}, STF32)::SymTridiagonal{Float64} == STF32
    @test convert(AbstractMatrix{Float64}, STF32)::SymTridiagonal{Float64} == STF32
end

@testset "constructors from matrix" begin
    @test SymTridiagonal([1 2 3; 2 5 6; 0 6 9]) == [1 2 0; 2 5 6; 0 6 9]
    @test Tridiagonal([1 2 3; 4 5 6; 7 8 9]) == [1 2 0; 4 5 6; 0 8 9]
end

@testset "constructors with range and other abstract vectors" begin
    @test SymTridiagonal(1:3, 1:2) == [1 1 0; 1 2 2; 0 2 3]
    @test Tridiagonal(4:5, 1:3, 1:2) == [1 1 0; 4 2 2; 0 5 3]
end

@testset "Issue #26994 (and the empty case)" begin
    T = SymTridiagonal([1.0],[3.0])
    x = ones(1)
    @test T*x == ones(1)
    @test SymTridiagonal(ones(0), ones(0)) * ones(0, 2) == ones(0, 2)
end

@testset "Issue 29630" begin
    function central_difference_discretization(N; dfunc = x -> 12x^2 - 2N^2,
                                               dufunc = x -> N^2 + 4N*x,
                                               dlfunc = x -> N^2 - 4N*x,
                                               bfunc = x -> 114ℯ^-x * (1 + 3x),
                                               b0 = 0, bf = 57/ℯ,
                                               x0 = 0, xf = 1)
        h = 1/N
        d, du, dl, b = map(dfunc, (x0+h):h:(xf-h)), map(dufunc, (x0+h):h:(xf-2h)),
                       map(dlfunc, (x0+2h):h:(xf-h)), map(bfunc, (x0+h):h:(xf-h))
        b[1] -= dlfunc(x0)*b0     # subtract the boundary term
        b[end] -= dufunc(xf)*bf   # subtract the boundary term
        Tridiagonal(dl, d, du), b
    end

    A90, b90 = central_difference_discretization(90)

    @test A90\b90 ≈ inv(A90)*b90
end

@testset "singular values of SymTridiag" begin
    @test svdvals(SymTridiagonal([-4,2,3], [0,0])) ≈ [4,3,2]
    @test svdvals(SymTridiagonal(collect(0.:10.), zeros(10))) ≈ reverse(0:10)
    @test svdvals(SymTridiagonal([1,2,1], [1,1])) ≈ [3,1,0]
    # test that dependent methods such as `cond` also work
    @test cond(SymTridiagonal([1,2,3], [0,0])) ≈ 3
end

@testset "sum, mapreduce" begin
    T = Tridiagonal([1,2], [1,2,3], [7,8])
    Tdense = Matrix(T)
    S = SymTridiagonal([1,2,3], [1,2])
    Sdense = Matrix(S)
    @test sum(T) == 24
    @test sum(S) == 12
    @test_throws ArgumentError sum(T, dims=0)
    @test sum(T, dims=1) == sum(Tdense, dims=1)
    @test sum(T, dims=2) == sum(Tdense, dims=2)
    @test sum(T, dims=3) == sum(Tdense, dims=3)
    @test typeof(sum(T, dims=1)) == typeof(sum(Tdense, dims=1))
    @test mapreduce(one, min, T, dims=1) == mapreduce(one, min, Tdense, dims=1)
    @test mapreduce(one, min, T, dims=2) == mapreduce(one, min, Tdense, dims=2)
    @test mapreduce(one, min, T, dims=3) == mapreduce(one, min, Tdense, dims=3)
    @test typeof(mapreduce(one, min, T, dims=1)) == typeof(mapreduce(one, min, Tdense, dims=1))
    @test mapreduce(zero, max, T, dims=1) == mapreduce(zero, max, Tdense, dims=1)
    @test mapreduce(zero, max, T, dims=2) == mapreduce(zero, max, Tdense, dims=2)
    @test mapreduce(zero, max, T, dims=3) == mapreduce(zero, max, Tdense, dims=3)
    @test typeof(mapreduce(zero, max, T, dims=1)) == typeof(mapreduce(zero, max, Tdense, dims=1))
    @test_throws ArgumentError sum(S, dims=0)
    @test sum(S, dims=1) == sum(Sdense, dims=1)
    @test sum(S, dims=2) == sum(Sdense, dims=2)
    @test sum(S, dims=3) == sum(Sdense, dims=3)
    @test typeof(sum(S, dims=1)) == typeof(sum(Sdense, dims=1))
    @test mapreduce(one, min, S, dims=1) == mapreduce(one, min, Sdense, dims=1)
    @test mapreduce(one, min, S, dims=2) == mapreduce(one, min, Sdense, dims=2)
    @test mapreduce(one, min, S, dims=3) == mapreduce(one, min, Sdense, dims=3)
    @test typeof(mapreduce(one, min, S, dims=1)) == typeof(mapreduce(one, min, Sdense, dims=1))
    @test mapreduce(zero, max, S, dims=1) == mapreduce(zero, max, Sdense, dims=1)
    @test mapreduce(zero, max, S, dims=2) == mapreduce(zero, max, Sdense, dims=2)
    @test mapreduce(zero, max, S, dims=3) == mapreduce(zero, max, Sdense, dims=3)
    @test typeof(mapreduce(zero, max, S, dims=1)) == typeof(mapreduce(zero, max, Sdense, dims=1))

    T = Tridiagonal(Int[], Int[], Int[])
    Tdense = Matrix(T)
    S = SymTridiagonal(Int[], Int[])
    Sdense = Matrix(S)
    @test sum(T) == 0
    @test sum(S) == 0
    @test_throws ArgumentError sum(T, dims=0)
    @test sum(T, dims=1) == sum(Tdense, dims=1)
    @test sum(T, dims=2) == sum(Tdense, dims=2)
    @test sum(T, dims=3) == sum(Tdense, dims=3)
    @test typeof(sum(T, dims=1)) == typeof(sum(Tdense, dims=1))
    @test_throws ArgumentError sum(S, dims=0)
    @test sum(S, dims=1) == sum(Sdense, dims=1)
    @test sum(S, dims=2) == sum(Sdense, dims=2)
    @test sum(S, dims=3) == sum(Sdense, dims=3)
    @test typeof(sum(S, dims=1)) == typeof(sum(Sdense, dims=1))

    T = Tridiagonal(Int[], Int[2], Int[])
    Tdense = Matrix(T)
    S = SymTridiagonal(Int[2], Int[])
    Sdense = Matrix(S)
    @test sum(T) == 2
    @test sum(S) == 2
    @test_throws ArgumentError sum(T, dims=0)
    @test sum(T, dims=1) == sum(Tdense, dims=1)
    @test sum(T, dims=2) == sum(Tdense, dims=2)
    @test sum(T, dims=3) == sum(Tdense, dims=3)
    @test typeof(sum(T, dims=1)) == typeof(sum(Tdense, dims=1))
    @test_throws ArgumentError sum(S, dims=0)
    @test sum(S, dims=1) == sum(Sdense, dims=1)
    @test sum(S, dims=2) == sum(Sdense, dims=2)
    @test sum(S, dims=3) == sum(Sdense, dims=3)
    @test typeof(sum(S, dims=1)) == typeof(sum(Sdense, dims=1))
end

@testset "Issue #28994 (sum of Tridigonal and UniformScaling)" begin
    dl = [1., 1.]
    d = [-2., -2., -2.]
    T = Tridiagonal(dl, d, dl)
    S = SymTridiagonal(T)

    @test diag(T + 2I) == zero(d)
    @test diag(S + 2I) == zero(d)
end

@testset "convert Tridiagonal to SymTridiagonal error" begin
    du = rand(Float64, 4)
    d  = rand(Float64, 5)
    dl = rand(Float64, 4)
    T = Tridiagonal(dl, d, du)
    @test_throws ArgumentError SymTridiagonal{Float32}(T)
end

# Issue #38765
@testset "Eigendecomposition with different lengths" begin
    # length(A.ev) can be either length(A.dv) or length(A.dv) - 1
    A = SymTridiagonal(fill(1.0, 3), fill(-1.0, 3))
    F = eigen(A)
    A2 = SymTridiagonal(fill(1.0, 3), fill(-1.0, 2))
    F2 = eigen(A2)
    test_approx_eq_modphase(F.vectors, F2.vectors)
    @test F.values ≈ F2.values ≈ eigvals(A) ≈ eigvals(A2)
    @test eigvecs(A) ≈ eigvecs(A2)
    @test eigvecs(A, eigvals(A)[1:1]) ≈ eigvecs(A2, eigvals(A2)[1:1])
end

@testset "non-commutative algebra (#39701)" begin
    for A in (SymTridiagonal(Quaternion.(randn(5), randn(5), randn(5), randn(5)), Quaternion.(randn(4), randn(4), randn(4), randn(4))),
              Tridiagonal(Quaternion.(randn(4), randn(4), randn(4), randn(4)), Quaternion.(randn(5), randn(5), randn(5), randn(5)), Quaternion.(randn(4), randn(4), randn(4), randn(4))))
        c = Quaternion(1,2,3,4)
        @test A * c ≈ Matrix(A) * c
        @test A / c ≈ Matrix(A) / c
        @test c * A ≈ c * Matrix(A)
        @test c \ A ≈ c \ Matrix(A)
    end
end

@testset "adjoint of LDLt" begin
    Sr = SymTridiagonal(randn(5), randn(4))
    Sc = SymTridiagonal(complex.(randn(5)) .+ 1im, complex.(randn(4), randn(4)))
    b = ones(size(Sr, 1))

    F = ldlt(Sr)
    @test F\b == F'\b

    F = ldlt(Sc)
    @test copy(Sc')\b == F'\b
end

@testset "symmetric and hermitian tridiagonals" begin
    A = [im 0; 0 -im]
    @test issymmetric(A)
    @test !ishermitian(A)

    # real
    A = SymTridiagonal(randn(5), randn(4))
    @test issymmetric(A)
    @test ishermitian(A)

    A = Tridiagonal(A.ev, A.dv, A.ev .+ 1)
    @test !issymmetric(A)
    @test !ishermitian(A)

    # complex
    # https://github.com/JuliaLang/julia/pull/41037#discussion_r645524081
    S = SymTridiagonal(randn(5) .+ 0im, randn(5) .+ 0im)
    S.ev[end] = im
    @test issymmetric(S)
    @test ishermitian(S)

    S = SymTridiagonal(randn(5) .+ 1im, randn(4) .+ 1im)
    @test issymmetric(S)
    @test !ishermitian(S)

    S = Tridiagonal(S.ev, S.dv, adjoint.(S.ev))
    @test !issymmetric(S)
    @test !ishermitian(S)

    S = Tridiagonal(S.dl, real.(S.d) .+ 0im, S.du)
    @test !issymmetric(S)
    @test ishermitian(S)
end

isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
using .Main.ImmutableArrays

@testset "Conversion to AbstractArray" begin
    # tests corresponding to #34995
    v1 = ImmutableArray([1, 2])
    v2 = ImmutableArray([3, 4, 5])
    v3 = ImmutableArray([6, 7])
    T = Tridiagonal(v1, v2, v3)
    Tsym = SymTridiagonal(v2, v1)

    @test convert(AbstractArray{Float64}, T)::Tridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == T
    @test convert(AbstractMatrix{Float64}, T)::Tridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == T
    @test convert(AbstractArray{Float64}, Tsym)::SymTridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Tsym
    @test convert(AbstractMatrix{Float64}, Tsym)::SymTridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Tsym
end

@testset "dot(x,A,y) for A::Tridiagonal or SymTridiagonal" begin
    for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)
        x = fill(convert(elty, 1), 0)
        T = Tridiagonal(x, x, x)
        Tsym = SymTridiagonal(x, x)
        @test dot(x, T, x) == 0.0
        @test dot(x, Tsym, x) == 0.0
    end
end

isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "SizedArrays.jl"))
using .Main.SizedArrays
@testset "non-number eltype" begin
    @testset "sum for SymTridiagonal" begin
        dv = [SizedArray{(2,2)}(rand(1:2048,2,2)) for i in 1:10]
        ev = [SizedArray{(2,2)}(rand(1:2048,2,2)) for i in 1:10]
        S = SymTridiagonal(dv, ev)
        Sdense = Matrix(S)
        @test Sdense == collect(S)
        @test sum(S) == sum(Sdense)
        @test sum(S, dims = 1) == sum(Sdense, dims = 1)
        @test sum(S, dims = 2) == sum(Sdense, dims = 2)
    end
    @testset "issymmetric/ishermitian for Tridiagonal" begin
        @test !issymmetric(Tridiagonal([[1 2;3 4]], [[1 2;2 3], [1 2;2 3]], [[1 2;3 4]]))
        @test !issymmetric(Tridiagonal([[1 3;2 4]], [[1 2;3 4], [1 2;3 4]], [[1 2;3 4]]))
        @test issymmetric(Tridiagonal([[1 3;2 4]], [[1 2;2 3], [1 2;2 3]], [[1 2;3 4]]))

        @test ishermitian(Tridiagonal([[1 3;2 4].+im], [[1 2;2 3].+0im, [1 2;2 3].+0im], [[1 2;3 4].-im]))
        @test !ishermitian(Tridiagonal([[1 3;2 4].+im], [[1 2;2 3].+0im, [1 2;2 3].+0im], [[1 2;3 4].+im]))
        @test !ishermitian(Tridiagonal([[1 3;2 4].+im], [[1 2;2 3].+im, [1 2;2 3].+0im], [[1 2;3 4].-im]))
    end
    @testset "== between Tridiagonal and SymTridiagonal" begin
        dv = [SizedArray{(2,2)}([1 2;3 4]) for i in 1:4]
        ev = [SizedArray{(2,2)}([3 4;1 2]) for i in 1:4]
        S = SymTridiagonal(dv, ev)
        Sdense = Matrix(S)
        @test S == Tridiagonal(diag(Sdense, -1), diag(Sdense),  diag(Sdense, 1)) == S
        @test S !== Tridiagonal(diag(Sdense, 1), diag(Sdense),  diag(Sdense, 1)) !== S
    end
end
end # module TestTridiagonal

max_stars_repo_name: BSnelling/julia

max_issues_repo_name: BSnelling/julia

max_forks_repo_name: BSnelling/julia


text:
include("utils.jl")

@testset "Generalized Stiefel" begin
    @testset "Real" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedStiefel(3,2,B)
        x = [1.0 0.0; 0.0 0.5; 0.0 0.0]
        @testset "Basics" begin
            @test repr(M) == "GeneralizedStiefel(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℝ)"
            @test representation_size(M) == (3,2)
            @test manifold_dimension(M) == 3
            @test base_manifold(M) === M
            @test_throws DomainError is_manifold_point(M, [1., 0., 0., 0.],true)
            @test_throws DomainError is_manifold_point(M, 1im*[1.0 0.0; 0.0 1.0; 0.0 0.0],true)
            @test !is_tangent_vector(M, x, [0., 0., 1., 0.])
            @test_throws DomainError is_tangent_vector(M, x, [0., 0., 1., 0.], true)
            @test_throws DomainError is_tangent_vector(M, x, 1 * im * zero_tangent_vector(M,x), true)
        end
        @testset "Embedding and Projection" begin
            y = similar(x)
            z = embed(M,x)
            @test z==x
            embed!(M,y,x)
            @test y==z
            a = [1.0 0.0; 0.0 2.0; 0.0 0.0]
            @test !is_manifold_point(M,a)
            b = similar(a)
            c = project(M,a)
            @test c==x
            project!(M,b,a)
            @test b==x
            X = [ 0.0 0.0; 0.0 0.0; -1.0 1.0]
            Y = similar(X)
            Z = embed(M,x,X)
            embed!(M,Y,x,X)
            @test Y == X
            @test Z == X
        end

        types = [Matrix{Float64}, ]
        TEST_STATIC_SIZED && push!(types, MMatrix{3, 2, Float64})
        X = [0.0 0.0; 0.0 0.0; 1.0 1.0]
        Y = [ 0.0 0.0; 0.0 0.0; -1.0 1.0]
        @test inner(M,x,X,Y) == 0
        y = retract(M, x, X)
        z = retract(M, x, Y)
        @test is_manifold_point(M,y)
        @test is_manifold_point(M,z)
        a = project(M, x+X)
        b = retract(M, x, X)
        c = retract(M, x, X, ProjectionRetraction())
        d = retract(M, x, X, PolarRetraction())
        @test a == b
        @test c == d
        @test b == c
        e = similar(a)
        retract!(M, e, x, X)
        @test e == a
        @test vector_transport_to(M, x, X, y, ProjectionTransport()) ==
              project(M, y, X)
        @testset "Type $T" for T in types
            pts = convert.(T, [x, y, z])
            @test !is_manifold_point(M, 2 * x)
            @test_throws DomainError !is_manifold_point(M, 2 * x, true)
            @test !is_tangent_vector(M, x, y)
            @test_throws DomainError is_tangent_vector(M, x, y, true)
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = nothing,
                default_retraction_method = ProjectionRetraction(),
                test_injectivity_radius = false,
                test_is_tangent = true,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                projection_atol_multiplier = 15.0,
                retraction_atol_multiplier = 10.0,
                is_tangent_atol_multiplier = 4 * 10.0^2,
                retraction_methods = [PolarRetraction(), ProjectionRetraction()],
            )
        end
    end

    @testset "Complex" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedStiefel(3, 2, B, ℂ)
        @testset "Basics" begin
            @test repr(M) == "GeneralizedStiefel(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℂ)"
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 8
            @test !is_manifold_point(M, [1.0, 0.0, 0.0, 0.0])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0.0, 0.0, 1.0, 0.0])
            x = [1.0 0.0; 0.0 0.5; 0.0 0.0]

            x = [1im 0.0; 0.0 0.5im; 0.0 0.0]
            @test is_manifold_point(M, x)
            @test !is_manifold_point(M, 2 * x)
        end
    end

    @testset "Quaternion" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedStiefel(3, 2, B, ℍ)
        @test repr(M) == "GeneralizedStiefel(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℍ)"
        @testset "Basics" begin
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 18
        end
    end
end


text:
"""
    area(poly)

Calculate the signed area via the shoelace formula. If the points are ordered
counterclockwise, the result will be positive. If the points are ordered clockwise,
the result will be negative.
"""
function area(poly::T) where T
    a = zero(eltype(eltype(T)))
    @inbounds for i in 1:length(poly)-1
        p1 = poly[i]
        p2 = poly[i+1]
        a += p1[1]*p2[2]-p2[1]*p1[2]
    end
    return a/2
end

function centroid(poly::T) where T
    ZT = zero(eltype(eltype(T)))
    a = ZT
    c = zero(eltype(poly))

    @inbounds for i in 1:length(poly)-1
        p = poly[i]
        n = poly[i+1]
        k = p[1]*n[2]-n[1]*p[2]
        a += k
        c = c .+ (p.+n).*k
    end
    a /= 2
    return c./(6*a)
end


text:
# Autogenerated wrapper script for GAP_pkg_datastructures_jll for aarch64-apple-darwin
export datastructures

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_datastructures")
JLLWrappers.@declare_file_product(datastructures)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        datastructures,
        "lib/gap/datastructures.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export subgraph_density, subgraph_relative_density, common_subgraph_density, dcs_intersection_graph, dcs_union_graph, dcs_avg_graph
export in_degree, out_degree, degree

function getBytes(x)
  total = 0;
  fieldNames = typeof(x).names;
  if fieldNames == ()
    return sizeof(x);
  else
    for fieldName in fieldNames
      try
        total += getBytes(getfield(x,fieldName));
        println("getBytes()$fieldName, $total")
      catch e
        continue
      end
    end
    return total;
  end
end

function in_degree{T}(G::SparseMatrixCSC{T, Int64}, S::Array{Int64, 1}=Int64[])
  n = size(G, 1)
  if length(S) == 0
    S = 1:n
  end
  ns = length(S)
  result = zeros(Float64, ns , 1)
  H = G[S, S]
  (rr, cc) = findn(H)
  ne = length(rr)
  for i=1:ne
    crow = rr[i]
    ccol = cc[i]
    result[ccol] = result[ccol] + H[crow, ccol]
  end
  return result
end

function out_degree{T}(G::SparseMatrixCSC{T, Int64}, S::Array{Int64, 1}=Int64[])
  n = size(G, 1)
  if length(S) == 0
    S = 1:n
  end
  ns = length(S)
  result = zeros(Float64, ns , 1)
  H = G[S, S]
  (rr, cc) = findn(H)
  ne = length(rr)
  for i=1:ne
    crow = rr[i]
    ccol = cc[i]
    result[crow] = result[crow] + H[crow, ccol]
  end
  return result
end

function degree{T}(G::SparseMatrixCSC{T, Int64}, S::Array{Int64, 1}=Int64[]; is_directed=false, direction="both" )
  n = size(G, 1)
  if length(S) == 0
    S = 1:n
  end
  ns = length(S)
  result = zeros(Float64, ns , 1)
  H = G[S, S]
  (rr, cc) = findn(H)
  ne = length(rr)
  @inbounds for i=1:ne
    crow = rr[i]
    ccol = cc[i]
    cval = H[crow, ccol]
    if is_directed==false
      result[crow] = result[crow] + cval
    elseif direction == "both"
      result[crow] = result[crow] + cval
      result[ccol] = result[ccol] + cval
    elseif direction == "in"
      result[ccol] = result[ccol] + cval
    elseif direction == "out"
      result[crow] = result[crow] + cval
    else
      error("Invalid direction option!")
    end
  end
  return result
end

function subgraph_density{T}(G::SparseMatrixCSC{T, Int64}, S::Array{Int64, 1}=Int64[]; is_directed=false)
  n = size(G, 1)
  ns = length(S)
  result = 0.0::Float64
  if ns == 0 # return the graph density
    result= sum(G)/n
  else
    @assert minimum(S) >= 1
    @assert maximum(S) <= n
    result = sum(G[S, S])/ length(S)
  end

  if is_directed == false
    result /= 2.0
  end
  return result
end

function subgraph_relative_density{T}(G::SparseMatrixCSC{T, Int64}, S::Array{Int64, 1}=Int64[]; is_directed=false)
  relative=0.0
  deg = subgraph_density(G, S; is_directed=is_directed)
  n = size(G, 1)
  ns = length(S)
  if ns == 0
    relative = deg/(n-1)
  else
    relative = deg/(ns - 1)
  end

  if is_directed == false
    relative = relative * 2.0
  end

  return relative
end

function common_subgraph_density{T}(G::Array{SparseMatrixCSC{T, Int64}, 1}, S::Array{Int64, 1}; is_directed=false)
  min_density = 0.0
  m = length(G)
  densities = Float64[]
  for j=1:m
    curr_density = subgraph_density(G[j], S; is_directed=is_directed)
    push!(densities, curr_density)
  end
  (min_density, min_idx) = findmin(densities)
  return min_density
end

function dcs_intersection_graph{T}(H::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(H)
  n = size(H[1], 1)
  G = spzeros(T, n, n)

  ii, jj = findn(H[1])
  for k=1:length(ii)
    cr = ii[k]
    cc = jj[k]
    G[cr, cc] = H[1][cr, cc]
    for j=2:m
      G[cr, cc] = min(G[cr, cc], H[j][cr, cc])
    end
  end
  return G
end

function dcs_union_graph{T}(H::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(H)
  n = size(H[1], 1)
  G = spzeros(T, n, n)

  ii, jj = findn(H[1])
  for k=1:length(ii)
    cr = ii[k]
    cc = jj[k]
    G[cr, cc] = H[1][cr, cc]
    for j=2:m
      G[cr, cc] = max(G[cr, cc], H[j][cr, cc])
    end
  end
  return G
end

function dcs_avg_graph{T}(H::Array{SparseMatrixCSC{T, Int64}, 1})
  m = length(H)
  n = size(H[1], 1)
  G = spzeros(Float64, n, n)

  for j=1:m
    ii, jj= findn(H[j])
    for k=1:length(ii)
      cr = ii[k]
      cc = jj[k]
      G[cr, cc] += H[j][cr, cc]/m
    end
  end
  return G
end


text:
function model_filtration(
    conv::Function = expconv;
    t::AbstractVector,
    params::NamedTuple,
    ca::AbstractVector,
)
    if all(haskey(params, key) for key in (:Tp, :Te))
        @extract (ve, vp) params
        fp = vp / params.Tp
        ps = ve / params.Te
    else
        @extract (fp, ps, ve, vp) params
    end
    ct = _model_filtration(conv, t, [fp, ps, ve, vp], ca)
    return ct
end

function _model_filtration(conv, t, params, ca)
    fp, ps, ve, vp = params
    Tminus = vp / fp
    Tplus = ve / ps
    T = (vp + ve) / fp
    ct = ((T - Tminus) / (Tplus - Tminus)) .* expconv(ca, 1 / Tplus, t)
    ct .+= ((Tplus - T) / (Tplus - Tminus)) .* expconv(ca, 1 / Tminus, t)
    ct .*= fp
    return ct
end

function fit_filtration_nls(
    conv::Function = expconv;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    fp, ps, vp, ve = (fill(NaN, volume_size...) for _ = 1:4)
    model(x, p) = _model_filtration(conv, x, p, ca)
    lls_est = fit_filtration_lls(; t, ca, ct, mask).est
    init_fp, init_ps, init_ve, init_vp = select(lls_est, (:fp, :ps, :ve, :vp))
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalue = [init_fp[idx], init_ps[idx], init_ve[idx], init_vp[idx]]
        fp[idx], ps[idx], ve[idx], vp[idx] =
            curve_fit(model, t, ct[idx, :], initialvalue).param
    end
    T = @. (vp + ve) / fp
    Tp = @. vp / fp
    Te = @. ve / ps
    @. T[fp==0] = 0
    @. Tp[fp==0] = 0
    @. Te[ps==0] = 0
    return (; est = (; fp, ps, ve, vp, T, Tp, Te))
end

function fit_filtration_lls(;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    α, β, γ, fp = (fill(NaN, volume_size...) for _ = 1:4)

    M = zeros(num_timepoints, 4)
    M[:, 4] .= cumul_integrate(t, ca, TrapezoidalFast())
    M[:, 3] .= cumul_integrate(t, M[:, 4], TrapezoidalFast())
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 2] .= -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        M[:, 1] .= cumul_integrate(t, M[:, 2], TrapezoidalFast())
        α[idx], β[idx], γ[idx], fp[idx] = M \ ct[idx, :]
    end
    root_term = @. β^2 - 4 * α
    @. root_term[root_term<0] = 0
    Tp = @. (β - sqrt(root_term)) / (2 * α)
    Te = @. (β + sqrt(root_term)) / (2 * α)
    T = @. γ / (α * fp)
    vp = @. fp * Tp
    ve = @. fp * (T - Tp)
    ps = @. ve / Te
    return (; est = (; fp, ps, ve, vp, T, Te, Tp))
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>
#
#  This script is mainly for testing the model input/output works correctly
#
using HMMGradients, FiniteStateTransducers
using Random, Statistics, FileIO
using BSON

setup="2a"

include("WFSTs.jl")
include("Models.jl")
include("Utils.jl")
include("conf/$(setup)/feat_conf.jl")
include("conf/$(setup)/model_conf.jl")

# get transition matrix
lexicon, ilexicon = get_lexicon()
H, L = get_HL(lexicon)
a, A = get_aA(H)
Ns = size(A,1)

# init model
modely = get_convnet(Nf,Ns;
                     Nks=Nks,
                     Nhs=Nhs,
                     strides=strides,
                     dilations=dilations,
                     dropout=dropout,
                     fout=fout)

Nt,Nb = rand(500:1000),4
x = zeros(Float32,Nt,Nf,Nb)
t,b = 500,1
x[t,:,b] .= 1.0

y = modely(x)
z = sum(y[:,:,b],dims=2)
Nt2 = ceil(Int,Nt/3)
@assert Nt2 == size(y,1)
println("Setup = $setup")
println("Num of parameters = $(sum(prod.(size.(params(modely)))))")
println("Context bins = $(subsample*sum( (!).(z .≈ z[100]) ))")


text:
export Activationtrain1st!, activationupdate, collectweights, activationpre

using Flux
using Flux: onehotbatch
using Juno
using Base.Iterators
using StatsBase
using Zygote
using Zygote: @nograd
using LinearAlgebra
using CUDA
using Plots

function collectweights(model)
    for (i,cell) in enumerate(model.cells)
        for (j,mixedop) in enumerate(cell.mixedops)
            for (k,op) in enumerate(mixedop.ops)
                for (l,layer) in enumerate(op.op)
                    if typeof(layer) <: Flux.Conv
                        @show (i,j,k,l,norm(layer.weight),size(layer.weight))
                    end
                end
            end
        end
    end
end

connectstrings = vcat([[string(j,"-",i) for j = 1:i-1] for i = 3:6]...)
axlab = vcat([string("n",e) for e in connectstrings],[string("r",e) for e in connectstrings])

edgedict = Dict{String, Int64}()
for (i,strconn) in enumerate(connectstrings)
    edgedict[strconn] = i
end
primdict = Dict{String, Int64}()
for (i,prim) in enumerate(PRIMITIVES)
    primdict[prim] = i
end
function activationupdatema(model)
    update = all_αs(model).*0
    reducecells = Array{Int64, 1}()
    normalcells = Array{Int64, 1}()
    for (i,cell) in enumerate(model.cells)
        if cell.reduction
            push!(reducecells, i)
        else
            push!(normalcells, i)
        end
    end
    for key in keys(model.activations.currentacts)
        acts = model.activations.currentacts[key]
        if size(acts,2) == 3
            acts = acts[:,2:3,:]
        elseif size(acts,2) == 6
            acts = acts[:,[2,3,5,6],:]
        end
        reduce = acts[reducecells,:,:]
        normal = acts[normalcells,:,:]
        statreduce = x -> mean(abs.(x))
        #statreduce = x -> std(x)
        update[edgedict[key[1:3]]][primdict[key[5:length(key)]]] = statreduce(normal)
        update[edgedict[key[1:3]]+14][primdict[key[5:length(key)]]] = statreduce(reduce)
    end
    update
end
function activationupdatesd(model)
    update = all_αs(model).*0
    reducecells = Array{Int64, 1}()
    normalcells = Array{Int64, 1}()
    for (i,cell) in enumerate(model.cells)
        if cell.reduction
            push!(reducecells, i)
        else
            push!(normalcells, i)
        end
    end
    for key in keys(model.activations.currentacts)
        acts = model.activations.currentacts[key]
        if size(acts,2) == 3
            acts = acts[:,2:3,:]
        elseif size(acts,2) == 6
            acts = acts[:,[2,3,5,6],:]
        end
        reduce = acts[reducecells,:,:]
        normal = acts[normalcells,:,:]
        #statreduce = x -> mean(abs.(x))
        statreduce = x -> std(x)
        update[edgedict[key[1:3]]][primdict[key[5:length(key)]]] = statreduce(normal)
        update[edgedict[key[1:3]]+14][primdict[key[5:length(key)]]] = statreduce(reduce)
    end
    update
end

all_αs(model::DARTSModel) = Flux.params([model.normal_αs, model.reduce_αs])
all_ws(model::DARTSModel) = Flux.params([model.stem, model.cells..., model.global_pooling, model.classifier])

i = 1

function activationpre(loss, model, val)
    acts = []
    for val_batch in CuIterator(val)
        val_loss = loss(model, val_batch...)
        foreach(CUDA.unsafe_free!, val_batch)
        CUDA.reclaim()
        GC.gc()
        ac1 = activationupdatesd(model)
        push!(acts, ac1)
        CUDA.reclaim()
    end
    acts
end

function Activationtrain1st!(loss, model, train, val, opt_α, opt_w, baseacts, losses=[0.0,0.0]; cbepoch = () -> (), cbbatch = () -> ())
    local train_loss
    local val_loss
    local gsα
    w = all_ws(model)
    α = all_αs(model)
    fakeg = all_αs(model).*0
    fakea = all_αs(model).*0
    global i
    for (train_batch, val_batch, base_act) in zip(TrainCuIterator(train), TrainCuIterator(val), acts)
        gsw = gradient(w) do
            train_loss = loss(model, train_batch...)
            return train_loss
        end
        losses[1] = train_loss
        foreach(CUDA.unsafe_free!, train_batch)
        Flux.Optimise.update!(opt_w, w, gsw)
        CUDA.reclaim()

        gsα = gradient(α) do
            val_loss = loss(model, val_batch...)
            return val_loss
        end
        losses[2] = val_loss
        foreach(CUDA.unsafe_free!, val_batch)
        CUDA.reclaim()
        GC.gc()
        gs = [gsα[x] for x in α]
        for (f,g) in zip(fakeg, gs)
            f .-= g
        end
        CUDA.reclaim()
        GC.gc()
        curracts = activationupdatesd(model)
        actupdate = [log.(curr./base) for (curr, base) in zip(curracts, baseacts)]
        for (f,a) in zip(fakea, actupdate)
            f .+= a
        end
        CUDA.reclaim()
        cbbatch()
    end

    #ac1 = activationupdatesd(model)
    p1 = heatmap(axlab, PRIMITIVES, hcat((fakea |> cpu)...), xrotation = 90, title = "Activation Standard Deviation")
    #ac2 = activationupdatema(model)
    #p2 = heatmap(axlab, PRIMITIVES, hcat((ac2 |> cpu)...), xrotation = 90, title = "Activation Mean of Absolute Value")
    #gs = [gsα[x] for x in α]
    p3 = heatmap(axlab, PRIMITIVES, hcat((fakeg |> cpu)...), xrotation = 90 , title = "Alpha Gradient Update")
    plot(p1,p3,layout = (2,1),size = (1200, 800))
    savefig(string("test/models/fig_gsvsacts", i, ".png"));

    i += 1
    cbepoch()
end

max_stars_repo_name: SuReLI/DifferentiableNAS

max_issues_repo_name: SuReLI/DifferentiableNAS

max_forks_repo_name: SuReLI/DifferentiableNAS


text:
function gausschebyshev( n::Int32, kind::Int32=1 )
# GAUSS-CHEBYSHEV NODES AND WEIGHTS.

x = (Array(Float64,n), Array(Float64,n))
# Use known explicit formulas. Complexity O(n).
if kind == 1
    # Gauss-ChebyshevT quadrature, i.e., w(x) = 1/sqrt(1-x^2)
    x = (cos((2*[n:-1:1]-1)*pi/2n), pi./n*ones(n))
elseif kind == 2
    # Gauss-ChebyshevU quadrature, i.e., w(x) = sqrt(1-x^2)
    x = (cos([n:-1:1]*pi./(n+1)), pi/(n+1)*sin([n:-1:1]./(n+1)*pi).^2 )
elseif kind == 3
    # Gauss-ChebyshevV quadrature, i.e., w(x) = sqrt((1+x)/(1-x))
    x = (cos(([n:-1:1]-.5)*pi/(n+.5)), 2*pi/(n+.5)*cos(([n:-1:1]-.5)*pi/(2(n+.5))).^2)
elseif kind == 4
    # Gauss-ChebyshevW quadrature, i.e., w(x) = sqrt((1-x)/(1+x))
    x = (cos([n:-1:1]*pi/(n+.5)), 2*pi/(n+.5)*sin([n:-1:1]*pi/(2(n+.5))).^2)
else
   throw(ArgumentError("Chebyshev kind should be 1, 2, 3, or 4"))
end
return x
end