text:
infile = "./src/Day14/input.txt"
input = readlines(infile)

function part1(input)
    sp = split.(input, " = ")
    instructions = first.(sp)
    data = last.(sp)
    ln = maximum(length.(data)) # 36
    mask = "X"^ln
    memory = Dict{Int,Int}()
    for i in 1:length(sp)
        if instructions[i] == "mask"
            mask = data[i]
        else
            pointer = parse(Int, match(r"\d+", instructions[i]).match)
            towrite = applymask(mask, parse(Int, data[i]))
            memory[pointer] = towrite
        end
    end
    return memory
end

function applymask(mask, data)
    # convert data to binary
    d = reverse(digits(data, base = 2, pad = 36))
    for (i,v) in enumerate(mask)
        if v != 'X'
            d[i] = parse(Int, v)
        end
    end
    # convert binary to int
    val = parse(Int, join(d,""), base = 2)
    return val
end

sp = split.(input, " = ")
instructions = first.(sp)
data = last.(sp)

function part2(instructions, data)
    mask = "X"^36
    memory = Dict{Int,Int}()
    for i in 1:length(instructions)
        if instructions[i] == "mask"
            mask = data[i]
        else
            pointer = parse(Int, match(r"\d+", instructions[i]).match)
            allpointers = applymemmask(mask, pointer)
            d = parse(Int, data[i])
            for p in allpointers
                memory[p] = d
            end
        end
    end
    return memory
end

function applymemmask(mask, pointer)
    allpointers = Int[]
    push!(allpointers, 0)
    l = length(mask)
    pbin = reverse(digits(pointer, base = 2, pad = l))
    for (i, v) in enumerate(mask)
        size = 2^(l - i)
        if v == 'X'
            # double up the entries
            added = allpointers .+ size
            push!(allpointers, added...)
        elseif v == '1'
            allpointers .+= size
        elseif v == '0'
            p = pbin[i]
            if p == 1
                allpointers .+= size
            else
                # nothing
            end
        end
    end
    return allpointers
end

applymemmask("000000000000000000000000000000X1001X", 42)

mem2 = part2(instructions, data)
sum(values(mem2))

using BenchmarkTools
@benchmark part2(instructions, data)


text:
"""
    PowerNetwork <: AbstractPowerSystem

This structure contains constant parameters that define the topology and
physics of the power network.

The object `PowerNetwork` uses its own contiguous indexing for the buses.
The indexing is independent from those specified in the Matpower or the
PSSE input file. However, a correspondence between the two indexing
(Input indexing *to* `PowerNetwork` indexing) is stored inside the
attribute `bus_to_indexes`.

## Note
The object `PowerNetwork` is created in the host memory.
Use a `AbstractFormulation` object to move data to the target device.

"""
struct PowerNetwork <: AbstractPowerSystem
    vbus::Vector{Complex{Float64}}
    # Admittance matrix
    Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}
    # Lines
    lines::Branches{Complex{Float64}}
    # Data
    buses::Array{Float64, 2}
    branches::Array{Float64, 2}
    generators::Array{Float64, 2}
    costs::Union{Array{Float64, 2}, Nothing}
    baseMVA::Float64

    nbus::Int64
    ngen::Int64

    bustype::Vector{Int64}
    bus_to_indexes::Dict{Int, Int}
    ref::Vector{Int64}
    pv::Vector{Int64}
    pq::Vector{Int64}

    sbus::Vector{Complex{Float64}}
    sload::Vector{Complex{Float64}}

    function PowerNetwork(data::Dict{String, Array}; remove_lines=Int[])
        # Parsed data indexes
        BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
        LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()

        # retrive required data
        bus = data["bus"]
        gen = data["gen"]
        lines = data["branch"]
        SBASE = data["baseMVA"][1]
        costs = Base.get(data, "cost", nothing)

        bus_id_to_indexes = get_bus_id_to_indexes(data["bus"])

        # Remove specified lines
        lines = get_active_branches(lines, remove_lines)

        # size of the system
        nbus = size(bus, 1)
        ngen = size(gen, 1)

        # obtain V0 from raw data
        vbus = zeros(Complex{Float64}, nbus)
        for i in 1:nbus
            vbus[i] = bus[i, VM]*exp(1im * pi/180 * bus[i, VA])
        end

        # form Y matrix
        topology = makeYbus(bus, lines, SBASE, bus_id_to_indexes)

        branches = Branches{Complex{Float64}}(
            topology.yff, topology.yft, topology.ytf, topology.ytt,
            topology.from_buses, topology.to_buses,
        )

        # bus type indexing
        ref, pv, pq, bustype = bustypeindex(bus, gen, bus_id_to_indexes)

        sbus, sload = assembleSbus(gen, bus, SBASE, bus_id_to_indexes)
        Ybus = topology.ybus

        new(vbus, Ybus, branches, bus, lines, gen, costs, SBASE, nbus, ngen, bustype, bus_id_to_indexes, ref, pv, pq, sbus, sload)
    end
end

function PowerNetwork(datafile::String)
    data = import_dataset(datafile)
    return PowerNetwork(data)
end

# Getters
## Network attributes
get(pf::PowerNetwork, ::NumberOfBuses) = pf.nbus
get(pf::PowerNetwork, ::NumberOfLines) = size(pf.branches, 1)
get(pf::PowerNetwork, ::NumberOfGenerators) = pf.ngen
get(pf::PowerNetwork, ::NumberOfPVBuses) = length(pf.pv)
get(pf::PowerNetwork, ::NumberOfPQBuses) = length(pf.pq)
get(pf::PowerNetwork, ::NumberOfSlackBuses) = length(pf.ref)

## Loads
get(pf::PowerNetwork, ::ActiveLoad) = real.(pf.sload)
get(pf::PowerNetwork, ::ReactiveLoad) = imag.(pf.sload)
function get(pf::PowerNetwork, ::ActivePower)
    GEN_BUS, PG, QG, _ = IndexSet.idx_gen()
    return pf.generators[:, PG] ./ pf.baseMVA
end
function get(pf::PowerNetwork, ::ReactivePower)
    GEN_BUS, PG, QG, _ = IndexSet.idx_gen()
    return pf.generators[:, QG] ./ pf.baseMVA
end

## Indexing
function get(pf::PowerNetwork, ::GeneratorIndexes)
    GEN_BUS = IndexSet.idx_gen()[1]
    gens = pf.generators
    ngens = size(gens, 1)
    # Create array on host memory
    indexing = zeros(Int, ngens)
    # Here, we keep the same ordering as specified in Matpower.
    for i in 1:ngens
        indexing[i] = pf.bus_to_indexes[gens[i, GEN_BUS]]
    end
    return indexing
end
get(pf::PowerNetwork, ::PVIndexes) = pf.pv
get(pf::PowerNetwork, ::PQIndexes) = pf.pq
get(pf::PowerNetwork, ::SlackIndexes) = pf.ref

# Pretty printing
function Base.show(io::IO, pf::PowerNetwork)
    println("Power Network characteristics:")
    @printf("\tBuses: %d. Slack: %d. PV: %d. PQ: %d\n", pf.nbus, length(pf.ref),
            length(pf.pv), length(pf.pq))
    println("\tGenerators: ", pf.ngen, ".")
    # Print system status
    @printf("\t==============================================\n")
    @printf("\tBUS \t TYPE \t VMAG \t VANG \t P \t Q\n")
    @printf("\t==============================================\n")

    for i=1:pf.nbus
        type = pf.bustype[i]
        vmag = abs(pf.vbus[i])
        vang = angle(pf.vbus[i])*(180.0/pi)
        pinj = real(pf.sbus[i])
        qinj = imag(pf.sbus[i])
        @printf("\t%d \t  %d \t %1.3f\t%3.2f\t%3.3f\t%3.3f\n", i,
                type, vmag, vang, pinj, qinj)
    end
end

# Some utils function

function bounds(pf::PowerNetwork, ::Buses, ::VoltageMagnitude)
    BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
    LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()

    bus = pf.buses
    v_min = convert.(Float64, bus[:, VMIN])
    v_max = convert.(Float64, bus[:, VMAX])
    return v_min, v_max
end

function bounds(pf::PowerNetwork, ::Generators, ::ActivePower)
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()

    gens = pf.generators
    baseMVA = pf.baseMVA

    p_min = convert.(Float64, gens[:, PMIN] / baseMVA)
    p_max = convert.(Float64, gens[:, PMAX] / baseMVA)
    return p_min, p_max
end

function bounds(pf::PowerNetwork, ::Generators, ::ReactivePower)
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()

    gens = pf.generators
    baseMVA = pf.baseMVA

    q_min = convert.(Float64, gens[:, QMIN] / baseMVA)
    q_max = convert.(Float64, gens[:, QMAX] / baseMVA)
    return q_min, q_max
end

function bounds(pf::PowerNetwork, ::Lines, ::ActivePower)
    RATE_A = IndexSet.idx_branch()[6]
    n_lines = get(pf, NumberOfLines())
    # Flow min is always set equal to 0
    flow_min = zeros(n_lines)
    flow_max = (pf.branches[:, RATE_A] ./ pf.baseMVA).^2
    # According to the spec, if RATE_A is equal to 0, then the flow
    # is unconstrained.
    unlimited = findall(isequal(0), flow_max)
    flow_max[unlimited] .= Inf
    return (flow_min, flow_max)
end

"""
    get_costs_coefficients(pf::PowerNetwork)

Return coefficients for costs function.

"""
function get_costs_coefficients(pf::PowerNetwork)
    # indexes
    BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, VA, BASE_KV, ZONE, VMAX, VMIN,
    LAM_P, LAM_Q, MU_VMAX, MU_VMIN = IndexSet.idx_bus()
    GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, PC1, PC2, QC1MIN,
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF, MU_PMAG, MU_PMIN, MU_QMAX,
    MU_QMIN = IndexSet.idx_gen()
    MODEL, STARTUP, SHUTDOWN, NCOST, COST = IndexSet.idx_cost()

    ref = pf.ref
    pv = pf.pv
    pq = pf.pq
    b2i = pf.bus_to_indexes

    # Matpower assumes gens are ordered. Generator in row i has its cost on row i
    # of the cost table.
    gens = pf.generators
    baseMVA = pf.baseMVA
    bus = pf.buses
    ngens = size(gens)[1]
    nbus = size(bus)[1]

    # initialize cost
    # store coefficients in a Float64 array, with 4 columns:
    # - 1st column: bus type
    # - 2nd column: constant coefficient c0
    # - 3rd column: coefficient c1
    # - 4th column: coefficient c2
    coefficients = zeros(ngens, 4)
    # If cost is not specified, we return the array coefficients as is
    if isnothing(pf.costs)
        @warn("PowerSystem: cost is not specified in PowerNetwork dataset")
        return coefficients
    end

    cost_data = pf.costs
    # iterate generators and check if pv or ref.
    for i = 1:ngens
        # only 2nd degree polynomial implemented for now.
        @assert cost_data[i, MODEL] == 2
        @assert cost_data[i, NCOST] == 3
        genbus = b2i[gens[i, GEN_BUS]]
        bustype = bus[genbus, BUS_TYPE]

        # polynomial coefficients
        c0 = cost_data[i, COST+2]
        c1 = cost_data[i, COST+1] * baseMVA
        c2 = cost_data[i, COST] * baseMVA^2
        coefficients[i, :] .= (bustype, c0, c1, c2)
    end
    return coefficients
end

ext: jl


text:
using ExtremumPlanes

@testset "Basics (a.k.a. placing a single box)" begin
  c = Container(10, 10, 10)
  @test c.corners == [Point(0, 0, 0)]
  b1 = Box(5, 5, 5)
  place(b, c, 1)
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using InteractiveDynamics
using DynamicalSystems
using GLMakie

using OrdinaryDiffEq: Vern9
ds = Systems.henonheiles()

u0s = (
    [0.0, -0.25, 0.42, 0.0], # chaotic
    [0.0, 0.1, 0.5, 0.0], # quasiperiodic
    [0.0, 0.30266571044921875, 0.4205654433900762, 0.0], # periodic
)

diffeq = (alg = Vern9(), dtmax = 0.05)
idxs = (1, 2, 4)

fig, obs = interactive_evolution(
    ds, u0s; idxs, tail = 2000, diffeq, colors = COLORS[1:3],
)

ax = content(fig[1,1])
ax.xlabel = "q₁"
ax.ylabel = "q₂"
ax.zlabel = "v₂"
# main.figure[AbstractPlotting.Axis][:names, :axisnames] = ("q₁", "q₂", "v₂")

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 15)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# ============================================================
# Abstract Types
# ============================================================

abstract type AbstractPhysicalSpace end
abstract type AbstractStructPhysicalSpace <: AbstractPhysicalSpace end
abstract type AbstractUnstructPhysicalSpace <: AbstractPhysicalSpace end
abstract type AbstractPhysicalSpace1D <: AbstractPhysicalSpace end
abstract type AbstractPhysicalSpace2D <: AbstractPhysicalSpace end
abstract type AbstractPhysicalSpace3D <: AbstractPhysicalSpace end
abstract type AbstractVelocitySpace end
abstract type AbstractVelocitySpace1D <: AbstractVelocitySpace end
abstract type AbstractVelocitySpace2D <: AbstractVelocitySpace end
abstract type AbstractVelocitySpace3D <: AbstractVelocitySpace end
abstract type AbstractSetup end
abstract type AbstractProperty end
abstract type AbstractGas <: AbstractProperty end
abstract type AbstractPlasma <: AbstractProperty end
abstract type AbstractCondition end
abstract type AbstractSolverSet end

abstract type AbstractControlVolume end
abstract type AbstractInterface end
abstract type AbstractControlVolume1D <: AbstractControlVolume end
abstract type AbstractControlVolume2D <: AbstractControlVolume end
abstract type AbstractUnstructControlVolume <: AbstractControlVolume end
abstract type AbstractInterface1D <: AbstractInterface end
abstract type AbstractInterface2D <: AbstractInterface end

abstract type AbstractSolution end
abstract type AbstractFlux end
abstract type AbstractSolution1D <: AbstractSolution end
abstract type AbstractSolution2D <: AbstractSolution end
abstract type AbstractFlux1D <: AbstractFlux end
abstract type AbstractFlux2D <: AbstractFlux end

abstract type AbstractParticle end
abstract type AbstractParticle1D <: AbstractParticle end
abstract type AbstractParticle2D <: AbstractParticle end


text:
module SpectralLearning

using Quadrature
using HCubature
using ApproxFun
using Flux
using ChainRulesCore
using Zygote
using ForwardDiff

include("spectral_layer.jl")

export Spectral, spectralize

include("utils/integrals.jl")

include("parametric_functions/ParametricFunctions.jl")
using .ParametricFunctions

include("kernel_operators/inner_product.jl")
include("kernel_operators/kernel_operators.jl")

export InnerProduct
export AbstractOperator, KernelOperator
export FourierKernel

include("neural_operators/spectral_neural_operator.jl")
include("neural_operators/finite_dim_operator.jl")

export SpectralNeuralOperator
export FiniteDimOperator

end


text:
# Model m4.4s.jl

using Pkg, DrWatson

begin
    using StanQuap
    using StatisticalRethinking
end

begin
    df = CSV.read(sr_datadir("Howell1.csv"), DataFrame; delim=';')
    df = filter(row -> row[:age] >= 18, df);
    mean_weight = mean(df.weight)
    df.weight_c = df.weight .- mean_weight
end;

stan4_4 = "
data {
    int < lower = 1 > N;               // Sample size
    vector[N] height;                  // Outcome
    vector[N] weight_c;                // Predictor

    int N_new;                         // Number of predictions
    vector[N_new] x_new;               // Predict for x_new
}

parameters {
    real alpha;                        // Intercept
    real beta;                         // Slope (regression coefficients)
    real < lower = 0 > sigma;          // Error SD
}

model {
    height ~ normal(alpha + weight_c * beta , sigma);
}

generated quantities {
    vector[N_new] y_tilde;
    for (n in 1:N_new)
        y_tilde[n] = normal_rng(alpha + beta * x_new[n], sigma);
}
";

data = (
    N = length(df.height), :N_new => 5,
    weight_c = df.weight_c, :height => df.height,
    x_new = [-30, -10, 0, +10, +30]
)
init = (alpha = 180.0, beta = 1.0, sigma = 10.0)
q4_4s, m4_4s, o4_4s = stan_quap("m4.2s", stan4_4; data, init);

if !isnothing(q4_4s)
    quap4_4s_df = sample(q4_4s)          # DataFrame with samples
    precis(quap4_4s_df)
    post4_4s = read_samples(m4_4s)
    post4_4s |> display
end

# End of m4.4s.jl


text:
export convert_options!

"""
converts options for the PARSDMM solver to the right float types
"""
function convert_options!(options,TF)

#options.maxit                 = TI(options.maxit)
options.evol_rel_tol          = TF(options.evol_rel_tol)
options.feas_tol              = TF(options.feas_tol)
options.obj_tol               = TF(options.obj_tol )
options.rho_ini               = convert(Vector{TF},options.rho_ini)
#options.rho_update_frequency  = TI(options.rho_update_frequency)
options.gamma_ini             = TF(options.gamma_ini)
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
#transientplot.jl: Template for transient plots
#-------------------------------------------------------------------------------

#==
===============================================================================#
function fnbuild_transient(data)
	_title = "Transient Response"
	plot = cons(:plot, nstrips=data.nstrips, title=_title, legend=false,
		xaxis = set(scale=:lin, label=LBLAX_TIME),
	)
	set(plot, data.plotattr)

	for (i, sig) in enumerate(data[:siglist])
		(x,y) = sig.idpos
		push!(plot,
			cons(:wfrm, sig.wfrm, sig.line, label="", strip=sig.strip),
			cons(:atext, sig.id, x=x, y=y, align=sig.posalign, strip=sig.strip),
		)
	end

	pcoll = push!(cons(:plotcoll, title=""), plot)
	return pcoll
end

#==Register plot builder
===============================================================================#
push!(plot_builders, :tran => EasyPlot.EasyPlotBuilder(fnbuild_transient))

#Last line


text:
abstract type EnergyBC end

using ..TurbulenceClosures

"""
    Insulating() :: EnergyBC

No energy flux across the boundary.
"""
struct Insulating <: EnergyBC end
function atmos_energy_boundary_state!(nf, bc_energy::Insulating, atmos, args...) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::Insulating,
    atmos,
    args...,
) end

"""
    PrescribedTemperature(fn) :: EnergyBC

Prescribe the temperature at the boundary by `fn`, a function with signature
`fn(state, aux, t)` returning the temperature (in K).
"""
struct PrescribedTemperature{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    FT = eltype(aux⁻)
    _T_0::FT = T_0(atmos.param_set)
    _cv_d::FT = cv_d(atmos.param_set)

    T = bc_energy.fn(state⁻, aux⁻, t)
    E_int⁺ = state⁺.ρ * _cv_d * (T - _T_0)
    state⁺.ρe =
        E_int⁺ + state⁺.ρ * gravitational_potential(atmos.orientation, aux⁻)
end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diffusive⁻,
    hyperdiffusive⁻,
    aux⁻,
    state⁺,
    diffusive⁺,
    hyperdiffusive⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    # TODO: figure out a better way...
    ν, D_t, _ = turbulence_tensors(atmos, state⁻, diffusive⁻, aux⁻, t)
    d_h_tot = -D_t .* diffusive⁻.∇h_tot
    nd_h_tot = dot(n⁻, d_h_tot)
    # both sides involve projections of normals, so signs are consistent
    fluxᵀn.ρe += nd_h_tot * state⁻.ρ
end

"""
    PrescribedEnergyFlux(fn) :: EnergyBC

Prescribe the net inward energy flux across the boundary by `fn`, a function
with signature `fn(state, aux, t)`, returning the flux (in W/m^2).
"""
struct PrescribedEnergyFlux{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    args...,
) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diffusive⁻,
    hyperdiffusive⁻,
    aux⁻,
    state⁺,
    diffusive⁺,
    hyperdiffusive⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    # DG normal is defined in the outward direction
    # we want to prescribe the inward flux
    fluxᵀn.ρe -= bc_energy.fn(state⁻, aux⁻, t)
end

"""
    BulkFormulaEnergy(fn) :: EnergyBC

Calculate the net inward energy flux across the boundary. The drag
coefficient is `C_h = fn_C_h(state, aux, t, normu_int_tan)`. The surface
temperature and q_tot are `T, q_tot = fn_T_and_q_tot(state, aux, t)`.

Return the flux (in W m^-2).
"""
struct BulkFormulaEnergy{FNX, FNTM} <: EnergyBC
    fn_C_h::FNX
    fn_T_and_q_tot::FNTM
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::BulkFormulaEnergy,
    atmos,
    args...,
) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::BulkFormulaEnergy,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diffusive⁻,
    hyperdiffusive⁻,
    aux⁻,
    state⁺,
    diffusive⁺,
    hyperdiffusive⁺,
    aux⁺,
    bctype,
    t,
    state_int⁻,
    diffusive_int⁻,
    aux_int⁻,
)

    u_int⁻ = state_int⁻.ρu / state_int⁻.ρ
    u_int⁻_tan = projection_tangential(atmos, aux_int⁻, u_int⁻)
    normu_int⁻_tan = norm(u_int⁻_tan)
    C_h = bc_energy.fn_C_h(state⁻, aux⁻, t, normu_int⁻_tan)
    T, q_tot = bc_energy.fn_T_and_q_tot(state⁻, aux⁻, t)

    # calculate MSE from the states at the boundary and at the interior point
    ts = PhaseEquil_ρTq(atmos.param_set, state⁻.ρ, T, q_tot)
    ts_int = recover_thermo_state(atmos, atmos.moisture, state_int⁻, aux_int⁻)
    e_pot = gravitational_potential(atmos.orientation, aux⁻)
    e_pot_int = gravitational_potential(atmos.orientation, aux_int⁻)
    MSE = moist_static_energy(ts, e_pot)
    MSE_int = moist_static_energy(ts_int, e_pot_int)

    # TODO: use the correct density at the surface
    ρ_avg = average_density(state⁻.ρ, state_int⁻.ρ)
    # NOTE: difference from design docs since normal points outwards
    fluxᵀn.ρe -= C_h * ρ_avg * normu_int⁻_tan * (MSE - MSE_int)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct CubicHermiteBasis{T} <: AbstractBasis{1,4}
    funcs::SP.PolynomialSystem{4,1}
    points::Vector{T}
    function CubicHermiteBasis(
        funcs::Vector{DP.Polynomial{C,R}},
        points::Vector{T}
    ) where {C,R,T}

        @assert length(funcs) == 4
        @assert length(points) == 2
        @assert points[1] < points[2]
        @assert isapprox(funcs[1] + funcs[3],one(R))

        polysystem = SP.PolynomialSystem(funcs)
        new{T}(polysystem,points)
    end
end

function Base.show(io::IO,basis::CubicHermiteBasis{T}) where {T}
    print(io,"1-D CubicHermiteBasis")
end

function order(basis::CubicHermiteBasis{T}) where {T}
    return 3
end

function type_of_interpolation_points(B::CubicHermiteBasis{T}) where {T}
    return T
end

function CubicHermiteBasis(;start = -1.0, stop = 1.0)
    @assert start < stop

    DP.@polyvar x
    funcs = hermite_polynomials(x,start,stop)
    return CubicHermiteBasis(funcs,[start,stop])
end

function (B::CubicHermiteBasis{T})(x) where {T}
    return SP.evaluate(B.funcs,[x])
end

function derivative(B::CubicHermiteBasis{T},x) where {T}
    vals = SP.jacobian(B.funcs,[x])
    return vals
end

function gradient(B::CubicHermiteBasis{T},x) where {T}
    return derivative(B,x)
end

function hessian(B::CubicHermiteBasis{T},x) where {T}
    h = zeros(4,1,1)
    SP.hessian!(h,B.funcs,[x])
    return h[:,1,1]
end


text:
Φ(p) = p[1:2]
Φ⁻¹(p) = [p[1]; p[2]; (-p[1] * exp(-p[1]^2 - p[2]^2)) * 4]
p = rand(3)
θ = rand() * 2pi
v = [cos(θ); sin(θ)]
covector = dΦ(Φ, p)

@test_broken size(Φ⁻¹(covector)) != size(p)
@test_broken typeof(Φ⁻¹(covector)) != typeof(p)
@test_broken typeof(dΦξ(Φ, v, p)) != Float64


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

ISC_TIME_SECONDS_PRECISION = 10000

# Protocol Types (accept_type)
const ptype_batch_send  = 3 # Batch sends, no asynchrony
const ptype_out_of_band = 4 # Batch sends w/ out of band notification
const ptype_lazy_send   = 5 # Deferred packets delivery

# Protocol Version
const PROTOCOL_VERSION13 = 13
const PROTOCOL_VERSION16 = 16

const CNCT_user              = 1
const CNCT_passwd            = 2
const CNCT_host              = 4
const CNCT_group             = 5
const CNCT_user_verification = 6
const CNCT_specific_data     = 7
const CNCT_plugin_name       = 8
const CNCT_login             = 9
const CNCT_plugin_list       = 10
const CNCT_client_crypt      = 11

const DSQL_close = 1
const DSQL_drop = 2

const isc_info_end            = 1
const isc_info_truncated      = 2
const isc_info_error          = 3
const isc_info_data_not_ready = 4
const isc_info_length         = 126
const isc_info_flag_end       = 127

const isc_info_db_id                 = 4
const isc_info_reads                 = 5
const isc_info_writes                = 6
const isc_info_fetches               = 7
const isc_info_marks                 = 8
const isc_info_implementation        = 11
const isc_info_version               = 12
const isc_info_base_level            = 13
const isc_info_page_size             = 14
const isc_info_num_buffers           = 15
const isc_info_limbo                 = 16
const isc_info_current_memory        = 17
const isc_info_max_memory            = 18
const isc_info_window_turns          = 19
const isc_info_license               = 20
const isc_info_allocation            = 21
const isc_info_attachment_id         = 22
const isc_info_read_seq_count        = 23
const isc_info_read_idx_count        = 24
const isc_info_insert_count          = 25
const isc_info_update_count          = 26
const isc_info_delete_count          = 27
const isc_info_backout_count         = 28
const isc_info_purge_count           = 29
const isc_info_expunge_count         = 30
const isc_info_sweep_interval        = 31
const isc_info_ods_version           = 32
const isc_info_ods_minor_version     = 33
const isc_info_no_reserve            = 34
const isc_info_logfile               = 35
const isc_info_cur_logfile_name      = 36
const isc_info_cur_log_part_offset   = 37
const isc_info_num_wal_buffers       = 38
const isc_info_wal_buffer_size       = 39
const isc_info_wal_ckpt_length       = 40
const isc_info_wal_cur_ckpt_interval = 41
const isc_info_wal_prv_ckpt_fname    = 42
const isc_info_wal_prv_ckpt_poffset  = 43
const isc_info_wal_recv_ckpt_fname   = 44
const isc_info_wal_recv_ckpt_poffset = 45
const isc_info_wal_grpc_wait_usecs   = 47
const isc_info_wal_num_io            = 48
const isc_info_wal_avg_io_size       = 49
const isc_info_wal_num_commits       = 50
const isc_info_wal_avg_grpc_size     = 51
const isc_info_forced_writes         = 52
const isc_info_user_names            = 53
const isc_info_page_errors           = 54
const isc_info_record_errors         = 55
const isc_info_bpage_errors          = 56
const isc_info_dpage_errors          = 57
const isc_info_ipage_errors          = 58
const isc_info_ppage_errors          = 59
const isc_info_tpage_errors          = 60
const isc_info_set_page_buffers      = 61
const isc_info_db_sql_dialect        = 62
const isc_info_db_read_only          = 63
const isc_info_db_size_in_pages      = 64
const isc_info_att_charset           = 101
const isc_info_db_class              = 102
const isc_info_firebird_version      = 103
const isc_info_oldest_transaction    = 104
const isc_info_oldest_active         = 105
const isc_info_oldest_snapshot       = 106
const isc_info_next_transaction      = 107
const isc_info_db_provider           = 108
const isc_info_active_transactions   = 109
const isc_info_active_tran_count     = 110
const isc_info_creation_date         = 111
const isc_info_db_file_size          = 112

# isc_info_sql_records items
const isc_info_req_select_count = 13
const isc_info_req_insert_count = 14
const isc_info_req_update_count = 15
const isc_info_req_delete_count = 16

const isc_info_svc_svr_db_info        = 50
const isc_info_svc_get_license        = 51
const isc_info_svc_get_license_mask   = 52
const isc_info_svc_get_config         = 53
const isc_info_svc_version            = 54
const isc_info_svc_server_version     = 55
const isc_info_svc_implementation     = 56
const isc_info_svc_capabilities       = 57
const isc_info_svc_user_dbpath        = 58
const isc_info_svc_get_env            = 59
const isc_info_svc_get_env_lock       = 60
const isc_info_svc_get_env_msg        = 61
const isc_info_svc_line               = 62
const isc_info_svc_to_eof             = 63
const isc_info_svc_timeout            = 64
const isc_info_svc_get_licensed_users = 65
const isc_info_svc_limbo_trans        = 66
const isc_info_svc_running            = 67
const isc_info_svc_get_users          = 68

const isc_tpb_version1         = 1
const isc_tpb_version3         = 3
const isc_tpb_consistency      = 1
const isc_tpb_concurrency      = 2
const isc_tpb_shared           = 3
const isc_tpb_protected        = 4
const isc_tpb_exclusive        = 5
const isc_tpb_wait             = 6
const isc_tpb_nowait           = 7
const isc_tpb_read             = 8
const isc_tpb_write            = 9
const isc_tpb_lock_read        = 10
const isc_tpb_lock_write       = 11
const isc_tpb_verb_time        = 12
const isc_tpb_commit_time      = 13
const isc_tpb_ignore_limbo     = 14
const isc_tpb_read_committed   = 15
const isc_tpb_autocommit       = 16
const isc_tpb_rec_version      = 17
const isc_tpb_no_rec_version   = 18
const isc_tpb_restart_requests = 19
const isc_tpb_no_auto_undo     = 20
const isc_tpb_lock_timeout     = 21

# Service Parameter Block parameter
const isc_spb_version1              = 1
const isc_spb_current_version       = 2
const isc_spb_version               = isc_spb_current_version
const isc_spb_user_name             = 28 # isc_dpb_user_name
const isc_spb_sys_user_name         = 19 # isc_dpb_sys_user_name
const isc_spb_sys_user_name_enc     = 31 # isc_dpb_sys_user_name_enc
const isc_spb_password              = 29 # isc_dpb_password
const isc_spb_password_enc          = 30 # isc_dpb_password_enc
const isc_spb_command_line          = 105
const isc_spb_dbname                = 106
const isc_spb_verbose               = 107
const isc_spb_options               = 108
const isc_spb_address_path          = 109
const isc_spb_process_id            = 110
const isc_spb_trusted_auth          = 111
const isc_spb_process_name          = 112
const isc_spb_trusted_role          = 113
const isc_spb_connect_timeout       = 57 # isc_dpb_connect_timeout
const isc_spb_dummy_packet_interval = 58 # isc_dpb_dummy_packet_interval
const isc_spb_sql_role_name         = 60 # isc_dpb_sql_role_name

# Database Parameter Block Types
const isc_dpb_version1              = 1
const isc_dpb_page_size             = 4
const isc_dpb_num_buffers           = 5
const isc_dpb_force_write           = 24
const isc_dpb_user_name             = 28
const isc_dpb_password              = 29
const isc_dpb_password_enc          = 30
const isc_dpb_lc_ctype              = 48
const isc_dpb_overwrite             = 54
const isc_dpb_connect_timeout       = 57
const isc_dpb_dummy_packet_interval = 58
const isc_dpb_sql_role_name         = 60
const isc_dpb_set_page_buffers      = 61
const isc_dpb_sql_dialect           = 63
const isc_dpb_set_db_charset        = 68
const isc_dpb_process_id            = 71
const isc_dpb_no_db_triggers        = 72
const isc_dpb_trusted_auth          = 73
const isc_dpb_process_name          = 74
const isc_dpb_utf8_filename         = 77
const isc_dpb_specific_auth_data    = 84
const isc_dpb_auth_plugin_list      = 85
const isc_dpb_auth_plugin_name      = 86
const isc_dpb_config                = 87
const isc_dpb_nolinger              = 88
const isc_dpb_reset_icu             = 89
const isc_dpb_map_attach            = 90
const isc_dpb_session_time_zone     = 91

# isc_info_svc_db_stats params
const isc_spb_sts_data_pages      = 0x01
const isc_spb_sts_db_log          = 0x02
const isc_spb_sts_hdr_pages       = 0x04
const isc_spb_sts_idx_pages       = 0x08
const isc_spb_sts_sys_relations   = 0x10
const isc_spb_sts_record_versions = 0x20
const isc_spb_sts_table           = 0x40
const isc_spb_sts_nocreation      = 0x80

# Transaction informatino items
const isc_info_tra_id                 = 4
const isc_info_tra_oldest_interesting = 5
const isc_info_tra_oldest_snapshot    = 6
const isc_info_tra_oldest_active      = 7
const isc_info_tra_isolation          = 8
const isc_info_tra_access             = 9
const isc_info_tra_lock_timeout       = 10

# SQL information items
const isc_info_sql_select        = 4
const isc_info_sql_bind          = 5
const isc_info_sql_num_variables = 6
const isc_info_sql_describe_vars = 7
const isc_info_sql_describe_end  = 8
const isc_info_sql_sqlda_seq     = 9
const isc_info_sql_message_seq   = 10
const isc_info_sql_type          = 11
const isc_info_sql_sub_type      = 12
const isc_info_sql_scale         = 13
const isc_info_sql_length        = 14
const isc_info_sql_null_ind      = 15
const isc_info_sql_field         = 16
const isc_info_sql_relation      = 17
const isc_info_sql_owner         = 18
const isc_info_sql_alias         = 19
const isc_info_sql_sqlda_start   = 20
const isc_info_sql_stmt_type     = 21
const isc_info_sql_get_plan      = 22
const isc_info_sql_records       = 23
const isc_info_sql_batch_fetch   = 24

const isc_info_sql_stmt_select         = 1
const isc_info_sql_stmt_insert         = 2
const isc_info_sql_stmt_update         = 3
const isc_info_sql_stmt_delete         = 4
const isc_info_sql_stmt_ddl            = 5
const isc_info_sql_stmt_get_segment    = 6
const isc_info_sql_stmt_put_segment    = 7
const isc_info_sql_stmt_exec_procedure = 8
const isc_info_sql_stmt_start_trans    = 9
const isc_info_sql_stmt_commit         = 10
const isc_info_sql_stmt_rollback       = 11
const isc_info_sql_stmt_select_for_upd = 12
const isc_info_sql_stmt_set_generator  = 13
const isc_info_sql_stmt_savepoint      = 14

const isc_arg_end         = 0
const isc_arg_gds         = 1
const isc_arg_string      = 2
const isc_arg_cstring     = 3
const isc_arg_number      = 4
const isc_arg_interpreted = 5
const isc_arg_vms         = 6
const isc_arg_unix        = 7
const isc_arg_domain      = 8
const isc_arg_dos         = 9
const isc_arg_mpexl       = 10
const isc_arg_mpexl_ipc   = 11
const isc_arg_next_mach   = 15
const isc_arg_netware     = 16
const isc_arg_win32       = 17
const isc_arg_warning     = 18
const isc_arg_sql_state   = 19

const op_connect            = 1
const op_exit               = 2
const op_accept             = 3
const op_reject             = 4
const op_protocrol          = 5
const op_disconnect         = 6
const op_response           = 9
const op_attach             = 19
const op_create             = 20
const op_detach             = 21
const op_transaction        = 29
const op_commit             = 30
const op_rollback           = 31
const op_open_blob          = 35
const op_get_segment        = 36
const op_put_segment        = 37
const op_close_blob         = 39
const op_info_database      = 40
const op_info_transaction   = 42
const op_batch_segments     = 44
const op_que_events         = 48
const op_cancel_events      = 49
const op_commit_retaining   = 50
const op_event              = 52
const op_connect_request    = 53
const op_aux_connect        = 53
const op_create_blob2       = 57
const op_allocate_statement = 62
const op_execute            = 63
const op_exec_immediate     = 64
const op_fetch              = 65
const op_fetch_response     = 66
const op_free_statement     = 67
const op_prepare_statement  = 68
const op_info_sql           = 70
const op_dummy              = 71
const op_execute2           = 76
const op_sql_response       = 78
const op_drop_database      = 81
const op_service_attach     = 82
const op_service_detach     = 83
const op_service_info       = 84
const op_service_start      = 85
const op_rollback_retaining = 86
# FB3
const op_update_account_info  = 87
const op_authenticate_user    = 88
const op_partial              = 89
const op_trusted_auth         = 90
const op_cancel               = 91
const op_cont_auth            = 92
const op_ping                 = 93
const op_accept_data          = 94
const op_abort_aux_connection = 95
const op_crypt                = 96
const op_crypt_key_callback   = 97
const op_cond_accept          = 98


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletcirclemarker

"""
    ''_leafletcirclemarker(;kwargs...)
    ''_leafletcirclemarker(children::Any;kwargs...)
    ''_leafletcirclemarker(children_maker::Function;kwargs...)

A LeafletCircleMarker component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `center` (required): . center has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `pathOptions` (optional)
- `radius` (Real; optional)
"""
function ''_leafletcirclemarker(; kwargs...)
        available_props = Symbol[:children, :id, :center, :loading_state, :pathOptions, :radius]
        wild_props = Symbol[]
        return Component("''_leafletcirclemarker", "LeafletCircleMarker", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

''_leafletcirclemarker(children::Any; kwargs...) = ''_leafletcirclemarker(;kwargs..., children = children)
''_leafletcirclemarker(children_maker::Function; kwargs...) = ''_leafletcirclemarker(children_maker(); kwargs...)

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
"""
SymbolicMatrix
"""

# Matrix{<:AbstractSymbolic{Tv}}
# SparseCSR{<:AbstractSymbolic{Tv},Ti}
# SparseCSC{<:AbstractSymbolic{Tv},Ti}
# SparseCOO{<:AbstractSymbolic{Tv},Ti}
# SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}

mutable struct SymbolicMatrix{Tv,MatrixT,Ti} <: AbstractMatrixSymbolic{Tv}
    params::Set{Symbol}
    mat #::MaT{<:AbstractNumberSymbolic{Tv},Ti}
    dim
end

# mutable struct SymbolicMatrixExpression{Tv} <: AbstractVectorSymbolic{Tv}
#     params::Set{Symbol}
#     op::Symbol
#     args::Vector{<:AbstractSymbolic{Tv}}
#     dim::Int
#     others
# end

"""
convert
"""

function Base.convert(::Type{<:AbstractSymbolic{Tv}}, v::Matrix{<:AbstractNumberSymbolic{Tv}}) where {Tv<:Number}
    s = union([u.params for u = v]...)
    SymbolicMatrix{Tv,Matrix,Int}(s, v, size(v))
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicMatrix{S,Matrix,Int}) where {T<:Number,S<:Number}
    v = [convert(AbstractNumberSymbolic{T}, x) for x = s.mat]
    SymbolicMatrix{T,Matrix,Int}(s.params, v, length(v))
end

##

function Base.convert(::Type{<:AbstractSymbolic{Tv}}, v::SparseMatrixCSC{<:AbstractNumberSymbolic{Tv},Ti}) where {Tv<:Number,Ti}
    s = union([u.params for u = v.nzval]...)
    SymbolicMatrix{Tv,SparseMatrixCSC,Ti}(s, v, size(v))
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicMatrix{S,SparseMatrixCSC,Ti}) where {T<:Number,S<:Number,Ti}
    f = s.mat
    v = SparseMatrixCSC(f.m, f.n, f.colptr, f.rowval, [convert(AbstractNumberSymbolic{T}, x) for x = f.nzval])
    SymbolicMatrix{T,SparseMatrixCSC,Ti}(s.params, v, size(v))
end

##

for m in [:SparseCSR, :SparseCSC, :SparseCOO]
    @eval function Base.convert(::Type{<:AbstractSymbolic{Tv}}, v::$m{<:AbstractNumberSymbolic{Tv},Ti}) where {Tv<:Number,Ti}
        s = union([u.params for u = v.val]...)
        SymbolicMatrix{Tv,$m,Ti}(s, v, size(v))
    end
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicMatrix{S,SparseCSR,Ti}) where {T<:Number,S<:Number,Ti}
    f = s.mat
    v = SparseCSR(f.m, f.n, [convert(AbstractNumberSymbolic{T}, x) for x = f.val], f.rowptr, f.colind)
    SymbolicMatrix{T,SparseCSR,Ti}(s.params, v, size(v))
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicMatrix{S,SparseCSC,Ti}) where {T<:Number,S<:Number,Ti}
    f = s.mat
    v = SparseCSC(f.m, f.n, [convert(AbstractNumberSymbolic{T}, x) for x = f.val], f.colptr, f.rowind)
    SymbolicMatrix{T,SparseCSC,Ti}(s.params, v, size(v))
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicMatrix{S,SparseCOO,Ti}) where {T<:Number,S<:Number,Ti}
    f = s.mat
    v = SparseCOO(f.m, f.n, [convert(AbstractNumberSymbolic{T}, x) for x = f.val], f.rowind, f.colind)
    SymbolicMatrix{T,SparseCOO,Ti}(s.params, v, size(v))
end

"""
promotion
"""

function Base.promote_rule(::Type{<:AbstractMatrixSymbolic{T}}, ::Type{S}) where {T<:Number,S<:Number}
    AbstractMatrixSymbolic{promote_type(T,S)}
end

function Base.promote_rule(::Type{<:AbstractMatrixSymbolic{T}}, ::Type{<:AbstractMatrixSymbolic{S}}) where {T<:Number,S<:Number}
    AbstractMatrixSymbolic{promote_type(T,S)}
end

"""
IO for structures
"""

function Base.show(io::IO, x::AbstractMatrixSymbolic)
    expr = _toexpr(x)
    Base.show(io, expr)
end

function _toexpr(x::SymbolicMatrix{Tv,MatrixT,Ti}) where {Tv<:Number,MatrixT,Ti}
    Expr(:call, Symbol(MatrixT), size(x.mat)[1], size(x.mat)[2])
end

# function _toexpr(x::SymbolicVectorExpression)
#     args = [_toexpr(e) for e = x.args]
#     Expr(:call, x.op, args...)
# end

"""
seval(f, env, cache)
Return the value for expr f
"""

function seval(s::SymbolicMatrix{Tv,Matrix,Int}, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    [seval(x, env, cache) for x = s.mat]
end

function seval(s::SymbolicMatrix{Tv,SparseMatrixCSC,Ti}, env::SymbolicEnv, cache::SymbolicCache)::SparseMatrixCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, env, cache) for x = f.nzval])
end

function seval(s::SymbolicMatrix{Tv,SparseCSR,Ti}, env::SymbolicEnv, cache::SymbolicCache)::SparseCSR{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSR{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.rowptr, f.colind)
end

function seval(s::SymbolicMatrix{Tv,SparseCSC,Ti}, env::SymbolicEnv, cache::SymbolicCache)::SparseCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSC{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.colptr, f.rowind)
end

function seval(s::SymbolicMatrix{Tv,SparseCOO,Ti}, env::SymbolicEnv, cache::SymbolicCache)::SparseCOO{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCOO{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.rowind, f.colind)
end

"""
seval(f, dvar, env, cache)
Return the first derivative of expr f with respect to dvar
"""

function seval(s::SymbolicMatrix{Tv,Matrix,Int}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    [seval(x, dvar, env, cache) for x = s.mat]
end

function seval(s::SymbolicMatrix{Tv,SparseMatrixCSC,Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::SparseMatrixCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, dvar, env, cache) for x = f.nzval])
end

function seval(s::SymbolicMatrix{Tv,SparseCSR,Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::SparseCSR{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSR{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowptr, f.colind)
end

function seval(s::SymbolicMatrix{Tv,SparseCSC,Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::SparseCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSC{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.colptr, f.rowind)
end

function seval(s::SymbolicMatrix{Tv,SparseCOO,Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::SparseCOO{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCOO{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowind, f.colind)
end

"""
seval(f, dvar, env, cache)
Return the second derivative of expr f with respect to dvar1 and dvar2
"""

function seval(s::SymbolicMatrix{Tv,Matrix,Int}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    [seval(x, dvar, env, cache) for x = s.mat]
end

function seval(s::SymbolicMatrix{Tv,SparseMatrixCSC,Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::SparseMatrixCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, dvar, env, cache) for x = f.nzval])
end

function seval(s::SymbolicMatrix{Tv,SparseCSR,Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::SparseCSR{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSR{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowptr, f.colind)
end

function seval(s::SymbolicMatrix{Tv,SparseCSC,Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::SparseCSC{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCSC{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.colptr, f.rowind)
end

function seval(s::SymbolicMatrix{Tv,SparseCOO,Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::SparseCOO{Tv,Ti} where {Tv <: Number, Ti}
    f = s.mat
    SparseCOO{Tv,Ti}(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowind, f.colind)
end

"""
seval for vec
"""

function seval(v::Matrix{<:AbstractNumberSymbolic{Tv}}, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    seval(convert(SymbolicMatrix{Tv}, v), env, cache)
end

for m in [:SparseCSR, :SparseCSC, :SparseCOO, :SparseMatrixCSC]
    @eval function seval(v::$m{<:AbstractNumberSymbolic{Tv},Ti}, env::SymbolicEnv, cache::SymbolicCache)::$m{Tv,Ti} where {Tv <: Number,Ti}
        seval(convert(SymbolicMatrix{Tv}, v), env, cache)
    end
end

##

function seval(v::Matrix{<:AbstractNumberSymbolic{Tv}}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    seval(convert(SymbolicMatrix{Tv}, v), dvar, env, cache)
end

for m in [:SparseCSR, :SparseCSC, :SparseCOO, :SparseMatrixCSC]
    @eval function seval(v::$m{<:AbstractNumberSymbolic{Tv},Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::$m{Tv,Ti} where {Tv <: Number,Ti}
        seval(convert(SymbolicMatrix{Tv}, v), dvar, env, cache)
    end
end

##

function seval(v::Matrix{<:AbstractNumberSymbolic{Tv}}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Matrix{Tv} where {Tv <: Number}
    seval(convert(SymbolicMatrix{Tv}, v), dvar, env, cache)
end

for m in [:SparseCSR, :SparseCSC, :SparseCOO, :SparseMatrixCSC]
    @eval function seval(v::$m{<:AbstractNumberSymbolic{Tv},Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::$m{Tv,Ti} where {Tv <: Number,Ti}
        seval(convert(SymbolicMatrix{Tv}, v), dvar, env, cache)
    end
end

# function seval(f::AbstractVectorSymbolic{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
#     get(cache, f) do
#         retval = _eval(Val(f.op), f, env, cache)
#         cache[f] = retval
#     end
# end

# """
# convert
# """

# function Base.convert(::Type{<:AbstractSymbolic{T}}, m::Matrix{<:AbstractSymbolic{S}}) where {T<:Number,S<:Number}
#     [convert(AbstractSymbolic{T}, x) for x = m]
# end

# function Base.convert(::Type{<:AbstractSymbolic{T}}, m::SparseCSR{<:AbstractSymbolic{S},Ti}) where {T<:Number,S<:Number,Ti}
#     SparseCSR{T,Ti}(f.m, f.n, [convert(AbstractSymbolic{T}, x) for x = m.val], f.rowptr, f.colind)
# end

# function Base.convert(::Type{<:AbstractSymbolic{T}}, m::SparseCSC{<:AbstractSymbolic{S},Ti}) where {T<:Number,S<:Number,Ti}
#     SparseCSC{T,Ti}(f.m, f.n, [convert(AbstractSymbolic{T}, x) for x = m.val], f.rowptr, f.colind)
# end

# function Base.convert(::Type{<:AbstractSymbolic{T}}, m::SparseCOO{<:AbstractSymbolic{S},Ti}) where {T<:Number,S<:Number,Ti}
#     SparseCOO{T,Ti}(f.m, f.n, [convert(AbstractSymbolic{T}, x) for x = m.val], f.rowptr, f.colind)
# end

# function Base.convert(::Type{<:AbstractSymbolic{T}}, m::SparseMatrixCSC{<:AbstractSymbolic{S},Ti}) where {T<:Number,S<:Number,Ti}
#     SparseMatrixCSC{T,Ti}(f.m, f.n, f.colptr, f.rowval, [convert(AbstractSymbolic{T}, x) for x = m.nzval])
# end

# """
# seval(f, env, cache)
# Return the value for expr f
# """

# function seval(f::Matrix{<:AbstractSymbolic{Tv}}, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     Tv[seval(x, env, cache) for x = f]
# end

# function seval(f::SparseCSR{<:AbstractSymbolic{Tv},Ti}, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseCSR{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.rowptr, f.colind)
# end

# function seval(f::SparseCSC{<:AbstractSymbolic{Tv},Ti}, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseCSC{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.colptr, f.rowind)
# end

# function seval(f::SparseCOO{<:AbstractSymbolic{Tv},Ti}, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseCOO{Tv,Ti}(f.m, f.n, Tv[seval(x, env, cache) for x = f.val], f.rowind, f.colind)
# end

# function seval(f::SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, env, cache) for x = f.nzval])
# end

# """
# seval(f, dvar, env, cache)
# Return the first derivative of expr f with respect to dvar
# """

# function seval(f::Matrix{<:AbstractSymbolic{Tv}}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     Tv[seval(x, dvar, env, cache) for x = f]
# end

# function seval(f::SparseCSR{<:AbstractSymbolic{Tv},Int}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCSR(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowptr, f.colind)
# end

# function seval(f::SparseCSC{<:AbstractSymbolic{Tv},Int}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCSC(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.colptr, f.rowind)
# end

# function seval(f::SparseCOO{<:AbstractSymbolic{Tv},Int}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCOO(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowind, f.colind)
# end

# function seval(f::SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, dvar, env, cache) for x = f.nzval])
# end

# """
# seval(f, dvar, env, cache)
# Return the second derivative of expr f with respect to dvar1 and dvar2
# """

# function seval(f::Matrix{<:AbstractSymbolic{Tv}}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     Tv[seval(x, dvar, env, cache) for x = f]
# end

# function seval(f::SparseCSR{<:AbstractSymbolic{Tv},Int}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCSR(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowptr, f.colind)
# end

# function seval(f::SparseCSC{<:AbstractSymbolic{Tv},Int}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCSC(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.colptr, f.rowind)
# end

# function seval(f::SparseCOO{<:AbstractSymbolic{Tv},Int}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
#     SparseCOO(f.m, f.n, Tv[seval(x, dvar, env, cache) for x = f.val], f.rowind, f.colind)
# end

# function seval(f::SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where {Tv,Ti}
#     SparseMatrixCSC{Tv,Ti}(f.m, f.n, f.colptr, f.rowval, Tv[seval(x, dvar, env, cache) for x = f.nzval])
# end


text:
using Spandex

n = 3
a = from_csr(
    n,
    n,
    [1, 3, 6, 8],
    [2, 3, 1, 2, 3, 2, 3],
    [1.0, 1.0, 2.0, 4.0, -2.0, 3.0, 15.0],
)

solver = Spandex.CholeskySolver{Float64}(n, n)
# # by default
# solver.use_permutation = true
# solver.use_normalization = true

solve_sym(solver, a)

x = solve(solver, a, [17.0, 2.89, -3.3])

m = SparseArray{Float64}(n)
m[1] = 7.0
m[2] = -5.0
m[3] = 1.0

u = update!(solver, m, 9.0)
@assert sum(abs.(x - u)) > 0

v = downdate!(solver, m, 9.0)
@assert sum(abs.(x - v)) < 1e-8

println("ok")


text:
module OneDmin

using Polynomials
using LinearAlgebra

include("Local/LocalIterations.jl")

using NLPModels
using SolverTools
using SolverCore
using Logging
using Stopping

include("LineSearch.jl")

include("bracket/Bracket.jl")

include("TR/TR1D.jl")

end


text:
# interfacing with Random standard library

using Random

using GPUArrays

mutable struct RNG <: Random.AbstractRNG
    handle::curandGenerator_t
    ctx::CuContext
    typ::Int

    function RNG(typ=CURAND_RNG_PSEUDO_DEFAULT)
        handle = curandCreateGenerator(typ)
        curandSetStream(handle, stream())   # XXX: duplicate with default_rng
        obj = new(handle, context(), typ)
        finalizer(unsafe_destroy!, obj)
        return obj
    end

    # TODO: this design doesn't work nicely in the presence of streams.
    #       if the user switches streams, a local RNG object won't be updated.
end

function unsafe_destroy!(rng::RNG)
    CUDA.isvalid(rng.ctx) || return
    context!(rng.ctx) do
        curandDestroyGenerator(rng)
    end
end

Base.unsafe_convert(::Type{curandGenerator_t}, rng::RNG) = rng.handle

## seeding

function Random.seed!(rng::RNG, seed=Base.rand(UInt64), offset=0)
    curandSetPseudoRandomGeneratorSeed(rng, seed)
    curandSetGeneratorOffset(rng, offset)
    @check unsafe_curandGenerateSeeds(rng) CURAND_STATUS_PREEXISTING_FAILURE
    return
end

Random.seed!(rng::RNG, ::Nothing) = Random.seed!(rng)

## in-place

# uniform
const UniformType = Union{Type{Float32},Type{Float64},Type{UInt32}}
const UniformArray = DenseCuArray{<:Union{Float32,Float64,UInt32}}
function Random.rand!(rng::RNG, A::DenseCuArray{UInt32})
    curandGenerate(rng, A, length(A))
    return A
end
function Random.rand!(rng::RNG, A::DenseCuArray{Float32})
    curandGenerateUniform(rng, A, length(A))
    return A
end
function Random.rand!(rng::RNG, A::DenseCuArray{Float64})
    curandGenerateUniformDouble(rng, A, length(A))
    return A
end

# some functions need pow2 lengths: use a padded array and copy back to the original one
function inplace_pow2(A, f)
    len = length(A)
    if len > 1 && ispow2(len)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        B = similar(A, padlen)
        f(B)
        copyto!(A, 1, B, 1, len)
        CUDA.unsafe_free!(B)
    end
    A
end

# normal
const NormalType = Union{Type{Float32},Type{Float64}}
const NormalArray = DenseCuArray{<:Union{Float32,Float64}}
function Random.randn!(rng::RNG, A::DenseCuArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->curandGenerateNormal(rng, B, length(B), mean, stddev))
    return A
end
function Random.randn!(rng::RNG, A::DenseCuArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->curandGenerateNormalDouble(rng, B, length(B), mean, stddev))
    return A
end

# log-normal
const LognormalType = Union{Type{Float32},Type{Float64}}
const LognormalArray = DenseCuArray{<:Union{Float32,Float64}}
function rand_logn!(rng::RNG, A::DenseCuArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->curandGenerateLogNormal(rng, B, length(B), mean, stddev))
    return A
end
function rand_logn!(rng::RNG, A::DenseCuArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->curandGenerateLogNormalDouble(rng, B, length(B), mean, stddev))
    return A
end

# poisson
const PoissonType = Union{Type{Cuint}}
const PoissonArray = DenseCuArray{Cuint}
function rand_poisson!(rng::RNG, A::DenseCuArray{Cuint}; lambda=1)
    curandGeneratePoisson(rng, A, length(A), lambda)
    return A
end

## out of place

# some functions need pow2 lengths: construct a compatible array and return part of it
function outofplace_pow2(shape, ctor, f)
    len = prod(shape)
    if ispow2(len)
        A = ctor(shape)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        A = ctor(padlen)
        f(A)
        B = reshape(A[1:len], shape)
        return B
    end
end

# arrays
Random.rand(rng::RNG, T::UniformType, dims::Dims) =
    Random.rand!(rng, CuArray{T}(undef, dims))
Random.randn(rng::RNG, T::NormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->CuArray{T}(undef, dims), A->randn!(rng, A; kwargs...))
rand_logn(rng::RNG, T::LognormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->CuArray{T}(undef, dims), A->rand_logn!(rng, A; kwargs...))
rand_poisson(rng::RNG, T::PoissonType, dims::Dims; kwargs...) =
    rand_poisson!(rng, CuArray{T}(undef, dims); kwargs...)

# specify default types
Random.rand(rng::RNG, dims::Dims; kwargs...) = Random.rand(rng, Float32, dims; kwargs...)
Random.randn(rng::RNG, dims::Dims; kwargs...) = Random.randn(rng, Float32, dims; kwargs...)
rand_logn(rng::RNG, dims::Dims; kwargs...) = rand_logn(rng, Float32, dims; kwargs...)
rand_poisson(rng::RNG, dims::Dims; kwargs...) = rand_poisson(rng, Cuint, dims; kwargs...)

# support all dimension specifications
Random.rand(rng::RNG, T::UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(rng, T, Dims((dim1, dims...)))
Random.randn(rng::RNG, T::NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_logn(rng::RNG, T::LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_poisson(rng::RNG, T::PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(rng, T, Dims((dim1, dims...)); kwargs...)


text:
let
    x, y = variables(["x", "y"])
    @test (x*y).alpha == [1 1]
    @test (x*y).c == [1]

    @test (1 + 2*x*y).alpha == [0 0; 1 1]
    @test (1 + 2*x*y).c == [1; 2]

    @test ( (x+y)*(x-y) ).alpha == [2 0; 0 2]
    @test ( (x+y)*(x-y) ).c == [1; -1]

    @test ( (1+x+y)^2 * y ).deg == 3
    @test ( (1+x+y)^2 * y ).alpha == [0 1; 1 1; 0 2; 2 1; 1 2; 0 3]
    @test ( (1+x+y)^2 * y ).c == [1; 2; 2; 1; 2; 1]
end

let
    x, y, z = variables(["x", "y", "z"])
    @test Polyopt.monomials(3,[x,y,z]) == [1,x,y,z,x^2,x*y,x*z,y^2,y*z,z^2,x^3,x^2*y,x^2*z,x*y^2,x*y*z,x*z^2,y^3,y^2*z,y*z^2,z^3]
end


text:
using GridWorld

gfn() = Grid(
    6,
    7,
    (x,y)->string(rand(1:9)),
    line_color="red",
    title=Label("Title", color="#666", size=32),
    title_right=Label("Title Right", color="green", size=32),
    header=Label("Header", color="orange", size=32),
    header_right=Label("Header Right")
)
gfn()


text:
@testset "chainrules.jl" begin
    x = rand()
    test_frule(xlogx, x)
    test_rrule(xlogx, x)
    for x in (-x, 0.0, x)
        y = rand()
        test_frule(xlogy, x, y)
        test_rrule(xlogy, x, y)

        for z in (-y, y)
            test_frule(xlog1py, x, z)
            test_rrule(xlog1py, x, z)
        end
    end

    test_frule(logit, x)
    test_rrule(logit, x)

    for x in (-randexp(), randexp())
        test_frule(log1psq, x)
        test_rrule(log1psq, x)
    end

    # test all `Float64` and `Float32` branches of `logistic`
    for x in (-821.4, -23.5, 12.3, 41.2)
        test_frule(logistic, x)
        test_rrule(logistic, x)
    end
    for x in (-123.2f0, -21.4f0, 8.3f0, 21.5f0)
        test_frule(logistic, x; rtol=1f-3, atol=1f-3)
        test_rrule(logistic, x; rtol=1f-3, atol=1f-3)
    end

    # test all branches of `log1pexp`
    for x in (-20.9, 15.4, 41.5)
        test_frule(log1pexp, x)
        test_rrule(log1pexp, x)
    end
    for x in (8.3f0, 12.5f0, 21.2f0)
        test_frule(log1pexp, x; rtol=1f-3, atol=1f-3)
        test_rrule(log1pexp, x; rtol=1f-3, atol=1f-3)
    end

    for x in (-10.2, -3.3, -0.3)
        test_frule(log1mexp, x)
        test_rrule(log1mexp, x)
    end

    for x in (-10.2, -3.3, -0.3, 0.5)
        test_frule(log2mexp, x)
        test_rrule(log2mexp, x)
    end

    # test all branches of `logexpm1`
    for x in (5.2, 21.4, 41.5)
        test_frule(logexpm1, x)
        test_rrule(logexpm1, x)
    end
    for x in (4.3f0, 12.5f0, 21.2f0)
        test_frule(logexpm1, x; rtol=1f-3, atol=1f-3)
        test_rrule(logexpm1, x; rtol=1f-3, atol=1f-3)
    end

    for x in (-randexp(), randexp()), y in (-randexp(), randexp())
        test_frule(logaddexp, x, y)
        test_rrule(logaddexp, x, y)

        test_frule(logsubexp, x, y)
        test_rrule(logsubexp, x, y)
    end

    for x in (randn(10), randn(10, 8)), dims in (:, 1, 1:2, 2)
        dims isa Colon || all(d <= ndims(x) for d in dims) || continue
        test_frule(logsumexp, x; fkwargs=(dims=dims,))
        test_rrule(logsumexp, x; fkwargs=(dims=dims,))
    end

    for x in (randn(10), randn(10, 8))
        test_frule(softmax, x)
        test_rrule(softmax, x)

        for dims in (1, 1:2, 2)
            all(d <= ndims(x) for d in dims) || continue
            test_frule(softmax, x; fkwargs=(dims=dims,))
            test_rrule(softmax, x; fkwargs=(dims=dims,))
        end
    end
end


text:
import Base: +, -, *, /, ==, transpose, zero
import LinearAlgebra: dot, norm, normalize!
import SimilaritySearch: cosine_distance, angle_distance
export BOW, compute_bow, add!

# const BOW = Dict{Symbol,Int}
const BOW = Dict{Symbol,Float64}

"""
    compute_bow(config::TextConfig, text::String, voc::BOW)

Computes a bag of words and the maximum frequency of the bag
"""
function compute_bow(config::TextConfig, text::String, voc::BOW)
    maxfreq = 0.0
    for token in tokenize(config, text)
        sym = Symbol(token)
        m = get(voc, sym, 0.0) + 1.0
        voc[sym] = m
        maxfreq = max(m, maxfreq)
    end

    voc, maxfreq
end

"""
    compute_bow(tokenlist::AbstractVector{Symbol}, voc::BOW)::Tuple{BOW,Float64}

Updates a BOW using the given list of tokens
"""
function compute_bow(tokenlist::AbstractVector{Symbol}, voc::BOW)
    maxfreq = 0.0
    for sym in tokenlist
        m = get(voc, sym, 0.0) + 1.0
        voc[sym] = m
        maxfreq = max(m, maxfreq)
    end

    voc, maxfreq
end

# these are needed to call `compute_bow` for symbol's list but also for simplicity of the API
compute_bow(tokenlist::AbstractVector{Symbol}) = compute_bow(tokenlist, BOW())
compute_bow(config::TextConfig, tokenlist::AbstractVector{Symbol}, bow::BOW) = compute_bow(tokenlist, bow)
compute_bow(config::TextConfig, tokenlist::AbstractVector{Symbol}) = compute_bow(tokenlist, BOW())

"""
    compute_bow(config::TextConfig, text::String)

Compute a bag of words and the maximum frequency of the bag
"""
compute_bow(config::TextConfig, text::String) = compute_bow(config, text, BOW())

"""
    compute_bow(config::TextConfig, arr::AbstractVector{String})

Computes a bag of word and the maximum frequency of the bag; the input is an array of strings that represent a single object
"""
function compute_bow(config::TextConfig, arr::AbstractVector{String})
    D = BOW()
    maxfreq = 0.0
    for text in arr
       _, maxfreq = compute_bow(config, text, D)
    end

    D, maxfreq
end

"""
    normalize!(bow::BOW)

Inplace normalization of `bow`
"""
function normalize!(bow::BOW)
    s = 1.0 / norm(bow)
    for (k, v) in bow
        bow[k] = v * s
    end

    bow
end

function normalize!(matrix::AbstractVector{BOW})
    for bow in matrix
        normalize!(bow)
    end
end

function dot(a::BOW, b::BOW)
    if length(b) < length(a)
        a, b = b, a  # a must be the smallest bow
    end

    s = 0.0
    for (k, v) in a
        w = get(b, k, 0.0)
        s += v * w
    end

    s
end

function norm(a::BOW)::Float64
    s = 0.0
    for w in values(a)
        s += w * w
    end

    sqrt(s)
end

function zero(::Type{BOW})
    BOW()
end

## inplace sum
function add!(a::BOW, b::BOW)
    for (k, w) in b
        if w != 0.0
            a[k] = get(a, k, 0.0) + w
        end
    end

    a
end

## sum
function +(a::BOW, b::BOW)
    if length(a) < length(b)
        a, b = b, a  # a must be the largest bow
    end

    c = copy(a)
    for (k, w) in b
        if w != 0.0
            c[k] = get(c, k, 0.0) + w
        end
    end

    c
end

## definitions for substraction
function -(a::BOW, b::BOW)
    c = copy(a)
    for (k, w) in b
        if w != 0.0
            c[k] = get(c, k, 0.0) - w
        end
    end

    c
end

## definitions for product

function *(a::BOW, b::BOW)
    if length(b) < length(a)
        a, b = b, a  # a must be the smallest bow
    end

    c = copy(a)
    for k in keys(a)
        w = get(b, k, 0.0)
        if w == 0.0
            delete!(c, k)
        else
            c[k] *= w
        end
    end

    c
end

function *(a::BOW, b::F) where F <: Real
    c = copy(a)
    for (k, v) in a
        c[k] = v * b
    end

    c
end

function *(b::F, a::BOW) where F <: Real
    a * b
end

function /(a::BOW, b::F) where F <: Real
    c = copy(a)
    for (k, v) in a
        c[k] = v / b
    end

    c
end

"""
cosine_distance

Computes the cosine_distance between two weighted bags

It supposes that bags are normalized (see `normalize!` function)

"""
function cosine_distance(a::BOW, b::BOW)::Float64
    return 1.0 - dot(a, b)
end

"""
angle_distance

Computes the angle  between two weighted bags

It supposes that all bags are normalized (see `normalize!` function)

"""
function angle_distance(a::BOW, b::BOW)
    d = dot(a, b)

    if d <= -1.0
        return π
    elseif d >= 1.0
        return 0.0
    elseif d == 0  # turn around for zero vectors, in particular for denominator=0
        return π_2
    else
        return acos(d)
    end
end

function cosine(a::BOW, b::BOW)::Float64
    return dot(a, b) # * a.invnorm * b.invnorm # it is already normalized
end


text:
include("../imports.jl")

@testset "ImageClassification" begin
    method = BlockMethod(
        (Image{2}(), Label(1:2)),
        (
            ProjectiveTransforms((16, 16), inferencefactor=8),
            ImagePreprocessing(),
            OneHot()
        )
    )

    testencoding(method.encodings, method.blocks)
    DLPipelines.checkmethod_core(method)
    @test_nowarn methodlossfn(method)
    @test_nowarn methodmodel(method, Models.xresnet18())

    @testset "`encodeinput`" begin
        image = rand(RGB, 32, 48)

        xtrain = encodeinput(method, Training(), image)
        @test size(xtrain) == (16, 16, 3)
        @test eltype(xtrain) == Float32

        xinference = encodeinput(method, Inference(), image)
        @test size(xinference) == (16, 24, 3)
        @test eltype(xinference) == Float32
    end
    @testset "`encodetarget`" begin
        category = 1
        y = encodetarget(method, Training(), category)
        @test y ≈ [1, 0]
        # depends on buffered interface for `BlockMethod`s and `Encoding`s
        #encodetarget!(y, method, Training(), 2)
        #@test y ≈ [0, 1]
    end
    FastAI.checkmethod_plot(method)
end


text:
module InactivePathsExample

include("../src/Jaynes.jl")
using .Jaynes

model0 = @jaynes function bar(x::Float64)
    sleep(5.0)
end

model1 = @jaynes function foo()
    x ~ Normal(0.0, 1.0)
    y ~ Normal(x, 1.0)
    z ~ Normal(0.0, 1.0)
    q ~ Normal(x, 1.0)
    m ~ bar(x)
end

# Actually see here that the logpdf has a dependency structure:
# x -> y
# x -> q
# z
# x -> m
# which means that the choice at z is "inactive" with respect to the grad of logpdf(tr) for sample traces generated from this model.

tr = simulate(model1, ())
display(tr)

# If you want ∇ logpdf wrt z, is the AD smart enough to ignore bar?

end # module


text:
using Test
using DynamicalBilliards
using DynamicalBilliards.Testing

Floats = [Float16, Float32, Float64, BigFloat]

function alltypes(f, args...)
    r = 0.25; x = y = 1.0
    for T ∈ Floats
        bd = billiard_sinai(T(r), T(x), T(y))
        @test eltype(bd) == T

        p = randominside(bd)
        @test eltype(p) == T
        f(p, bd, args...)

        p = randominside(bd, T(1.0))
        @test eltype(p) == T
        f(p, bd, args...)
    end
end

function type_stability(p, bd)
    @testset "$(tag(p))" begin
        T = eltype(bd)
        for o in bd
            t, cp = collision(p, o)
            @test typeof(t) == T
            @test eltype(cp) == T
        end
        i, t, cp = next_collision(p, bd)
        @test typeof(t) == T
        @test eltype(cp) == T

        ct, poss, vels = evolve(p, bd, 5)
        @test eltype(ct) == T
        @test eltype(poss[1]) == eltype(vels[1]) == T
        ct, poss, vels = evolve(p, bd, 5.0)
        @test eltype(ct) == T
        @test eltype(poss[1]) == eltype(vels[1]) == T

        xt, yt, vxt, vyt, tt = timeseries(p, bd, T(10.0); dt = T(0.1))

        for x in (xt, yt, vxt, vyt, tt)
            @test eltype(x) == T
        end

    end
end

billiards_testset("type stability", type_stability; caller = alltypes)


text:
const CSNES = Ptr{Cvoid}
const CSNESType = Cstring

mutable struct SNES{T}
    ptr::CSNES
    comm::MPI.Comm
    opts::Options{T}
    fn!
    fn_vec
    update_jac!
    jac_A
    jac_P
end

scalartype(::SNES{T}) where {T} = T

Base.cconvert(::Type{CSNES}, obj::SNES) = obj.ptr
Base.unsafe_convert(::Type{Ptr{CSNES}}, obj::SNES) =
    convert(Ptr{CSNES}, pointer_from_objref(obj))

Base.eltype(::SNES{T}) where {T} = T

# How to handle Jacobians?
#  - https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESComputeJacobianDefault.html
#  - https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESComputeJacobianDefaultColor.html
#  -

struct SNESFn{T}
end

struct SNESJac{T}
end

#=
function _snesfn(csnes::CSNES, cx::CVec, cfx::CVec, ctx::Ptr{Cvoid})
    snes = unsafe_pointer_to_objref(ctx)
    snes.Feval(cfx, cx)
end

function _snesjac(csnes::CSNES, cx::CVec, cAmat::CMat, cPmat::CMat, ctx::Ptr{Cvoid})
    snes = unsafe_pointer_to_objref(ctx)
    snes.Jeval(cAmat, cPmat, cx)
end
=#

@for_libpetsc begin

    function SNES{$PetscScalar}(comm::MPI.Comm; kwargs...)
        initialize($PetscScalar)
        opts = Options{$PetscScalar}(kwargs...)
        snes = SNES{$PetscScalar}(C_NULL, comm, opts, nothing, nothing, nothing, nothing, nothing)
        @chk ccall((:SNESCreate, $libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CSNES}), comm, snes)

        with(snes.opts) do
            setfromoptions!(snes)
        end

        if comm == MPI.COMM_SELF
            finalizer(destroy, snes)
        end
        return snes
    end

    function (::SNESFn{$PetscScalar})(csnes::CSNES, cx::CVec, cfx::CVec, ctx::Ptr{Cvoid})::$PetscInt
        snes = unsafe_pointer_to_objref(ctx)
        x = unsafe_localarray($PetscScalar, cx; write=false)
        fx = unsafe_localarray($PetscScalar, cfx; read=false)
        snes.fn!(fx, x)
        Base.finalize(x)
        Base.finalize(fx)
        return $PetscInt(0)
    end

    function setfunction!(snes::SNES{$PetscScalar}, fn!, vec::AbstractVec{$PetscScalar})
        ctx = pointer_from_objref(snes)
        fptr = @cfunction(SNESFn{$PetscScalar}(), $PetscInt, (CSNES, CVec, CVec, Ptr{Cvoid}))
        with(snes.opts) do
            @chk ccall((:SNESSetFunction, $libpetsc), PetscErrorCode,
                (CSNES, CVec, Ptr{Cvoid}, Ptr{Cvoid}),
                snes, vec, fptr, ctx)
        end
        snes.fn_vec = vec
        snes.fn! = fn!
        return nothing
    end

    function destroy(snes::SNES{$PetscScalar})
        finalized($PetscScalar) ||
            @chk ccall((:SNESDestroy, $libpetsc), PetscErrorCode, (Ptr{CSNES},), snes)
        return nothing
    end

    function setfromoptions!(snes::SNES{$PetscScalar})
        @chk ccall((:SNESSetFromOptions, $libpetsc), PetscErrorCode, (CSNES,), snes)
    end

    function gettype(snes::SNES{$PetscScalar})
        t_r = Ref{CSNESType}()
        @chk ccall((:SNESGetType, $libpetsc), PetscErrorCode, (CSNES, Ptr{CSNESType}), snes, t_r)
        return unsafe_string(t_r[])
    end

    function view(snes::SNES{$PetscScalar}, viewer::Viewer{$PetscScalar}=ViewerStdout{$PetscScalar}(snes.comm))
        @chk ccall((:SNESView, $libpetsc), PetscErrorCode,
                    (CSNES, CPetscViewer),
                snes, viewer);
        return nothing
    end

    function (::SNESJac{$PetscScalar})(csnes::CSNES, cx::CVec, cA::CMat, cP::CMat, ctx::Ptr{Cvoid})::$PetscInt
        snes = unsafe_pointer_to_objref(ctx)
        @assert snes.ptr == csnes
        @assert snes.jac_A.ptr == cA
        @assert snes.jac_P.ptr == cP
        x = unsafe_localarray($PetscScalar, cx; write=false)
        snes.update_jac!(x, snes.jac_A, snes.jac_P)
        Base.finalize(x)
        return $PetscInt(0)
    end

    function setjacobian!(snes::SNES{$PetscScalar}, update_jac!, A::AbstractMat{$PetscScalar}, P::AbstractMat{$PetscScalar}=A)
        ctx = pointer_from_objref(snes)
        jacptr = @cfunction(SNESJac{$PetscScalar}(), $PetscInt, (CSNES, CVec, CMat, CMat, Ptr{Cvoid}))

        with(snes.opts) do
            @chk ccall((:SNESSetJacobian, $libpetsc), PetscErrorCode,
                (CSNES, CMat, CMat, Ptr{Cvoid}, Ptr{Cvoid}),
                snes, A, P, jacptr, ctx)
        end
        snes.update_jac! = update_jac!
        snes.jac_A = A
        snes.jac_P = P
        return nothing
    end

    function solve!(x::AbstractVec{$PetscScalar}, snes::SNES{$PetscScalar}, b::AbstractVec{$PetscScalar})
        with(snes.opts) do
            @chk ccall((:SNESSolve, $libpetsc), PetscErrorCode,
            (CSNES, CVec, CVec), snes, b, x)
        end
        return x
    end
    function solve!(x::AbstractVec{$PetscScalar}, snes::SNES{$PetscScalar})
        with(snes.opts) do
            @chk ccall((:SNESSolve, $libpetsc), PetscErrorCode,
            (CSNES, CVec, CVec), snes, C_NULL, x)
        end
        return x
    end

end

solve!(x::AbstractVector{T}, snes::SNES{T}) where {T} = parent(solve!(AbstractVec(x), snes))


text:
module JuliaReachTemplatePkg

end  # module

max_stars_repo_name: JuliaReach/JuliaReachTemplatePkg

max_issues_repo_name: JuliaReach/JuliaReachTemplatePkg

max_forks_repo_name: JuliaReach/JuliaReachTemplatePkg


text:
export MILLISECONDS_IN_SECOND, NANOSECONDS_IN_SECOND, NANOSECONDS_IN_MILLISECOND
export setTimeout

const MILLISECONDS_IN_SECOND     = 1_000.0::Float64
const NANOSECONDS_IN_SECOND      = 1_000_000_000.0::Float64
const NANOSECONDS_IN_MILLISECOND = 1_000_000.0::Float64

"""
    setTimeout(f::Function, timeout::Int)

Creates a `Task` which will asynchornously invoke fucntion `f` after specified `timeout` time in milliseconds.

# Arguments
- `f`::Function, function to be invoked asynchronously
- `timeout`::Int, timeout in milliseconds

# Examples

```
using Rocket
using Dates

println("Before: ", Dates.format(now(), "MM:SS"))
setTimeout(1000) do
    println("Inside: ", Dates.format(now(), "MM:SS"))
end
println("Right after: ", Dates.format(now(), "MM:SS"))
;

# Logs
# Before: 20:59
# Right after: 20:59
# Inside: 21:00
```
"""
function setTimeout(f::Function, timeout::Int)
    @async begin
        sleep(timeout / MILLISECONDS_IN_SECOND)
        f()
    end
end

"""
    combined_type(sources)

Returns a Tuple el-type of observable el-types in `sources` argument in the same order
"""
combined_type(sources) = Tuple{ map(source -> subscribable_extract_type(source), sources)... }

"""
    union_type(sources)

Returns a Union el-type of observable el-types in `sources` argument
"""
union_type(sources) = Union{ map(source -> subscribable_extract_type(source), sources)... }


text:
include("test_common.jl")

@testset "contexts" begin
    NodeCall.clear_context()
    delete_context()
    @test !delete_context(nothing)
    ctx1 = current_context()
    ctx2 = new_context(be_current=false)
    node"const f_ctx1 = () => true"
    switch_context(ctx2)
    node"const f_ctx2 = () => 1"
    @test node"f_ctx2()" == 1
    @test_throws NodeError node"f_ctx1()"
    @test list_contexts() == [ctx1, ctx2]
    switch_context(1)
    @test node"f_ctx1()"
    delete_context(ctx1)
    @test_throws NodeError node"f_ctx1()"
    @test node"f_ctx2()" == 1
    ctx = node"require('vm').createContext()"
    switch_context(ctx)
    delete_context()
end


text:
abstract type AbstractWindFarmModel end

"""
WindFarm(windfarm, windresource, windfarmstates)

Struct defining a wind farm

# Arguments
- `turbine_x::Array{Float}(Nturbines)`: contains windturbine x coordinates in the global reference frame
- `turbine_y::Array{Float}(Nturbines)`: contains windturbine y coordinates in the global reference frame
- `turbine_z::Array{Float}(Nturbines)`: contains windturbine base/z coordinates in the global reference frame
- `turbine_definition_ids::Array{Int}(Nturbines)`: contains integers for each wind turbine specifying its definition
- `turbine_definitions::Array{AbstractTurbineDefinition}(Ntypes)`: contains structs defining each wind turbine definition (design) used in the farm
"""
struct WindFarm{AF1,AF2,AF3,AI,AS} <: AbstractWindFarmModel

    # farm design properties
    turbine_x::AF1
    turbine_y::AF2
    turbine_z::AF3
    turbine_definition_ids::AI
    turbine_definitions::AS

end

struct SingleWindFarmState{TI,AF1,AF2,AF3,AF4,AF5,AF6,AF7,AF8,AF9,AI} <: AbstractWindFarmModel

    # farm properties in rotated frame
    id::TI
    turbine_x::AF1
    turbine_y::AF2
    turbine_z::AF3
    turbine_yaw::AF4
    turbine_ct::AF5
    turbine_ai::AF6
    turbine_inflow_velcities::AF7
    turbine_generators_powers::AF8
    turbine_local_ti::AF9
    sorted_turbine_index::AI

end


text:
"""
    struct DemoSection <: Any
    DemoSection(root::String)

Constructs a demo section that holds either nested subsections or demo cards.

# Fields

Besides the root path to the demo section folder `root`, this struct has some other fields:

* `cards`: demo cards found in `root`
* `subsections`: nested subsections found in `root`

# Configuration

You can manage an extra `config.json` file to customize rendering of a demo section.
Supported items are:

* `order`: specify the cards order or subsections order. By default, it's case-insensitive alphabetic order.
* `title`: specify the title of this demo section. By default, it's the folder name of `root`.

The following is an example of `config.json`:

```json
{
    "title": "learn by examples"
    "order": [
        "quickstart.md",
        "array.md"
    ]
}
```

!!! warning

    You can't specify files or foldernames in other folders.

# Examples

The following is the simplest folder structure of a `DemoSection`:

```text
section
├── demo_1.md
├── demo_2.md
├── demo_3.md
├── demo_4.md
├── demo_5.md
├── demo_6.md
├── demo_7.md
└── demo_8.md
```

The following is a typical folder structure of a `DemoSection`:

```text
section
├── config.json
├── part1
│   ├── demo_21.md
│   └── demo_22.md
└── part2
    ├── config.json
    ├── demo_23.md
    └── demo_24.md
```

!!! warning

    A section should only hold either subsections or demo files. A folder that has both subfolders and demo files (e.g., `*.md`) is invalid.

See also: [`MarkdownDemoCard`](@ref DemoCards.MarkdownDemoCard), [`DemoPage`](@ref DemoCards.DemoPage)
"""
struct DemoSection
    root::String
    cards::Vector # Why would `Vector{<:AbstractDemoCard}` fail here?
    subsections::Vector{DemoSection}
    title::String
end

basename(sec::DemoSection) = basename(sec.root)

function DemoSection(root::String)::DemoSection
    root = replace(root, r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
    isdir(root) || throw(ArgumentError("section root should be a valid dir, instead it's $(root)"))

    path = joinpath.(root, filter(x->!startswith(x, "."), readdir(root))) # filter out hidden files
    card_paths = filter(x->isfile(x) && !endswith(x, config_filename), path)
    section_paths = filter(x->isdir(x)&&!(basename(x) in ignored_dirnames), path)

    if isempty(card_paths) && isempty(section_paths)
        throw(ArgumentError("emtpy section folder $(root)"))
    elseif !xor(isempty(card_paths), isempty(section_paths))
        throw(ArgumentError("section folder $(root) should only hold either cards or subsections"))
    end

    # first consturct an incomplete section
    # then load the config and reconstruct a new one
    section = DemoSection(root,
                          map(democard, card_paths),
                          map(DemoSection, section_paths),
                          "")

    ordered_paths = joinpath.(root, load_config(section, "order"))
    if !isempty(section.cards)
        cards = map(democard, ordered_paths)
        subsections = []
    else
        cards = []
        subsections = map(DemoSection, ordered_paths)
    end

    title = load_config(section, "title")
    DemoSection(root, cards, subsections, title)
end

function load_config(sec::DemoSection, key)
    path = joinpath(sec.root, config_filename)
    config = isfile(path) ? JSON.parsefile(path) : Dict()

    if key == "order"
        haskey(config, key) || return get_default_order(sec)

        order = config[key]
        validate_order(order, sec)
        return order
    elseif key == "title"
        get(config, key, basename(sec))
    else
        throw(ArgumentError("Unrecognized key $(key) for DemoSection"))
    end
end

"""return case-insensitive alphabetic order"""
function get_default_order(sec::DemoSection)
    order = isempty(sec.cards) ? basename.(sec.subsections) : basename.(sec.cards)
    sort(order, by = x->lowercase(x))
end


text:
export JITEventListener, GDBRegistrationListener, IntelJITEventListener,
       OProfileJITEventListener, PerfJITEventListener

@checked struct JITEventListener
    ref::API.LLVMJITEventListenerRef
end
Base.unsafe_convert(::Type{API.LLVMJITEventListenerRef}, listener::JITEventListener) = listener.ref

GDBRegistrationListener()  = JITEventListener(LLVM.API.LLVMCreateGDBRegistrationListener())
IntelJITEventListener()    = JITEventListener(LLVM.API.LLVMCreateIntelJITEventListener())
OProfileJITEventListener() = JITEventListener(LLVM.API.LLVMCreateOProfileJITEventListener())
PerfJITEventListener()     = JITEventListener(LLVM.API.LLVMCreatePerfJITEventListener())


text:
using MultiFloats

import StructArrays
import Base: getproperty

# minimal effort to get MultiFloat struct of arrays to work.
# todo: complex.
Base.getproperty(a::MultiFloat, i::Int) = a._limbs[i]
StructArrays.staticschema(::Type{MultiFloat{T,M}}) where {T,M} = NTuple{M,T}
StructArrays.createinstance(::Type{MultiFloat{T,M}}, args::Vararg{T,M}) where {T,M} = MultiFloat{T,M}(values(args))

# better idea is to just define a 3d array.
struct MFArray{T,M,N,TA} <: AbstractArray{MultiFloat{T,M},N}
    A::TA
end

import Base: size, getindex, setindex, view, IndexStyle
using Base.Cartesian: @ntuple, @nexprs
export MFArray

Base.size(A::MFArray) = reverse(Base.tail(reverse(size(A.A))))

Base.IndexStyle(x::MFArray) = IndexStyle(x.A)

@generated function Base.getindex(A::MFArray{T,M,N}, i::Int) where {T,M,N}
    quote
        $(Expr(:meta, :inline, :propagate_inbounds))
        dist = length(A)
        MultiFloat{T,M}(@ntuple $M j -> A.A[i + (j - 1) * dist])
    end
end

@generated function Base.setindex!(A::MFArray{T,M,N}, v::MultiFloat{T,M}, i::Int) where {T,M,N}
    quote
        $(Expr(:meta, :inline, :propagate_inbounds))
        dist = length(A)
        @nexprs $M j -> A.A[i + (j-1) * dist] = v._limbs[j]
        return v
    end
end

@generated function Base.getindex(A::MFArray{T,M,N}, I::Vararg{Int,N}) where {T,M,N}
    quote
        $(Expr(:meta, :inline, :propagate_inbounds))
        MultiFloat{T,M}(@ntuple $M j -> A.A[I..., j])
    end
end

@generated function Base.setindex!(A::MFArray{T,M,N}, v::MultiFloat{T,M}, I::Vararg{Int,N}) where {T,M,N}
    quote
        $(Expr(:meta, :inline, :propagate_inbounds))
        @nexprs $M j -> A.A[I..., j] = v._limbs[j]
        return v
    end
end

function MFArray(A::AbstractArray{MultiFloat{T,M},N}) where {T,M,N}
    A′ = permutedims(reshape(reinterpret(T, A), M, size(A)...), circshift(OneTo(N+1), -1))
    return MFArray{T,M,N,typeof(A′)}(A′)
end

@inline function view(A::MFArray{T,M,N}, I::Vararg{Any,N}) where {T,M,N}
    B = view(A.A, I..., :)
    return MFArray{T,M,ndims(B)-1,typeof(B)}(B)
end

using VectorizationBase: Unroll, VecUnroll

@generated function rmul!(A::MFArray, G::Rotation2{MultiFloat{T,M},MultiFloat{T,M}}) where {T<:Base.IEEEFloat,M}

    W, Wshift = pick_vector_width_shift(T)

    quote
        m = size(A.A, 1)
        i = 1
        pA = stridedpointer(A.A)

        for step = Base.OneTo(m >> $Wshift)

            I₁ = Unroll{3,1,$M,1,$W,zero(UInt)}((i, G.i  , 1))
            I₂ = Unroll{3,1,$M,1,$W,zero(UInt)}((i, G.i + 1, 1))

            a₁ = MultiFloat(vload(pA, I₁).data)
            a₂ = MultiFloat(vload(pA, I₂).data)

            a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
            a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

            vstore!(pA, VecUnroll(a₁′._limbs), I₁)
            vstore!(pA, VecUnroll(a₂′._limbs), I₂)

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        I₁ = Unroll{3,1,$M,1,$W,typemax(UInt)}((i, G.i  , 1))
        I₂ = Unroll{3,1,$M,1,$W,typemax(UInt)}((i, G.i + 1, 1))

        a₁ = MultiFloat(vload(pA, I₁, remainder_mask).data)
        a₂ = MultiFloat(vload(pA, I₂, remainder_mask).data)

        a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
        a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

        vstore!(pA, VecUnroll(a₁′._limbs), I₁, remainder_mask)
        vstore!(pA, VecUnroll(a₂′._limbs), I₂, remainder_mask)

        A
    end
end

@generated function lmul!(G::Rotation2{MultiFloat{T,M},MultiFloat{T,M}}, A::MFArray) where {T<:Base.IEEEFloat,M}

    W, Wshift = pick_vector_width_shift(T)

    quote
        m = size(A.A, 2)
        i = 1
        pA = stridedpointer(A.A)

        for step = Base.OneTo(m >> $Wshift)

            I₁ = Unroll{3,1,$M,2,$W,zero(UInt)}((G.i,     i, 1))
            I₂ = Unroll{3,1,$M,2,$W,zero(UInt)}((G.i + 1, i, 1))

            a₁ = MultiFloat(vload(pA, I₁).data)
            a₂ = MultiFloat(vload(pA, I₂).data)

            a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
            a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

            vstore!(pA, VecUnroll(a₁′._limbs), I₁)
            vstore!(pA, VecUnroll(a₂′._limbs), I₂)

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        I₁ = Unroll{3,1,$M,2,$W,typemax(UInt)}((G.i    , i, 1))
        I₂ = Unroll{3,1,$M,2,$W,typemax(UInt)}((G.i + 1, i, 1))

        a₁ = MultiFloat(vload(pA, I₁, remainder_mask).data)
        a₂ = MultiFloat(vload(pA, I₂, remainder_mask).data)

        a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
        a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

        vstore!(pA, VecUnroll(a₁′._limbs), I₁, remainder_mask)
        vstore!(pA, VecUnroll(a₂′._limbs), I₂, remainder_mask)

        A
    end
end

const MFVector{T,M} = MFArray{T,M,1}

using Base.Cartesian: @ntuple

import LinearAlgebra: axpy!, dot

@generated function LinearAlgebra.axpy!(a::MultiFloat{T,M}, xs::MFVector{T,M}, ys::MFVector{T,M}) where {T,M}

    W, Wshift = pick_vector_width_shift(T)

    quote
        m = size(xs.A, 1)
        i = 1

        px = stridedpointer(xs.A)
        py = stridedpointer(ys.A)

        av = MultiFloat(@ntuple $M j -> Vec{$W,$T}(a._limbs[j]))

        for step = Base.OneTo(m >> $Wshift)
            I = Unroll{2,1,$M,1,$W,zero(UInt)}((i, 1))

            xi = MultiFloat(vload(px, I).data)
            yi = MultiFloat(vload(py, I).data)

            yi′ = xi * av + yi

            vstore!(py, VecUnroll(yi′._limbs), I)

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return ys
        remainder_mask = mask($T, remainder)

        I = Unroll{2,1,$M,1,$W,typemax(UInt)}((i, 1))

        xi = MultiFloat(vload(px, I, remainder_mask).data)
        yi = MultiFloat(vload(py, I, remainder_mask).data)

        yi′ = xi * av + yi

        vstore!(py, VecUnroll(yi′._limbs), I, remainder_mask)

        return ys
    end
end

@generated function LinearAlgebra.dot(xs::MFVector{T,M}, ys::MFVector{T,M}) where {T,M}
    W, Wshift = pick_vector_width_shift(T)

    quote
        m = size(xs.A, 1)
        i = 1

        px = stridedpointer(xs.A)
        py = stridedpointer(ys.A)

        vec_total = zero(MultiFloat{Vec{$W,$T},$M})

        for step = Base.OneTo(m >> $Wshift)
            I = Unroll{2,1,$M,1,$W,zero(UInt)}((i, 1))

            xi = MultiFloat(vload(px, I).data)
            yi = MultiFloat(vload(py, I).data)

            vec_total += xi * yi

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && @goto sum_scalar
        remainder_mask = mask($T, remainder)

        I = Unroll{2,1,$M,1,$W,typemax(UInt)}((i, 1))

        xi = MultiFloat(vload(px, I, remainder_mask).data)
        yi = MultiFloat(vload(py, I, remainder_mask).data)

        vec_total += xi * yi

        @label sum_scalar
        return +(TupleOfMultiFloat(vec_total)...)
    end
end

const BUFFER = Vector{UInt8}(undef, 0)

const Q_WIDTH = 4
const Q_HEIGHT = 32

function A_panel_width(::Type{MultiFloat{T, M}}) where {T, M}
    W = pick_vector_width(T)
    return W * (((VectorizationBase.L₁CACHE.size ÷ sizeof(T) ÷ M - Q_WIDTH) ÷ Q_HEIGHT - Q_WIDTH) ÷ W)
end

function pack!(Q::AbstractMatrix{MultiFloat{Float64,M}}, τs::AbstractVector{MultiFloat{Float64,M}}, A::AbstractMatrix{MultiFloat{Float64,M}}) where {M}
    m, n = size(A)

    MF_T = MultiFloat{Float64, M}
    W = pick_vector_width(Float64)

    @assert size(Q, 1) == m == Q_HEIGHT
    @assert size(Q, 2) == length(τs) == Q_WIDTH
    @assert n % A_panel_width(MF_T) == 0

    A_size = M * Q_HEIGHT * n
    Q_size = M * Q_HEIGHT * Q_WIDTH
    τ_size = M * Q_WIDTH

    a = max(VectorizationBase.REGISTER_SIZE, VectorizationBase.L₁CACHE.linesize)
    resize!(BUFFER, (A_size + Q_size + τ_size) * sizeof(Float64) + 3a)

    p = VectorizationBase.align(convert(Ptr{Float64}, pointer(BUFFER)), a)
    τs′ = unsafe_wrap(Array, p, τ_size, own=false)
    p = VectorizationBase.align(p + sizeof(Float64) * τ_size, a)
    Q′ = unsafe_wrap(Array, p, Q_size, own=false)
    p = VectorizationBase.align(p + sizeof(Float64) * Q_size, a)
    A′ = unsafe_wrap(Array, p, A_size, own=false)

    # Copy τ, Q and A into the buffers.
    # not gonna order
    copyto!(τs′,  reinterpret(Float64, τs))

    # Q is simply column major.
    copyto!(Q′, reinterpret(Float64, Q))

    # And for A we access
    fast_pack!(
        reshape(A′, W, M * Q_HEIGHT, n ÷ W),
        reshape(reinterpret(Float64, A), M * Q_HEIGHT, W, n ÷ W),
        Val{W}(),
        Val{M * Q_HEIGHT}()
    )

    return Q′, τs′, A′
end

function unpack!(A::AbstractMatrix{MultiFloat{Float64,M}}, Abuf::Vector{Float64}) where {M}
    m, n = size(A)

    W = pick_vector_width(Float64)

    @assert size(A, 1) == Q_HEIGHT
    @assert length(A) == length(Abuf) ÷ M

    # And for A we acces
    fast_unpack!(
        reshape(reinterpret(Float64, A), M * Q_HEIGHT, W, n ÷ W),
        reshape(Abuf, W, M * Q_HEIGHT, n ÷ W),
        Val{W}(),
        Val{M * Q_HEIGHT}()
    )

    return A
end

function fast_pack!(B::AbstractArray{Float64, 3}, A::AbstractArray{Float64, 3}, ::Val{W}, ::Val{M}) where {W, M}
    pA = stridedpointer(A)
    pB = stridedpointer(B)

    for k = axes(A, 3)
        i = 1
        while i < M
            v = vload(pA, Unroll{2,1,W,1,W,typemax(UInt)}((i, 1, k)))
            vt = VectorizationBase.vtranspose(v, Val{W}())
            vstore!(pB, vt, Unroll{2,1,W,1,W,typemax(UInt)}((1, i, k)))
            i += W
        end
    end

    return nothing
end

function fast_unpack!(B::AbstractArray{Float64, 3}, A::AbstractArray{Float64, 3}, ::Val{W}, ::Val{M}) where {W, M}
    pA = stridedpointer(A)
    pB = stridedpointer(B)

    for k = axes(A, 3)
        i = 1
        while i < M
            v = vload(pA, Unroll{2,1,W,1,W,typemax(UInt)}((1, i, k)))
            vt = VectorizationBase.vtranspose(v, Val{W}())
            vstore!(pB, vt, Unroll{2,1,W,1,W,typemax(UInt)}((i, 1, k)))
            i += W
        end
    end

    return nothing
end

using Base.Cartesian: @nexprs, @ntuple

@generated function do_qr!(packedQ, packedτ, packedA, ::Type{MultiFloat{Float64,M}}) where {M}
    # minikernel for size.
    # τ is stored in buffer[1:4] (when interpreted as Float64xM)
    # Q is stored in buffer[5:end][:, 1:4] (interpreted as Matrix{Float64xM} of size 34xN)
    # A is stored in buffer[5:end][:, 5:end] (interpreted as Matrix{Float64xM} of size 34xN)

    A_width = A_panel_width(MultiFloat{Float64, M})
    W = pick_vector_width(Float64)

    quote
        pA = stridedpointer(packedA)
        panel = 0
        @inbounds while panel < length(packedA)
            τi = 1
            for q_col = OneTo(4)
                # broadcast τ's limbs to full-width registers.
                @nexprs $M j -> τ_j = Vec{$W}(packedτ[τi + j - 1])

                τ_mf = MultiFloat(@ntuple $M j -> τ_j)

                # Reduce columns [col_offset, col_offset+4)
                aj = panel + $M * $W * (q_col - 1) + 1

                # Loop over the columns in block of size W
                for col_block = OneTo($A_width ÷ $W)

                    qi = (q_col - 1) * $M * Q_HEIGHT + $M * (q_col - 1) + 1

                    start_ai = aj
                    start_qi = qi

                    # Accumulator. reflector[1] == 1, so no multiplication here.
                    @nexprs $M j -> a_j = vload(pA, (MM{$W}(aj + $W * (j - 1)),))

                    accumulator_mf = MultiFloat(@ntuple $M j -> a_j)

                    aj += $M * $W
                    qi += $M

                    # dot over the next 31 rows.
                    for row = OneTo(Q_HEIGHT - q_col)
                        # Load the Q-values and broadcast them to a single register
                        @nexprs $M j -> q_j = Vec{$W}(packedQ[qi + j - 1])
                        q_mf = MultiFloat(@ntuple $M j -> q_j)

                        # Load the values of A of 4 columns
                        @nexprs $M j -> a_j = vload(pA, (MM{$W}(aj + $W * (j - 1)),))
                        a_mf = MultiFloat(@ntuple $M j -> a_j)

                        accumulator_mf += q_mf * a_mf

                        qi += $M
                        aj += $W * $M
                    end

                    α_mf = τ_mf * accumulator_mf

                    # go back to the top.
                    aj = start_ai
                    qi = start_qi

                    @nexprs $M j -> a_j = vload(pA, (MM{$W}(aj + $W * (j - 1)),))

                    a_mf = MultiFloat(@ntuple $M j -> a_j)

                    a_mf -= α_mf

                    @nexprs $M j -> vstore!(pA, a_mf._limbs[j], (MM{$W}(aj + $W * (j - 1)),))

                    aj += $M * $W
                    qi += $M

                    # Finally we do the axpy a -= α * q
                    for row = OneTo(Q_HEIGHT - q_col)
                        # Load the Q-values and broadcast them to a single register
                        @nexprs $M j -> q_j = Vec{$W}(packedQ[qi + j - 1])

                        q_mf = MultiFloat(@ntuple $M j -> q_j)

                        # Load the values of A of 4 columns
                        @nexprs $M j -> a_j = vload(pA, (MM{$W}(aj + $W * (j - 1)),))
                        a_mf = MultiFloat(@ntuple $M j -> a_j)

                        a_mf -= α_mf * q_mf

                        @nexprs $M j -> vstore!(pA, a_mf._limbs[j], (MM{$W}(aj + $W * (j - 1)),))

                        aj += $M * $W
                        qi += $M
                    end

                    aj += $M * $W * (q_col - 1)
                end

                τi += $M
            end

            panel += Q_HEIGHT * $A_width * $M
        end
    end
end