text:
using Oceananigans.Operators

"""
    struct MixingLength{FT}

Contains mixing length parameters for CATKE vertical diffusivity.

The mixing length is the maximum of three mixing lengths:

```math
ℓᵩ = max(ℓᵟᵩ, ℓᴬᵩ, ℓ⋆ᵩ)
```

where ``ℓᵟ`` is a grid-spacing-proportional length,
``ℓʰᵩ`` is a convective mixing length, and ``ℓ⋆ᵩ`` is
a mixing length in stably-stratified mixing.

Grid-spacing-proportional mixing length
=======================================

```math
ℓᵟᵩ = Cᵟϕ * Δz
```

Convective mixing length
========================

```math
ℓᴬᵩ = convecting ? Cᴬϕ * e^3/2 / Qᵇ : 0
```

Stably-stratified mixing length
===============================

```math
ℓ⋆ᵩ = σᵩ * min(d, Cᵇ * √e / N)
```

where ``σᵩ(Ri)`` is a stability function that depends on
the local Richardson number, ...

The Richardson number is

    ``Ri = ∂z B / ( (∂z U)² + (∂z V)² )`` ,

where ``B`` is buoyancy and ``∂z`` denotes a vertical derviative.
The Richardson-number dependent diffusivities are multiplied by the stability
function

    1. ``σ(Ri) = σ⁻ * (1 + rσ * step(Ri, Riᶜ, Riʷ))``

where ``σ₀``, ``Δσ``, ``Riᶜ``, and ``Riʷ`` are free parameters,
and ``step`` is a smooth step function defined by

    ``step(x, c, w) = (1 + \tanh((x - c) / w)) / 2``.

The 8 free parameters in `RiDependentDiffusivityScaling` have been _experimentally_ calibrated
against large eddy simulations of ocean surface boundary layer turbulence in idealized
scenarios involving monotonic boundary layer deepening into variable stratification
due to constant surface momentum fluxes and/or destabilizing surface buoyancy flux.
See https://github.com/CliMA/LESbrary.jl for more information about the large eddy simulations.
The calibration was performed using a combination of Markov Chain Monte Carlo (MCMC)-based simulated
annealing and noisy Ensemble Kalman Inversion methods.
"""
Base.@kwdef struct MixingLength{FT}
    Cᵇ    :: FT = Inf # Inf is the "inert" value for this parameter
    Cˢ    :: FT = Inf # Inf is the "inert" value for this parameter
    Cᵇu   :: FT = 0.596
    Cᵇc   :: FT = 0.0723
    Cᵇe   :: FT = 0.637
    Cˢu   :: FT = 0.628
    Cˢc   :: FT = 0.426
    Cˢe   :: FT = 0.711
    Cᵟu   :: FT = 0.5
    Cᵟc   :: FT = 0.5
    Cᵟe   :: FT = 0.5
    Cᴬu   :: FT = 0.0
    Cᴬc   :: FT = 1.41
    Cᴬe   :: FT = 0.0
    Cᴬˢu  :: FT = 0.0
    Cᴬˢc  :: FT = 0.812
    Cᴬˢe  :: FT = 0.0
    Cᴷu⁻  :: FT = 0.343
    Cᴷuʳ  :: FT = -0.721
    Cᴷc⁻  :: FT = 0.343
    Cᴷcʳ  :: FT = -0.891
    Cᴷe⁻  :: FT = 1.42
    Cᴷeʳ  :: FT = 1.50
    CᴷRiʷ :: FT = 0.101
    CᴷRiᶜ :: FT = 2.03
end

#####
##### Mixing length
#####

@inline surface(i, j, k, grid)                = znode(Center(), Center(), Face(), i, j, grid.Nz+1, grid)
@inline bottom(i, j, k, grid)                 = znode(Center(), Center(), Face(), i, j, 1, grid)
@inline depthᶜᶜᶜ(i, j, k, grid)               = surface(i, j, k, grid) - znode(Center(), Center(), Center(), i, j, k, grid)
@inline height_above_bottomᶜᶜᶜ(i, j, k, grid) = znode(Center(), Center(), Center(), i, j, k, grid) - bottom(i, j, k, grid)

@inline wall_vertical_distanceᶜᶜᶜ(i, j, k, grid) = min(depthᶜᶜᶜ(i, j, k, grid), height_above_bottomᶜᶜᶜ(i, j, k, grid))

@inline function sqrt_∂z_b(i, j, k, grid, buoyancy, tracers)
    FT = eltype(grid)
    N² = ∂z_b(i, j, k, grid, buoyancy, tracers)
    N²⁺ = max(zero(FT), N²)
    return sqrt(N²⁺)
end

@inline function buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, e, tracers, buoyancy)
    FT = eltype(grid)
    N⁺ = ℑzᵃᵃᶜ(i, j, k, grid, sqrt_∂z_b, buoyancy, tracers)
    @inbounds e⁺ = max(zero(FT), e[i, j, k])
    return @inbounds ifelse(N⁺ == 0, FT(Inf), sqrt(e⁺) / N⁺)
end

@inline function shear_mixing_lengthᶜᶜᶜ(i, j, k, grid, e, velocities, tracers, buoyancy)
    FT = eltype(grid)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    S = sqrt(∂z_u² + ∂z_v²)
    @inbounds e⁺ = max(zero(FT), e[i, j, k])
    return ifelse(S == 0, FT(Inf), sqrt(e⁺) / S)
end

# TODO: Use types to distinguish between tracer, velocity, and TKE cases?
@inline function stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᵇ::Number, Cˢ::Number, e, velocities, tracers, buoyancy)
    d = wall_vertical_distanceᶜᶜᶜ(i, j, k, grid)
    ℓᵇ = Cᵇ * buoyancy_mixing_lengthᶜᶜᶜ(i, j, k, grid, e, tracers, buoyancy)
    ℓˢ = Cˢ * shear_mixing_lengthᶜᶜᶜ(i, j, k, grid, e, velocities, tracers, buoyancy)
    return min(d, ℓᵇ, ℓˢ)
end

@inline function convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᴬ::Number, Cᴬˢ::Number,
                                             velocities, tracers, buoyancy, clock, tracer_bcs)
    # Shear
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    S = sqrt(∂z_u² + ∂z_v²)

    # Surface buoyancy flux
    Qᵇ = top_buoyancy_flux(i, j, grid, buoyancy, tracer_bcs, clock, merge(velocities, tracers))

    # Strictly positive TKE
    FT = eltype(grid)
    e⁺ = @inbounds max(zero(FT), tracers.e[i, j, k])

    # "Sheared convection number"
    α = S * Qᵇ / e⁺

    # Mixing length
    ℓᴬ = sqrt(e⁺^3) / Qᵇ
    ℓʰ = Cᴬ * ℓᴬ * (1 - Cᴬˢ * α)

    # Are we convecting?
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    convecting = (N² < 0) & (Qᵇ > 0) & (e⁺ > 0)

    return ifelse(convecting, ℓʰ, zero(FT))
end

@inline ϕ²(i, j, k, grid, ϕ, args...) = ϕ(i, j, k, grid, args...)^2

@inline function Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    FT = eltype(grid)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return ifelse(N² == 0, zero(FT), N² / (∂z_u² + ∂z_v²))
end

@inline function Riᶜᶜᶠ(i, j, k, grid, velocities, tracers, buoyancy)
    FT = eltype(grid)
    ∂z_u² = ℑxᶜᵃᵃ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, velocities.u)
    ∂z_v² = ℑyᵃᶜᵃ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, velocities.v)
    N² = ∂z_b(i, j, k, grid, buoyancy, tracers)
    return ifelse(N² == 0, zero(FT), N² / (∂z_u² + ∂z_v²))
end

@inline step(x, c, w) = (1 + tanh(x / w - c)) / 2
@inline scale(Ri, σ⁻, rσ, c, w) = σ⁻ * (1 + rσ * step(Ri, c, w))

@inline function momentum_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷu⁻,
                 closure.mixing_length.Cᴷuʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function tracer_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷc⁻,
                 closure.mixing_length.Cᴷcʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function TKE_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)
    Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy)
    return scale(Ri,
                 closure.mixing_length.Cᴷe⁻,
                 closure.mixing_length.Cᴷeʳ,
                 closure.mixing_length.CᴷRiᶜ,
                 closure.mixing_length.CᴷRiʷ)
end

@inline function momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬ = closure.mixing_length.Cᴬu
    Cᴬˢ = closure.mixing_length.Cᴬˢu
    ℓʰ = convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᴬ, Cᴬˢ, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟu = closure.mixing_length.Cᵟu
    ℓᵟ = Δzᶜᶜᶜ(i, j, k, grid)

    Cᵇ = min(closure.mixing_length.Cᵇ, closure.mixing_length.Cᵇu)
    Cˢ = min(closure.mixing_length.Cˢ, closure.mixing_length.Cˢu)
    ℓ★ = stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᵇ, Cˢ, tracers.e, velocities, tracers, buoyancy)

    σu = momentum_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)

    return max(ℓʰ, σu * max(Cᵟu * ℓᵟ, ℓ★))
end

@inline function tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬ = closure.mixing_length.Cᴬc
    Cᴬˢ = closure.mixing_length.Cᴬˢc
    ℓʰ = convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᴬ, Cᴬˢ, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟc = closure.mixing_length.Cᵟc
    ℓᵟ = Δzᶜᶜᶜ(i, j, k, grid)

    Cᵇ = min(closure.mixing_length.Cᵇ, closure.mixing_length.Cᵇc)
    Cˢ = min(closure.mixing_length.Cˢ, closure.mixing_length.Cˢc)
    ℓ★ = stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᵇ, Cˢ, tracers.e, velocities, tracers, buoyancy)

    σc = tracer_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)

    return max(ℓʰ, σc * max(Cᵟc * ℓᵟ, ℓ★))
end

@inline function TKE_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴬ = closure.mixing_length.Cᴬe
    Cᴬˢ = closure.mixing_length.Cᴬˢe
    ℓʰ = convective_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᴬ, Cᴬˢ, velocities, tracers, buoyancy, clock, tracer_bcs)

    Cᵟe = closure.mixing_length.Cᵟe
    ℓᵟ = Δzᶜᶜᶜ(i, j, k, grid)

    Cᵇ = min(closure.mixing_length.Cᵇ, closure.mixing_length.Cᵇe)
    Cˢ = min(closure.mixing_length.Cˢ, closure.mixing_length.Cˢe)
    ℓ★ = stable_mixing_lengthᶜᶜᶜ(i, j, k, grid, Cᵇ, Cˢ, tracers.e, velocities, tracers, buoyancy)

    σe = TKE_stable_mixing_scale(i, j, k, grid, closure, velocities, tracers, buoyancy)

    return max(ℓʰ, σe * max(Cᵟe * ℓᵟ, ℓ★))
end

Base.show(io::IO, ML::MixingLength) =
    print(io, "MixingLength: \n" *
              "     Cᵇ   = $(ML.Cᵇ),    \n" *
              "     Cᵇu  = $(ML.Cᵇu),   \n" *
              "     Cᵇc  = $(ML.Cᵇe),   \n" *
              "     Cᵇe  = $(ML.Cᵇc),   \n" *
              "     Cˢ   = $(ML.Cˢ),    \n" *
              "     Cˢu  = $(ML.Cˢu),   \n" *
              "     Cˢc  = $(ML.Cˢe),   \n" *
              "     Cˢe  = $(ML.Cˢc),   \n" *
              "     Cᵟu  = $(ML.Cᵟu),   \n" *
              "     Cᵟc  = $(ML.Cᵟc),   \n" *
              "     Cᵟe  = $(ML.Cᵟe),   \n" *
              "     Cᴬu  = $(ML.Cᴬu),   \n" *
              "     Cᴬc  = $(ML.Cᴬc),   \n" *
              "     Cᴬe  = $(ML.Cᴬe),   \n" *
              "     Cᴷu⁻ = $(ML.Cᴷu⁻),  \n" *
              "     Cᴷc⁻ = $(ML.Cᴷc⁻),  \n" *
              "     Cᴷe⁻ = $(ML.Cᴷe⁻),  \n" *
              "     Cᴷuʳ = $(ML.Cᴷuʳ),  \n" *
              "     Cᴷcʳ = $(ML.Cᴷcʳ),  \n" *
              "     Cᴷeʳ = $(ML.Cᴷeʳ),  \n" *
              "    CᴷRiʷ = $(ML.CᴷRiʷ), \n" *
              "    CᴷRiᶜ = $(ML.CᴷRiᶜ)")


text:
@testset "multicolor_ordering" begin

end

ext: jl


text:
struct LTObsDist{M<:ObsModel}
    distances::ClearDistances
    model::M
end

obs_model(d::LTObsDist) = d.model

struct DESPOTEmu <: ObsModel
    std::Float64
    cdf::ReadingCDF
end

obs_type(::Type{DESPOTEmu}) = DMeas

gausscdf(mu, sigma, x) = (1+erf((x-mu)/(sigma*sqrt(2))))/2

function DESPOTEmu(f::Floor, std::Float64, maxread::Int=ceil(Int, max_diag(f)))
    maxclear = max(f.n_rows, f.n_cols) - 1
    cardcdf = Array{Float64}(undef, maxclear + 1, maxread + 1)

    for dist in 1:maxclear+1 # "dist = LaserRange" from cpp
        cdf = 0.0
        for reading in 0:maxread
            if reading < dist
                min_noise = reading - dist
                max_noise = min(dist, reading + 1) - dist
                cdf += 2*(gausscdf(0, std, max_noise)
                          - (reading > 0 ? gausscdf(0, std, min_noise) : 0.0))
                cardcdf[dist, reading+1] = cdf
            else
                cardcdf[dist, reading+1] = cdf
            end
        end
    end

    maxclear = min(f.n_rows, f.n_cols) - 1
    diagcdf = Array{Float64}(undef, maxclear + 1, maxread + 1)

    for c in 0:maxclear
        dist = sqrt(2)*(c+1)
        cdf = 0.0
        for reading in 0:maxread
            if reading < dist
                min_noise = reading - dist
                max_noise = min(dist, reading + 1) - dist
                cdf += 2*(gausscdf(0, std, max_noise)
                          - (reading > 0 ? gausscdf(0, std, min_noise) : 0.0))
                diagcdf[c+1, reading+1] = cdf
            else
                diagcdf[c+1, reading+1] = cdf
            end
        end
    end

    for c in 0:maxclear
        @assert abs(cardcdf[c+1, 1] + sum(diff(cardcdf[c+1, :])) - 1.0) < 1e-5 "cardinal cdf was $cardcdf"
        @assert abs(diagcdf[c+1, 1] + sum(diff(diagcdf[c+1, :])) - 1.0) < 1e-5 "diagonal cdf was $diagcdf"
    end

    return DESPOTEmu(std, ReadingCDF(cardcdf, diagcdf))
end

function Random.rand(rng::Random.AbstractRNG, d::LTObsDist{DESPOTEmu})
    if d.distances == CD_SAME_LOC
        return D_SAME_LOC
    end
    meas = MVector{8, Int}(undef)
    for i in 1:4
        meas[i] = max(0, floor(Int, n_clear_cells(d.distances, i)+1 - abs(d.model.std*randn(rng))))
    end
    for i in 5:8
        meas[i] = max(0, floor(Int, (n_clear_cells(d.distances, i)+1)*sqrt(2) - abs(d.model.std*randn(rng))))
    end
    return meas
end

function Distributions.pdf(d::LTObsDist{DESPOTEmu}, m::DMeas)
    a = d.distances == CD_SAME_LOC
    b = m == D_SAME_LOC
    if a || b
        return convert(Float64, a && b)
    end
    p = 1.0
    model = obs_model(d)
    for dir in 1:8
        if m[dir] == 0
            p *= cdf(model.cdf, dir, n_clear_cells(d.distances, dir), 0)
        else
            p *= (cdf(model.cdf, dir, n_clear_cells(d.distances, dir), m[dir]) -
                  cdf(model.cdf, dir, n_clear_cells(d.distances, dir), m[dir]-1))
        end
    end
    return p
end

function POMDPs.observation(p::LaserTagPOMDP{M}, sp::LTState) where M
    return LTObsDist(p.dcache[sp], p.obs_model)
end


text:
@testset "CoreTracker" begin
    @testset "General" begin
        F = equations(katsura(5))
        # test construction
        t1, start_sols = coretracker_startsolutions(F)

        test_show_juno(t1)
        @test_nowarn show(devnull, t1)
        @test !isempty(string(t1))
        test_show_juno(options(t1))
        @test_nowarn show(devnull, options(t1))
        test_show_juno(t1.state)
        @test_nowarn show(devnull, t1.state)
        @test !isempty(string(t1.state))

        @test_nowarn current_Δt(t1)

        caccuracy = accuracy(t1)
        @test_nowarn set_accuracy!(t1, 5e-5)
        @test accuracy(t1) == 5e-5
        set_accuracy!(t1, caccuracy)

        cmaxiter = max_corrector_iters(t1)
        @test_nowarn set_max_corrector_iters!(t1, 11)
        @test max_corrector_iters(t1) == 11
        set_max_corrector_iters!(t1, cmaxiter)

        @test t1 isa CoreTracker
        @test isa(current_x(t1), Vector)
        @test length(current_x(t1)) == 6

        t3, start_sols = coretracker_startsolutions(F, predictor = Euler())
        @test t3.predictor isa HC.EulerCache

        @test_deprecated setup!(t1, first(start_sols), 1.0, 0.4)
        init!(t1, first(start_sols), 1.0, 0.4)
        @test is_tracking(status(t1))
        @test current_t(t1) == 1.0

        setup!(t1, first(start_sols), 0.5, 0.4)
        @test is_invalid_startvalue(status(t1))
        @test is_terminated(status(t1))
        @test current_t(t1) == 0.5
        @test real(current_Δt(t1)) < 0
        @test cond(t1) ≥ 1

        R = track(t1, first(start_sols), 1.0, 0.0)
        @test R isa CoreTrackerResult
        @test is_success(R.returncode)
        @test R.accuracy < 1e-7
        @test_nowarn show(devnull, R)
        @test !isempty(string(R))
        test_show_juno(R)

        out = copy(first(start_sols))
        retcode = track!(out, t1, first(start_sols), 1.0, 0.0)
        @test is_success(retcode)
        @test out == solution(R)

        @test norm(t1) isa AbstractNorm
    end

    @testset "log_homotopy" begin
        F = equations(katsura(5))
        log_tracker, start_sols = coretracker_startsolutions(F; log_homotopy = true)
        @test log_tracker.homotopy.homotopy isa HC.LogHomotopy
        @test log_tracker.options.logarithmic_time_scale
        @test is_success(track!(log_tracker, first(start_sols), 0.0, 32.0))
    end

    @testset "deprecated" begin
        F = equations(katsura(5))
        # test construction
        t1 = coretracker(F)

        @test_warn(_ -> true, refinement_accuracy(t1))
        @test_warn(_ -> true, set_refinement_accuracy!(t1, 1e-4))
        @test_warn(_ -> true, max_refinement_iters(t1))
        @test_warn(_ -> true, set_max_refinement_iters!(t1, 11))
    end

    @testset "Affine tracking" begin
        @polyvar x a y b
        F = [x^2 - a, x * y - a + b]
        p = [a, b]

        tracker, starts = coretracker_startsolutions(
            F,
            [1.0, 1.0 + 0.0 * im],
            parameters = p,
            p₁ = [1, 0],
            p₀ = [2, 4],
            affine_tracking = true,
        )
        @test affine_tracking(tracker) == true
        res = track(tracker, starts[1], 1.0, 0.0)
        @test is_success(res)
        @test isa(res.x, Vector{ComplexF64})
        @test length(res.x) == 2

        x = @SVector [1.0, 1.0 + 0.0 * im]
        tracker, starts = coretracker_startsolutions(
            F,
            x,
            parameters = p,
            p₁ = [1, 0],
            p₀ = [2, 4],
            affine_tracking = true,
        )
        @test length(starts) == 1
        res = track(tracker, starts[1], 1.0, 0.0)
        @test isa(res.x, Vector{ComplexF64})
        @test length(res.x) == 2
    end

    # Broken on 1.4
    # @testset "Allocations" begin
    #     F = equations(cyclic(5))
    #     tracker, start_sols = coretracker_startsolutions(F; seed = 12356)
    #     s = first(start_sols)
    #     @allocated track!(tracker, s, 1.0, 0.01)
    #     # The path tracker should not allocate after the first tracked path
    #     @test (@allocated track!(tracker, s, 1.0, 0.01)) == 0
    #
    #     # log homotopy as well
    #     log_tracker, log_start_sols =
    #         coretracker_startsolutions(F; seed = 12356, log_homotopy = true)
    #     log_s = first(log_start_sols)
    #     @allocated track!(log_tracker, log_s, 0.0, -log(0.01))
    #     # The path tracker should not allocate after the first tracked path
    #     @test (@allocated track!(log_tracker, log_s, 0.0, -log(0.01))) == 0
    # end

    @testset "LinearSystem" begin
        A = rand(3, 3)
        b = rand(3)

        @polyvar x[1:3]
        F = A * x - b

        tracker, start_sols = coretracker_startsolutions(F; max_steps = 10)
        s = first(start_sols)
        result = track(tracker, s, 1.0, 0.0)
        @test is_success(result)
        @test result.t == 0.0
        x = result.x
        @test norm(x - A \ b) < 1e-6

        x_inter = copy(s)
        retcode = track!(x_inter, tracker, s, 1.0, 0.1)
        @test is_success(retcode)
        x_final = zero(x_inter)
        retcode = track!(x_final, tracker, x_inter, 0.1, 0.0)
        @test is_success(retcode)
        @test steps(tracker) < 4
        x = current_x(tracker)
        @test norm(x - A \ b) < 1e-6
    end

    @testset "FixedPatch continuation" begin
        F = homogenize(equations(katsura(5)))
        # test construction
        tracker, start_sols = coretracker_startsolutions(F, patch = OrthogonalPatch())
        patch = tracker.state.patch
        fixedtracker = coretracker(F, patch = FixedPatch())
        fixedpatch = fixedtracker.state.patch
        x1 = first(start_sols)

        r1 = track(fixedtracker, x1, 1.0, 0.1)
        v1 = copy(fixedpatch.v̄)
        track!(fixedtracker, r1.x, 0.1, 0.05, setup_patch = false)
        @test v1 == fixedpatch.v̄
        track!(fixedtracker, current_x(tracker), 0.05, 0.01, setup_patch = false)
        @test v1 == fixedpatch.v̄
        track!(fixedtracker, r1.x, 0.1, 0.05, setup_patch = true)
        @test v1 != fixedpatch.v̄

        r1 = track(tracker, x1, 1.0, 0.1)
        v1 = copy(patch.v̄)
        track!(tracker, r1.x, 0.1, 0.0)
        @test v1 != patch.v̄
    end

    @testset "coretracker_startsolutions" begin
        @polyvar x y
        tracker, S = coretracker_startsolutions([x^2 + y^2 - 4, x + y - 3])
        result = track(tracker, first(S), 1.0, 0.0)
        @test is_success(result)

        # test type promotion of startsolutions
        tracker, S = coretracker_startsolutions([x^2 + y^2 - 4, x + y - 3])
        result = track(tracker, [1, 1], 1.0, 0.0)
        @test is_success(result)

        tracker, S = coretracker_startsolutions(
            [x^2 + y^2 - 4, x + y - 3];
            simple_step_size_alg = true,
        )
        result = track(tracker, first(S), 1.0, 0.0)
        @test is_success(result)
    end

    @testset "iterator" begin
        @polyvar x a y b
        F = [x^2 - a, x * y - a + b]
        p = [a, b]
        s = [1.0, 1.0 + 0.0 * im]
        tracker = coretracker(F, parameters = p, p₁ = [1, 0], p₀ = [2, 4])
        init!(tracker, s, 1.0, 0.0)
        for tr in tracker
            # nothing
        end
        @test is_success(status(tracker))

        # path iterator
        typeof(first(iterator(tracker, s, 1.0, 0.0))) == Tuple{Vector{ComplexF64},Float64}

        @test max_step_size(tracker) == Inf
        set_max_step_size!(tracker, 0.01)
        @test max_step_size(tracker) == 0.01

        @test length(collect(iterator(tracker, s, 1.0, 0.0))) == 101

        @polyvar x a
        f = [x - a]
        ct = coretracker(
            f,
            [1.0],
            parameters = [a],
            p₁ = [1.0],
            p₀ = [2.0],
            max_step_size = 0.1,
        )
        Xs = Vector{ComplexF64}[]
        for (x, t) in iterator(ct, [1.0], 1.0, 0.0)
            push!(Xs, x)
        end
        @test round.(Int, real.(first.(Xs)) .* 10) == collect(10:20)
    end

    @testset "Limit Accuracy / Adaptive Precision" begin
        # We use Wilkinson polynomials of degree 12 and 15
        # 12 is fine for 64 bit precision, but in 15 is the evaluation error too large
        @polyvar x
        n = 12
        g = [x^n - 1]
        f = [prod(x - i for i = 1:n)]
        S = [[cis(i * 2π / n)] for i = 0:(n-1)]

        tracker = coretracker(
            g,
            f,
            S;
            log_homotopy = true,
            min_step_size = eps()^2,
            accuracy = 1e-7,
            precision = PRECISION_FIXED_64,
            seed = 131232,
        )
        results = map(s -> track(tracker, s, 0.0, 70), S)
        @test all(is_success, results)

        n = 15
        g = [x^n - 1]
        f = [prod(x - i for i = 1:n)]
        S = [[cis(i * 2π / n)] for i = 0:(n-1)]

        tracker = coretracker(
            g,
            f,
            S;
            log_homotopy = true,
            min_step_size = eps()^2,
            accuracy = 1e-7,
            precision = :double,
            seed = 131232,
        )
        results = map(s -> track(tracker, s, 0.0, 70), S)
        @test all(
            r -> is_success(
                r,
            ) || r.returncode == HC.CoreTrackerStatus.terminated_accuracy_limit,
            results,
        )

        tracker = coretracker(
            g,
            f,
            S;
            log_homotopy = true,
            min_step_size = eps()^2,
            accuracy = 1e-7,
            precision = :double_double,
            seed = 131232,
        )
        results = map(s -> track(tracker, s, 0.0, 70), S)
        @test all(is_success, results)

        tracker = coretracker(
            g,
            f,
            S;
            log_homotopy = true,
            min_step_size = eps()^2,
            accuracy = 1e-7,
            precision = :adaptive,
            seed = 131232,
        )
        results = map(s -> track(tracker, s, 0.0, 70), S)
        @test all(is_success, results)

        @test_throws ArgumentError coretracker(
            g,
            f,
            S;
            log_homotopy = true,
            precision = :NOT_A_PRECISION,
            seed = 131232,
        )
    end

    @testset "Change parameters" begin
        @polyvar x a y b
        F = SPSystem([x^2 - a, x * y - a + b]; parameters = [a, b])

        tracker, starts = coretracker_startsolutions(
            F,
            [1.0, 1.0 + 0.0 * im],
            generic_parameters = [2.2, 3.2],
        )
        start_parameters!(tracker, [1, 0])
        target_parameters!(tracker, [2, 4])
        @test is_success(track(tracker, starts[1], 1.0, 0.0))
    end

    @testset "(Multi-) projective" begin
        @polyvar x y u v
        f = [x * y - 6, x^2 - 5]
        tracker, starts = coretracker_startsolutions(
            f,
            variable_groups = [[x], [y]],
            seed = 123456,
            predictor = Pade21(),
        )
        current_x(tracker) isa ProjectiveVectors.PVector{ComplexF64,2}
        S = collect(starts)
        R1 = track(tracker, S[1], 1.0, 0.0)
        @test is_success(R1)
        @test solution(R1) isa ProjectiveVectors.PVector{ComplexF64,2}

        R2 = track(tracker, S[2], 1.0, 0.0)
        @test is_success(R2)
    end
end


text:
using Random
using Test
using FFTW
using CUDA
using TransformPlayground

Random.seed!(0)

FT = Float64
ε = eps(FT)

for N in (4, 5)
    A_cpu = randn(FT, N, N)
    A_gpu = CuArray(A_cpu)

    B_cpu = FFTW.r2r(A_cpu, FFTW.REDFT10)
    B_gpu = TransformPlayground.dct_makhoul_2d(A_gpu)

    @test isapprox.(B_cpu, Array(B_gpu), atol=16ε) |> all

    C_cpu = FFTW.r2r(B_cpu, FFTW.REDFT01)
    C_cpu = C_cpu ./ (2N)^2

    C_gpu = TransformPlayground.idct_makhoul_2d(B_gpu)

    @test isapprox.(C_cpu, Array(C_gpu), atol=ε) |> any

    @test isapprox.(A_cpu, C_cpu, atol=4ε) |> all
    @test isapprox.(Array(A_gpu), Array(C_gpu), atol=4ε) |> all
end


text:
function secret_handshake(code)
    handshake = []
    binarish = digits(code, base=2, pad=5)

    binarish[1] == 1 && (handshake = ["wink"])
    binarish[2] == 1 && (push!(handshake, "double blink"))
    binarish[3] == 1 && (push!(handshake, "close your eyes"))
    binarish[4] == 1 && (push!(handshake, "jump"))
    binarish[5] == 1 && (reverse!(handshake))
    return handshake

end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
"Call hid_vector() or hid_dict() to get references to human interface devices"
module WinControllerToFile
import REPL
using PyCall, REPL.TerminalMenus
export Hid, hid_dict, hid_vector
export documentation
export subscribe
import Base.show

struct Hid
    vendor_name::String
    vID::Int64
    pID::Int64
    device_path::String
    version_number::Int
    product_name::String
    object::PyObject
end
Hid() = Hid("", 0, 0, "", 0, "", PyObject(1) )
include("show_device.jl")

const SUBSCRIBETO = Array{String, 1}()
const ACTIVE_SUBSCRIPTIONS = Array{Hid, 1}()
const POLLINGTASKS = Array{Task, 1}()

"The minimum Δt is around 250 millisecond. The bottleneck seems to be PyWinUsb."
const ΔT = 0.25      # s
const TIMEOUT = 1800 # Half an hour

"Return a dictionary of usb human interface device references"
function hid_dict()
    pyhid = pyimport("pywinusb.hid")
    py_all_hids = pyhid.find_all_hid_devices()
    dic = Dict{String, Hid}()
    for h in py_all_hids
        pth = replace(h.device_path, r"[^a-zA-Z0-9_]" =>"")
        juhid = Hid(h.vendor_name,
            h.vendor_id,
            h.product_id,
            pth,
            h.version_number,
            h.product_name,
            h)
        nam = juhid.device_path
        push!(dic, nam => juhid)
    end
    dic
end

"Return a vector of usb human interface device references"
function hid_vector()
    collect(values(hid_dict()))
end

"Show full documentation for a device"
function documentation(io, h::Hid)
    pyhid = pyimport("pywinusb.hid")
    pyhid.core.show_hids(h.vID, h.pID, output = io)
    println(io, "\n\n")
    println(io, "Julia structure for device:")
    show(io, h)
    println(io, "\n\n")
    show(io, MIME("text/plain"), h)
end

open_hid(h::Hid) = h.object.open()
is_plugged(h::Hid) = h.object.is_plugged()
set_data_handler(i, h::Hid) = h.object.set_raw_data_handler(py"sample_handler$$i")

"The minimum Δt is around 250 millisecond. The bottleneck seems to be PyWinUsb."
function poller(timeout, Δt, h::Hid)
    t0 = time()
    while time()-t0 < timeout && is_plugged(h) # This triggers one sample written to file.
        # PyWinUsb tends to pile up threads. Not triggering too often may help avoiding that.
        sleep(Δt * 0.9)
    end
    @info "Exit polling of  $(h.vendor_name) $(h.product_name) after $timeout s "
end

"""
Define subscriptions and start reacting to events.
Subscriptions run asyncronously until this process is terminated,
although no events actually occur after timeout.
"""
function subscribe(;timeout = TIMEOUT, Δt = ΔT)
    if length(ACTIVE_SUBSCRIPTIONS) !=0
        error("Sorry, you need to restart this module since subscriptions are defined.")
    end
    fo = joinpath(homedir(), ".julia_hid")
    docfo = joinpath(fo, "DocumentationFo")
    hids = hid_dict()
    nohid = Hid()
    local count = 0
    for (i, su) in enumerate(SUBSCRIBETO)
        h = get(hids, su, nohid)
        fina = joinpath(fo, su * ".txt")
        if h == nohid
            @warn("The file $fina \n\t does not correspond to a connected Human Interface Device")
        else
            count += 1
            set_data_handler(i, h)
            open_hid(h)
            push!(ACTIVE_SUBSCRIPTIONS, h)
            task = @async poller(timeout, Δt, h)
            push!(POLLINGTASKS, task)
        end
        # Also write the full documentation for the subscribed device as a file
        docfina = joinpath(docfo, su * ".txt")
        @info "Writing documentation to $docfina"
        open(docfina, write=true) do io
            documentation(io, h)
        end
        println(stdout, h)
        @info "Updating $fina \n\t at $(Int(floor(1/Δt)))) Hz for $timeout s"
    end
    if count == 0
        @info "No subscription files X.txt found in $fo\n\t. "
        iob = IOBuffer()
        ioc = IOContext(iob, :color => true, :limit => false)
        show(ioc, "text/plain", values(hids))
        st=String(take!(iob))
        options = String.(split(st, "\n")[2:end])
        menu = MultiSelectMenu(options,  pagesize = length(options))
        choices = request("Select devices to subscribe to:", menu)
        vehids = collect(values(hids))
        for choice in choices
            device = vehids[choice]
            shfina = device.device_path * ".txt"
            fina = joinpath(fo, shfina)
            println(stdout, "\n", device, "\n\t\t\t\tStatus file ", shfina)
            open(fina, write = true) do io
                print(io, "\n")
            end
        end
        @info "Recompile WinControllerToFile to subscribe to updates based on connected devices and files in $fo"
    end
    POLLINGTASKS
end

function __init__()
    fo = joinpath(homedir(), ".julia_hid")
    if ispath(fo)
        filenames = filter(isfile, readdir(fo, join=true))
        if length(filenames) != 0
            docfo = joinpath(fo, "DocumentationFo")
            if ispath(docfo)
                local i = 0
                for fi in filenames
                    i += 1
                    subscription = replace(splitpath(fi)[end], ".txt" => "")
                    push!(SUBSCRIBETO, subscription)
                    # File name string enclosed by '"' and using forward slash
                    fina = "\"" * replace(fi, "\\" => "/") * "\""

                    # Define a Python handler for this object.
                    # The handler function name includes sequence 'i',
                    # so that we can handle multiple controllers.
                    py"""
                    import time
                    t0 = time.time()
                    t$$i = time.time()

                    def sample_handler$$i(data):
                        global t0
                        global t$$i
                        if (time.time() - t$$i) >= $$ΔT:
                            j = 1
                            headline = "#"
                            dataline = " "
                            for el in data:
                                headline +="chn{:<5d}".format(j)
                                dataline += "{:<8d}".format(el)
                                j +=1
                            timeline = "{} : time".format(time.time() - t0)
                            lines = headline + "\n" + dataline + "\n" + timeline
                            with open($$fina, 'w') as f:
                                f.write(lines)

                            t$$i = time.time()"""

                end
                @info("Generated $i Python handlers, run subscribe() when ready")
            else
                @warn("You need to create folder $docfo\n\tThen recompile WinControllerToFile")
            end
        else
            @warn "No subscription files X.txt found in $fo\n\t. Run subscribe() to pick devices and generate the corresponding files. "
        end
    else
        @warn("You need to create folder $fo\n\tThen recompile WinControllerToFile")
    end
end

end # module


text:
function _get_matrix_residue(msa::AbstractMultipleSequenceAlignment)::Matrix{Residue}
    getarray(namedmatrix(msa))
end
_get_matrix_residue(msa::NamedArray)::Matrix{Residue} = getarray(msa)
_get_matrix_residue(msa::Matrix{Residue})::Matrix{Residue} = msa

# Kernel function to fill ContingencyTable based on residues
function _mapfreq_kernel!(f::Function,
                          freqtable::Union{Probabilities{T,N,A},Counts{T,N,A}},
                          weights, pseudocounts, pseudofrequencies,
                          res::Vararg{AbstractVector{Residue},N}) where {T,N,A}
    table = freqtable.table
    _cleanup_table!(table) # count! calls _cleanup_temporal! and cleans marginals
    count!(table, weights, pseudocounts, res...) # count! calls apply_pseudocount! and  _update!
    if isa(freqtable, Probabilities{T,N,A})
        normalize!(table)
        apply_pseudofrequencies!(table, pseudofrequencies)
    end
    f(freqtable)
end

_mapfreq_kargs_doc = """
- `weights` (default: `NoClustering()`): Weights to be used for table counting.
- `pseudocounts` (default: `NoPseudocount()`): `Pseudocount` object to be applied to table.
- `pseudofrequencies` (default: `NoPseudofrequencies()`): `Pseudofrequencies` to be applied to the normalized (probabilities) table.
"""

_mappairfreq_kargs_doc = """
- `diagonalvalue` (default: `0`): Value to fill diagonal elements if `usediagonal` is `Val{false}`.
"""

# Residues: The output is a Named Vector
# --------------------------------------

function _mapfreq!(f::Function,
res_list::Vector{V}, # sequences or columns
table::Union{Probabilities{T,1,A},Counts{T,1,A}};
weights = NoClustering(),
pseudocounts::Pseudocount = NoPseudocount(),
pseudofrequencies::Pseudofrequencies = NoPseudofrequencies()) where {T,A,V<:AbstractArray{Residue}}
    scores = map(res_list) do res
        _mapfreq_kernel!(f, table, weights, pseudocounts, pseudofrequencies, res)
    end
    scores
end

# Map to each column

"""
It efficiently map a function (first argument) that takes a table of `Counts` or
`Probabilities` (third argument). The table is filled in place with the counts or
probabilities of each column from the `msa` (second argument).

$_mapfreq_kargs_doc
"""
function mapcolfreq!(f::Function, msa::AbstractMatrix{Residue},
                     table::Union{Probabilities{T,1,A},Counts{T,1,A}}; kargs...) where {T,A}
    N = ncolumns(msa)
    residues = _get_matrix_residue(msa)
    ncol = size(residues,2)
    columns = map(i -> view(residues,:,i), 1:ncol) # 2x faster than calling view inside the loop
    scores = _mapfreq!(f, columns, table; kargs...)
    name_list = columnnames(msa)
    NamedArray(reshape(scores, (1,N)),
               (OrderedDict{String,Int}(string(f) => 1),
                OrderedDict{String,Int}(name_list[i] => i for i in 1:N)),
               ("Function","Col"))
end

# Map to each sequence

"""
It efficiently map a function (first argument) that takes a table of `Counts` or
`Probabilities` (third argument). The table is filled in place with the counts or
probabilities of each sequence from the `msa` (second argument).

$_mapfreq_kargs_doc
"""
function mapseqfreq!(f::Function, msa::AbstractMatrix{Residue},
                     table::Union{Probabilities{T,1,A},Counts{T,1,A}}; kargs...) where {T,A}
    N = nsequences(msa)
    sequences = getresiduesequences(msa)
    name_list = sequencenames(msa)
    scores = _mapfreq!(f, sequences, table; kargs...)
    NamedArray(reshape(scores, (N,1)),
               (OrderedDict{String,Int}(name_list[i] => i for i in 1:N),
               OrderedDict{String,Int}(string(f) => 1)),
               ("Seq","Function"))
end

# Residue pairs: The output is a Named PairwiseListMatrix
# -------------------------------------------------------

function _mappairfreq!(f::Function,
res_list::Vector{V}, # sequences or columns
plm::PairwiseListMatrix{T,D,TV}, # output
table::Union{Probabilities{T,2,A},Counts{T,2,A}},
usediagonal::Type{Val{D}} = Val{true};
weights = NoClustering(),
pseudocounts::Pseudocount = NoPseudocount(),
pseudofrequencies::Pseudofrequencies = NoPseudofrequencies(),
diagonalvalue::T = zero(T)) where {T,D,TV,A,V<:AbstractArray{Residue}} # diagonalvalue used by mapcolpairfreq! ...
    @inbounds @iterateupper plm D begin
        list[k] = :($_mapfreq_kernel!)(:($f), :($table), :($weights),
                                       :($pseudocounts), :($pseudofrequencies),
                                       :($res_list)[i], :($res_list)[j])
    end
    plm
end

# Map to column pairs

"""
It efficiently map a function (first argument) that takes a table of `Counts` or
`Probabilities` (third argument). The table is filled in place with the counts or
probabilities of each pair of columns from the `msa` (second argument).
The fourth positional argument `usediagonal` indicates if the function should be applied
to identical element pairs (default to `Val{true}`).

$_mapfreq_kargs_doc
$_mappairfreq_kargs_doc
"""
function mapcolpairfreq!(f::Function, msa::AbstractMatrix{Residue},
                         table::Union{Probabilities{T,2,A},Counts{T,2,A}},
                         usediagonal::Type{Val{D}} = Val{true};
                         diagonalvalue::T = zero(T),
                         kargs...) where {T,A,D}
    ncol = ncolumns(msa)
    residues = _get_matrix_residue(msa)
    columns = map(i -> view(residues,:,i), 1:ncol) # 2x faster than calling view inside the loop
    scores = columnpairsmatrix(msa, T, Val{D}, diagonalvalue) # Named PairwiseListMatrix
    plm = getarray(scores)::PairwiseListMatrix{T,D,Vector{T}} # PairwiseListMatrix
    _mappairfreq!(f, columns, plm, table, Val{D}; kargs...)
    scores
end

# Map to sequence pairs

"""
It efficiently map a function (first argument) that takes a table of `Counts` or
`Probabilities` (third argument). The table is filled in place with the counts or
probabilities of each pair of sequences from the `msa` (second argument).
The fourth positional argument `usediagonal` indicates if the function should be applied
to identical element pairs (default to `Val{true}`).

$_mapfreq_kargs_doc
$_mappairfreq_kargs_doc
"""
function mapseqpairfreq!(f::Function, msa::AbstractMatrix{Residue},
                         table::Union{Probabilities{T,2,A},Counts{T,2,A}},
                         usediagonal::Type{Val{D}} = Val{true};
                         diagonalvalue::T = zero(T),
                         kargs...) where {T,A,D}
    sequences = getresiduesequences(msa)
    scores = sequencepairsmatrix(msa, T, Val{D}, diagonalvalue) # Named PairwiseListMatrix
    plm = getarray(scores)::PairwiseListMatrix{T,D,Vector{T}} # PairwiseListMatrix
    _mappairfreq!(f, sequences, plm, table, Val{D}; kargs...)
    scores
end

# cMI
# ===

"""
`cumulative` allows to calculate cumulative scores (i.e. cMI) as defined in Buslje et. al. 2010

*"We calculated a cumulative mutual information score (cMI) for each residue as the sum of
MI values above a certain threshold for every amino acid pair where the particular residue
appears. This value defines to what degree a given amino acid takes part in a mutual
information network."*
Buslje, Cristina Marino, Elin Teppa, Tomas Di Doménico, José María Delfino, and Morten
Nielsen. *Networks of high mutual information define the structural proximity of catalytic
sites: implications for catalytic residue identification.* PLoS Comput Biol 6, no. 11
(2010): e1000978.
"""
function cumulative(plm::PairwiseListMatrix{T,D,VT}, threshold::T) where {T,D,VT}
    N = size(plm, 1)
    out = zeros(T, N)
    @iterateupper plm false begin
        elem = list[k]
        if !isnan(elem) && elem >= :($threshold)
            :($out)[i] += elem
            :($out)[j] += elem
        end
    end
    reshape(out, (1,N))
end

function cumulative(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN},
                    threshold::T) where {T,D,TV,DN}
    N = size(nplm, 2)
    name_list = names(nplm, 2)
    NamedArray(cumulative(getarray(nplm), threshold),
               (OrderedDict{String,Int}("cumulative" => 1),
               OrderedDict{String,Int}(name_list[i] => i for i in 1:N)),
               ("Function", dimnames(nplm,2)))
end


text:
using ImmutableArrays

function ratio(f::Function, a1, a2, b1, b2, iters::Integer)
    if f(a1, a2) != f(b1, b2)
        return string("result mismatch:\n", f(a1, a2), "\n--- vs. ---\n", f(b1, b2))
    end
    tic()
    for i = 1:iters
        f(a1, a2)
    end
    atime = toq()
    tic()
    for i = 1:iters
        f(b1, b2)
    end
    btime = toq()
    return btime / atime
end
ratio(f::Function, a1, a2, b1, b2) = ratio(f, a1, a2, b1, b2, 1000000)
mat4ratio(f::Function, m1, m2) = ratio(f, m1, m2, Matrix4x4(m1), Matrix4x4(m2))

A = rand(4,4)
B = rand(4,4)

println("ImmutableArray operation speed relative to native Julia arrays:")
println("+ ", mat4ratio(+,A,B))
println("* ", mat4ratio(*,A,B))


text:
const HEADER_SIZE = 227
const DATA_OFFSET = 227
const SIZE_DATARECORD = 26

# """
# save file .las
# """
# function save_new_las(filename::String,h::LasIO.LasHeader,p::LasIO.Array{LasPoint,1})
#     if ispath(filename) #overwrite
#         rm(filename)
#     end
#     LasIO.FileIO.save(filename,h,p)
# end

"""
create a header
"""
function newHeader(aabb::AABB, software::String, sizePointRecord, npoints=0; scale=0.001)

	file_source_id=UInt16(0)
	global_encoding=UInt16(0)
	guid_1=UInt32(0)
	guid_2=UInt16(0)
	guid_3=UInt16(0)
	guid_4=""
	version_major=UInt8(1)
	version_minor=UInt8(2)
	system_id=""
	software_id = software
	creation_dayofyear = UInt16(Dates.dayofyear(today()))
	creation_year = UInt16(Dates.year(today()))
	header_size=UInt16(HEADER_SIZE) # valore fisso
	data_offset=UInt16(DATA_OFFSET) #valore fisso
	n_vlr=UInt32(0)
	data_format_id=UInt8(2)
	data_record_length=UInt16(sizePointRecord) #valore variabile
	records_count=UInt32(npoints)
	point_return_count=UInt32[0,0,0,0,0]
	x_scale=scale
	y_scale=scale
	z_scale=scale
	x_offset = aabb.x_min
	y_offset = aabb.y_min
	z_offset = aabb.z_min
	x_max = aabb.x_max
	x_min = aabb.x_min
	y_max = aabb.y_max
	y_min = aabb.y_min
	z_max = aabb.z_max
	z_min = aabb.z_min
	variable_length_records=Vector{LasVariableLengthRecord}()
	user_defined_bytes=Vector{UInt8}()

	return LasIO.LasHeader(file_source_id,
    global_encoding,
    guid_1,
    guid_2,
    guid_3,
    guid_4,
    version_major,
    version_minor,
    system_id,
    software_id,
    creation_dayofyear,
    creation_year,
    header_size,
    data_offset,
    n_vlr,
    data_format_id,
    data_record_length,
    records_count,
    point_return_count,
    x_scale,
    y_scale,
    z_scale,
    x_offset,
    y_offset,
    z_offset,
    x_max,
    x_min,
    y_max,
    y_min,
    z_max,
    z_min,
    variable_length_records,
    user_defined_bytes
	)
end

"""
create a new point record for las file.
"""
function newPointRecord(laspoint::LasIO.LasPoint, header::LasIO.LasHeader, type::DataType, mainHeader::LasIO.LasHeader; affineMatrix = Matrix{Float64}(Common.I,4,4))

	point = Common.apply_matrix(affineMatrix,[xcoord(laspoint,header), ycoord(laspoint,header), zcoord(laspoint,header)])
	x = LasIO.xcoord(point[1],mainHeader)
	y = LasIO.ycoord(point[2],mainHeader)
	z = LasIO.zcoord(point[3],mainHeader)
	intensity = laspoint.intensity
	flag_byte = laspoint.flag_byte
	raw_classification = laspoint.raw_classification
	scan_angle = laspoint.scan_angle
	user_data = laspoint.user_data
	pt_src_id = laspoint.pt_src_id

	if type == LasIO.LasPoint0
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id
					)

	elseif type == LasIO.LasPoint1
		gps_time = laspoint.gps_time
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id, gps_time
					)

	elseif type == LasIO.LasPoint2
		red = laspoint.red
		green = laspoint.green
		blue = laspoint.blue
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id,
					red, green, blue
					)

	elseif type == LasIO.LasPoint3
		gps_time = laspoint.gps_time
		red = laspoint.red
		green = laspoint.green
		blue = laspoint.blue
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id, gps_time,
					red, green, blue
					)

	end

end

function newPointRecord(point::Array{Float64,1}, rgb::Array{LasIO.N0f16,1} , type::LasIO.DataType, mainHeader::LasIO.LasHeader) #crea oggetto laspoint con vertici e colori

	x = LasIO.xcoord(point[1],mainHeader)
	y = LasIO.ycoord(point[2],mainHeader)
	z = LasIO.zcoord(point[3],mainHeader)
	intensity = UInt16(0)
	flag_byte = UInt8(0)
	raw_classification = UInt8(0)
	scan_angle = Int8(0)
	user_data = UInt8(0)
	pt_src_id = UInt16(0)

	if type == LasIO.LasPoint0
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id
					)

	elseif type == LasIO.LasPoint1
		gps_time = Float64(0)
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id, gps_time
					)

	elseif type == LasIO.LasPoint2
		red = rgb[1]
		green = rgb[2]
		blue = rgb[3]
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id,
					red, green, blue
					)

	elseif type == LasIO.LasPoint3
		gps_time = Float64(0)
		red = rgb[1]
		green = rgb[2]
		blue = rgb[3]
		return type(x, y, z,
					intensity, flag_byte, raw_classification,
					scan_angle, user_data, pt_src_id, gps_time,
					red, green, blue
					)

	end

end

"""
Save a point cloud in las format
"""
function save_pointcloud(filename::String, pc::PointCloud,  software::String)

	if pc.dimension == 2
		points = vcat(pc.coordinates,zeros(pc.n_points)')
	else
		points = pc.coordinates
	end

	rgbs = pc.rgbs
	npoints = pc.n_points

	aabb = AABB(points)
	header = newHeader(aabb,software,SIZE_DATARECORD,npoints)

	pvec = Array{LasPoint,1}(undef,npoints)
	for i in 1:npoints
		point = newPointRecord(points[:,i], rgbs[:,i], LasIO.LasPoint2, header)
		pvec[i] = point
	end

	LasIO.FileIO.save(filename,header,pvec)
end

#
# """
#  	createlasdata(p,h,header)
#
# create LasPoint from coordinates.
# """
# function createlasdata(p,h::LasIO.LasHeader,mainHeader::LasIO.LasHeader)
# 	type = pointformat(h)
#
# 	x = LasIO.xcoord(xcoord(p,h),mainHeader)
# 	y = LasIO.ycoord(ycoord(p,h),mainHeader)
# 	z = LasIO.zcoord(zcoord(p,h),mainHeader)
# 	intensity = p.intensity
# 	flag_byte = p.flag_byte
# 	raw_classification = p.raw_classification
# 	scan_angle = p.scan_angle
# 	user_data = p.user_data
# 	pt_src_id = p.pt_src_id
#
# 	if type == LasIO.LasPoint0
# 		return type(x, y, z,
# 					intensity, flag_byte, raw_classification,
# 					scan_angle, user_data, pt_src_id
# 					)
#
# 	elseif type == LasIO.LasPoint1
# 		gps_time = p.gps_time
# 		return type(x, y, z,
# 					intensity, flag_byte, raw_classification,
# 					scan_angle, user_data, pt_src_id, gps_time
# 					)
#
# 	elseif type == LasIO.LasPoint2
# 		red = p.red
# 		green = p.green
# 		blue = p.blue
# 		return type(x, y, z,
# 					intensity, flag_byte, raw_classification,
# 					scan_angle, user_data, pt_src_id,
# 					red, green, blue
# 					)
#
# 	elseif type == LasIO.LasPoint3
# 		gps_time = p.gps_time
# 		red = p.red
# 		green = p.green
# 		blue = p.blue
# 		return type(x, y, z,
# 					intensity, flag_byte, raw_classification,
# 					scan_angle, user_data, pt_src_id, gps_time,
# 					red, green, blue
# 					)
#
# 	end
# end

# """
# 	bbincremental!(coordpoint,bb)
#
# """
# function bbincremental!(coordpoint,bb)
#
# 	for i in 1:length(coordpoint)
# 		if coordpoint[i] < bb[1][i]
# 			bb[1][i] = coordpoint[i]
# 		end
# 		if coordpoint[i] > bb[2][i]
# 			bb[2][i] = coordpoint[i]
# 		end
# 	end
#
# 	return true
# end

# """
# .
# """
#
# function set_z_zero(points::Array{LasPoint2,1},header::LasIO.LasHeader)
# 	type = pointformat(header)
# 	pvec = Vector{LasPoint}()
# 	for p in points
# 		x = p.x
# 		y = p.y
# 		z = 0
# 		intensity = p.intensity
# 		flag_byte = p.flag_byte
# 		raw_classification = p.raw_classification
# 		scan_angle = p.scan_angle
# 		user_data = p.user_data
# 		pt_src_id = p.pt_src_id
#
# 		if type == LasIO.LasPoint0
# 			laspoint = type(x, y, z,
# 						intensity, flag_byte, raw_classification,
# 						scan_angle, user_data, pt_src_id
# 						)
#
# 		elseif type == LasIO.LasPoint1
# 			gps_time = p.gps_time
# 			laspoint = type(x, y, z,
# 						intensity, flag_byte, raw_classification,
# 						scan_angle, user_data, pt_src_id, gps_time
# 						)
#
# 		elseif type == LasIO.LasPoint2
# 			red = p.red
# 			green = p.green
# 			blue = p.blue
# 			laspoint = type(x, y, z,
# 						intensity, flag_byte, raw_classification,
# 						scan_angle, user_data, pt_src_id,
# 						red, green, blue
# 						)
#
# 		elseif type == LasIO.LasPoint3
# 			gps_time = p.gps_time
# 			red = p.red
# 			green = p.green
# 			blue = p.blue
# 			laspoint = type(x, y, z,
# 						intensity, flag_byte, raw_classification,
# 						scan_angle, user_data, pt_src_id, gps_time,
# 						red, green, blue
# 						)
#
# 		end
# 		push!(pvec,laspoint)
# 	end
# 	return pvec
# end

# """
# 	mergelas(headers,pointdata,bb,scale)
#
# Merge more file .las.
# """
# function mergelas(headers,pointdata)
# 	@assert length(headers) == length(pointdata) "mergelas: inconsistent data"
#
# 	# header of merging las
# 	hmerge = createheader(headers,pointdata)
# 	data = LasIO.LasPoint[]
#
# 	# Las point data merge
# 	for i in 1:length(pointdata)
# 		for p in pointdata[i]
# 			laspoint = createlasdata(p,headers[i],hmerge)
# 			push!(data,laspoint)
# 		end
# 	end
#
# 	return hmerge,data
# end

# """
#  	createheader(headers,pointdata,bb,scale)
#
# crea header coerente con i miei punti.
# """
# function createheader(headers,pointdata)
# 	type = pointformat(headers[1])
# 	h = deepcopy(headers[1])
# 	h.records_count = sum(length.(pointdata))
# 	return h
# end


text:
module TP

# O código do pacote

export fteste1
export fteste2

function fteste1(x)
    return "valor $x"
end

function fteste2()
	return "teste2 OK"
end

function fteste3()
	return "teste3 OK"
end

end

max_stars_repo_name: jmarcello-teste/TP

max_issues_repo_name: jmarcello-teste/TP

max_forks_repo_name: jmarcello-teste/TP


text:
module Models

export train!,
  predict,
  predict_response

# A model (m) must implement the following functions:
#
# train!(m, X[, y]) - train model with pattern X and a label y with it's associative
# predict(m, X) - returns the response of the model for pattern X

abstract type AbstractModel end

function train! end
function predict end
function predict_response end

# I can provide a default implementation of train! and predict for matrix inputs,
# assuming that their vector counterparts have already been implemented somewhere
# else

include("Nodes.jl")

include("Discriminator.jl")
export Discriminator,
  BitDiscriminator,
  BleachingDiscriminator,
  RegressionDiscriminator,
  GeneralizedRegressionDiscriminator

include("MultiDiscriminatorClassifier.jl")
export MultiDiscriminatorClassifier

end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
"""
   normalizationconstant

gives the normalization constants so that the jacobi matrix is symmetric,
that is, so we have orthonormal OPs:

    Q == P*normalizationconstant(P)
"""
function normalizationconstant(μ, P::AbstractQuasiMatrix{T}) where T
    X = jacobimatrix(P)
    c,b = subdiagonaldata(X),supdiagonaldata(X)
    # hide array type to avoid crazy compilation
    Accumulate{T,1,typeof(*),Vector{T},AbstractVector{T}}(*, T[μ], Vcat(zero(T),sqrt.(c ./ b)), 1, (1,))
end

normalizationconstant(P::AbstractQuasiMatrix) = normalizationconstant(inv(sqrt(sum(orthogonalityweight(P)))), P)

struct Normalized{T, OPs<:AbstractQuasiMatrix{T}, NL} <: OrthogonalPolynomial{T}
    P::OPs
    scaling::NL # Q = P * Diagonal(scaling)
end

Normalized(P::AbstractQuasiMatrix{T}) where T = Normalized(P, normalizationconstant(P))
Normalized(Q::Normalized) = Q
normalized(P) = Normalized(P)

struct NormalizedOPLayout{LAY<:AbstractBasisLayout} <: AbstractOPLayout end

MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedOPLayout{typeof(MemoryLayout(OPs))}()

struct QuasiQR{T, QQ, RR} <: Factorization{T}
    Q::QQ
    R::RR
end

QuasiQR(Q::AbstractQuasiMatrix{T}, R::AbstractMatrix{V}) where {T,V} =
    QuasiQR{promote_type(T,V),typeof(Q),typeof(R)}(Q, R)

Base.iterate(S::QuasiQR) = (S.Q, Val(:R))
Base.iterate(S::QuasiQR, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::QuasiQR, ::Val{:done}) = nothing

axes(Q::Normalized) = axes(Q.P)
==(A::Normalized, B::Normalized) = A.P == B.P

# There is no point in a Normalized OP thats ==, so just return false
==(A::Normalized, B::OrthogonalPolynomial) = false
==(A::OrthogonalPolynomial, B::Normalized) = false
==(A::Normalized, B::AbstractQuasiMatrix) = false
==(A::AbstractQuasiMatrix, B::Normalized) = false
==(A::Normalized, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = false
==(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::Normalized) = false

_p0(Q::Normalized) = Q.scaling[1]

# x * p[n] = c[n-1] * p[n-1] + a[n] * p[n] + b[n] * p[n+1]
# x * q[n]/h[n] = c[n-1] * q[n-1]/h[n-1] + a[n] * q[n]/h[n] + b[n] * q[n+1]/h[n+1]
# x * q[n+1] = c[n-1] * h[n]/h[n-1] * q[n-1] + a[n] * q[n] + b[n] * h[n]/h[n+1] * q[n+1]

# q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
# q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}

function symtridiagonalize(X)
    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
    SymTridiagonal(a, sqrt.(b .* c))
end
jacobimatrix(Q::Normalized) = symtridiagonalize(jacobimatrix(Q.P))

orthogonalityweight(Q::Normalized) = orthogonalityweight(Q.P)
singularities(Q::Normalized) = singularities(Q.P)

function demap(Q::Normalized)
    P,D =  arguments(ApplyLayout{typeof(*)}(), Q)
    demap(P) * D
end

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

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

ArrayLayouts.mul(Q::Normalized, C::AbstractArray) = ApplyQuasiArray(*, Q, C)

grid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = grid(view(parent(Q).P, parentindices(Q)...))
plotgrid(Q::SubQuasiArray{<:Any,2,<:Normalized,<:Tuple{Inclusion,Any}}) = plotgrid(view(parent(Q).P, parentindices(Q)...))

# transform_ldiv(Q::Normalized, C::AbstractQuasiArray) = Q.scaling .\ (Q.P \ C)
function transform_ldiv(Q::Normalized, C::AbstractQuasiArray)
    c = paddeddata(Q.P \ C)
    [Q.scaling[axes(c,1)] .\ c; zeros(eltype(c), ∞)]
end

function transform_ldiv(V::SubQuasiArray{<:Any,2,<:Normalized}, C::AbstractQuasiArray)
    Q = parent(V)
    P = Q.P
    kr, jr = parentindices(V)
    c = transform_ldiv(view(P, kr, jr), C)
    Q.scaling[axes(c,1)] .\ c
end

arguments(::ApplyLayout{typeof(*)}, Q::Normalized) = Q.P, Diagonal(Q.scaling)
_mul_arguments(Q::Normalized) = arguments(ApplyLayout{typeof(*)}(), Q)
_mul_arguments(Q::QuasiAdjoint{<:Any,<:Normalized}) = arguments(ApplyLayout{typeof(*)}(), Q)

# table stable identity if A.P == B.P
@inline _normalized_ldiv(An, C, Bn) = An \ (C * Bn)
@inline _normalized_ldiv(An, C::Eye{T}, Bn) where T = FillArrays.SquareEye{promote_type(eltype(An),T,eltype(Bn))}(ℵ₀)
@inline copy(L::Ldiv{<:NormalizedOPLayout,<:NormalizedOPLayout}) = _normalized_ldiv(Diagonal(L.A.scaling), L.A.P \ L.B.P, Diagonal(L.B.scaling))
@inline copy(L::Ldiv{Lay,<:NormalizedOPLayout}) where Lay = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,Lay}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{Lay,<:NormalizedOPLayout}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{Lay,<:NormalizedOPLayout}) where Lay<:AbstractLazyLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{ApplyLayout{typeof(*)},<:NormalizedOPLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))
for Lay in (:(ApplyLayout{typeof(*)}),:(BroadcastLayout{typeof(+)}),:(BroadcastLayout{typeof(-)}))
    @eval begin
        @inline copy(L::Ldiv{<:NormalizedOPLayout,$Lay}) = copy(Ldiv{ApplyLayout{typeof(*)},$Lay}(L.A, L.B))
        @inline copy(L::Ldiv{<:NormalizedOPLayout,$Lay,<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},$Lay}(L.A, L.B))
    end
end

copy(L::Ldiv{Lay,<:NormalizedOPLayout}) where Lay<:MappedBasisLayouts = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))

# want to use special re-expansion routines without expanding Normalized basis
@inline copy(L::Ldiv{<:NormalizedOPLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedOPLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))

# take out diagonal scaling for Weighted(::Normalized)
function _norm_expand_ldiv(A, w_B)
    w,B = w_B.args
    B̃,D = arguments(ApplyLayout{typeof(*)}(), B)
    (A \ (w .* B̃)) * D
end
copy(L::Ldiv{<:NormalizedOPLayout,<:WeightedBasisLayout{<:NormalizedOPLayout}}) = _norm_expand_ldiv(L.A, L.B)
copy(L::Ldiv{OPLayout,<:WeightedBasisLayout{<:NormalizedOPLayout}}) = _norm_expand_ldiv(L.A, L.B)

###
# show
###
show(io::IO, Q::Normalized) = print(io, "Normalized($(Q.P))")
show(io::IO, ::MIME"text/plain", Q::Normalized) = show(io, Q)

abstract type AbstractWeighted{T} <: Basis{T} end

getindex(Q::AbstractWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = weight(Q)[x] .* unweighted(Q)[x,jr]

MemoryLayout(::Type{<:AbstractWeighted}) = WeightedBasisLayout{OPLayout}()
convert(::Type{WeightedBasis}, Q::AbstractWeighted) = weight(Q) .* unweighted(Q)

# make act like WeightedBasisLayout
ContinuumArrays.sublayout(::WeightedOPLayout{Lay}, inds::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) where Lay = sublayout(WeightedBasisLayout{Lay}(), inds)
ContinuumArrays.sublayout(::WeightedOPLayout{Lay}, inds::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) where Lay = sublayout(WeightedBasisLayout{Lay}(), inds)

ContinuumArrays.unweighted(wP::AbstractWeighted) = wP.P
# function copy(L::Ldiv{WeightedOPLayout,WeightedOPLayout})
#     L.A.P == L.B.P && return Eye{eltype(L)}(∞)
#     convert(WeightedOrthogonalPolynomial, L.A) \ convert(WeightedOrthogonalPolynomial, L.B)
# end

# copy(L::Ldiv{WeightedOPLayout,<:BroadcastLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
# copy(L::Ldiv{WeightedOPLayout,<:AbstractLazyLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
# copy(L::Ldiv{WeightedOPLayout,<:AbstractBasisLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
# copy(L::Ldiv{<:AbstractLazyLayout,WeightedOPLayout}) = L.A \ convert(WeightedOrthogonalPolynomial, L.B)
# copy(L::Ldiv{<:AbstractBasisLayout,WeightedOPLayout}) = L.A \ convert(WeightedOrthogonalPolynomial, L.B)
# copy(L::Ldiv{WeightedOPLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
# copy(L::Ldiv{WeightedOPLayout,<:ExpansionLayout}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
# copy(L::Ldiv{WeightedOPLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))

copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
copy(L::Ldiv{Lay,<:WeightedOPLayout{<:NormalizedOPLayout}}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},<:WeightedOPLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
copy(L::Ldiv{Lay,<:WeightedBasisLayout{<:NormalizedOPLayout}}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},<:WeightedBasisLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},<:WeightedBasisLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},<:WeightedOPLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))

# function layout_broadcasted(::ExpansionLayout{WeightedOPLayout}, ::OPLayout, ::typeof(*), a, P)
#     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
#     wQ,c = arguments(a)
#     w,Q = arguments(wQ)
#     (w .* P) * Clenshaw(Q * c, P)
# end

"""
    OrthonormalWeighted(P)

is the orthonormal with respect to L^2 basis given by
`sqrt.(orthogonalityweight(P)) .* Normalized(P)`.
"""
struct OrthonormalWeighted{T, PP<:AbstractQuasiMatrix{T}} <: AbstractWeighted{T}
    P::Normalized{T, PP}
end

function OrthonormalWeighted(P)
    Q = normalized(P)
    OrthonormalWeighted{eltype(Q),typeof(P)}(Q)
end

axes(Q::OrthonormalWeighted) = axes(Q.P)
copy(Q::OrthonormalWeighted) = Q

==(A::OrthonormalWeighted, B::OrthonormalWeighted) = A.P == B.P

weight(Q::OrthonormalWeighted) = sqrt.(orthogonalityweight(Q.P))

broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))

"""
    Weighted(P)

is equivalent to `orthogonalityweight(P) .* P`
"""
struct Weighted{T, PP<:AbstractQuasiMatrix{T}} <: AbstractWeighted{T}
    P::PP
end

axes(Q::Weighted) = axes(Q.P)
copy(Q::Weighted) = Q

weight(wP::Weighted) = orthogonalityweight(wP.P)

_weightedmemorylayout(::PP) where PP<:AbstractOPLayout =  WeightedOPLayout{PP}()
_weightedmemorylayout(::PP) where PP =  WeightedBasisLayout{PP}()
MemoryLayout(::Type{<:Weighted{<:Any,PP}}) where PP = _weightedmemorylayout(MemoryLayout(PP))

function arguments(::ApplyLayout{typeof(*)}, Q::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{Weight,Normalized}})
    w,Q = Q.args
    P,D = arguments(*,Q)
    (w .* P),D
end
function arguments(::ApplyLayout{typeof(*)}, Q::Weighted{<:Any,<:Normalized})
    P,D = arguments(*, Q.P)
    Weighted(P),D
end
_mul_arguments(Q::Weighted{<:Any,<:Normalized}) = arguments(ApplyLayout{typeof(*)}(), Q)

# convert(::Type{WeightedOrthogonalPolynomial}, P::Weighted) = weight(P) .* unweighted(P)

broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::Weighted) = Q * (Q.P \ (x .* Q.P))

@simplify *(Ac::QuasiAdjoint{<:Any,<:Weighted}, wB::Weighted) =
    convert(WeightedBasis, parent(Ac))' * convert(WeightedBasis, wB)

@simplify function *(Ac::QuasiAdjoint{<:Any,<:Weighted}, B::AbstractQuasiVector)
    P = (Ac').P
    massmatrix(P) * (P\B)
end

@simplify function *(Ac::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:SubQuasiArray}}, B::Weighted{<:Any,<:SubQuasiArray})
    P = (Ac').P
    Q = B.P
    V = view(Weighted(parent(P)), parentindices(P)...)
    W = view(Weighted(parent(Q)), parentindices(Q)...)
    V'W
end

# Derivative overloads, see ContinuumArrays.jl/src/bases.jl

simplifiable(::typeof(*), A::Derivative, B::Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}) = simplifiable(*, Derivative(axes(parent(B),1)), Weighted(parent(B.P)))
simplifiable(::typeof(*), Ac::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = simplifiable(*, Bc', Ac')
function mul(A::Derivative, B::Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
    P = Weighted(parent(B.P))
    kr,jr = parentindices(B.P)
    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
end
mul(Ac::QuasiAdjoint{<:Any, Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}}}, Bc::QuasiAdjoint{<:Any,<:Derivative}) = mul(Bc', Ac')'

summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")

__sum(::NormalizedOPLayout, A, dims) = __sum(ApplyLayout{typeof(*)}(), A, dims)
function __sum(::WeightedOPLayout, A, dims)
    @assert dims == 1
    Hcat(sum(weight(A)), Zeros{eltype(A)}(1,∞))
end

_sum(p::SubQuasiArray{T,1,<:Weighted,<:Tuple{Inclusion,Int}}, ::Colon) where T =
    parentindices(p)[2] == 1 ? convert(T, sum(weight(parent(p)))) : zero(T)


text:
using MPI, Test
include(joinpath("..","..","testhelpers.jl"))

@testset "DGmethods" begin
    tests = [
        (1, "vars_test.jl")
        (1, "integral_test.jl")
        (1, "integral_test_sphere.jl")
        (1, "Euler/isentropicvortex.jl")
        (1, "Euler/isentropicvortex_imex.jl")
        (1, "Euler/isentropicvortex_multirate.jl")
        (1, "advection_diffusion/pseudo1D_advection_diffusion.jl")
        (1, "compressible_Navier_Stokes/ref_state.jl")
        (1, "horizontal_integral_test.jl")
        (2, "courant.jl")
    ]

    runmpi(tests, @__FILE__)
end


text:
using EconometricsLH
using Random, Test

function regression_test_test()
    @testset "Regression test" begin
        rng = MersenneTwister(44);
        n = Int(1e4);
        trueM = randn(rng, n);
        simM  = trueM .+ 0.3 .* randn(rng, n);
        wtM = 1.0 .+ rand(rng, n);

        for useWeights in [false, true]
            if useWeights
                dev, betaV, seV = regression_test(trueM, simM; weights = wtM);
            else
                dev, betaV, seV = regression_test(trueM, simM);
            end

            @test size(betaV) == (2,)
            @test size(seV) == (2,)
            @test all(seV .> 0.0)
            @test abs(betaV[1]) < 2.2 * seV[1]
            @test abs(betaV[2] - 1.0) < 2.2 * seV[2]
            @test (dev < 2.2);
        end
    end
end

@testset "Regression test" begin
    regression_test_test();
end

# ---------------


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestReverseAD

using Test
import LinearAlgebra
import MathOptInterface
import SparseArrays

const MOI = MathOptInterface
const Nonlinear = MOI.Nonlinear
const ReverseAD = Nonlinear.ReverseAD
const Coloring = ReverseAD.Coloring

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_objective_quadratic_univariate()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :($x^2 + 1))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2]) == 1.2^2 + 1
    g = [NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2])
    @test g == [2.4]
    @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1)]
    H = [NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, [1.2], 1.5, Float64[])
    @test H == 1.5 .* [2.0]
    MOI.eval_hessian_lagrangian_product(
        evaluator,
        H,
        [1.2],
        [1.2],
        1.5,
        Float64[],
    )
    @test H == [1.5 * 2.0 * 1.2]
    return
end

function test_objective_and_constraints_quadratic_univariate()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :($x^2 + 1))
    Nonlinear.add_constraint(model, :($x^2), MOI.LessThan(2.0))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2]) == 1.2^2 + 1
    g = [NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2])
    @test g == [2.4]
    @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1), (1, 1)]
    H = [NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, [1.2], 1.5, Float64[1.3])
    @test H == [1.5, 1.3] .* [2.0, 2.0]
    Hp = [NaN]
    MOI.eval_hessian_lagrangian_product(
        evaluator,
        Hp,
        [1.2],
        [1.2],
        1.5,
        Float64[1.3],
    )
    @test Hp == [1.5 * 2.0 * 1.2 + 1.3 * 2.0 * 1.2]
    return
end

function test_objective_quadratic_multivariate()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :($x^2 + $x * $y + $y^2))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2 + 1.2 * 2.3 + 2.3^2
    g = [NaN, NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
    @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (2, 1)]
    H = [NaN, NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, [1.2, 2.3], 1.5, Float64[])
    @test H == 1.5 .* [2.0, 2.0, 1.0]
    v = [0.3, 0.4]
    hv = [NaN, NaN]
    MOI.eval_hessian_lagrangian_product(
        evaluator,
        hv,
        [1.2, 2.3],
        v,
        1.5,
        Float64[],
    )
    @test hv ≈ 1.5 .* [2 1; 1 2] * v
    return
end

function test_objective_quadratic_multivariate_subexpressions()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    ex = Nonlinear.add_expression(model, :($x^2))
    ey = Nonlinear.add_expression(model, :($y^2))
    exy = Nonlinear.add_expression(model, :($ex + $x * $y))
    Nonlinear.set_objective(model, :($exy + $ey))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2 + 1.2 * 2.3 + 2.3^2
    g = [NaN, NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
    @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (2, 1)]
    H = [NaN, NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, [1.2, 2.3], 1.5, Float64[])
    @test H == 1.5 .* [2.0, 2.0, 1.0]
    v = [0.3, 0.4]
    hv = [NaN, NaN]
    MOI.eval_hessian_lagrangian_product(
        evaluator,
        hv,
        [1.2, 2.3],
        v,
        1.5,
        Float64[],
    )
    @test hv ≈ 1.5 .* [2 1; 1 2] * v
    return
end

function test_objective_ifelse_comparison()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(ifelse(1 <= $x <= 2, $x^2, $y^2)))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2
    @test MOI.eval_objective(evaluator, [2.2, 2.3]) == 2.3^2
    g = [NaN, NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
    @test g == [2 * 1.2, 0.0]
    MOI.eval_objective_gradient(evaluator, g, [2.2, 2.3])
    @test g == [0.0, 2 * 2.3]
    return
end

function test_objective_ifelse_logic()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(ifelse(1 <= $x && $x <= 2, $x^2, $y^2)))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.2, 2.3]) == 1.2^2
    @test MOI.eval_objective(evaluator, [2.2, 2.3]) == 2.3^2
    g = [NaN, NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
    @test g == [2 * 1.2, 0.0]
    MOI.eval_objective_gradient(evaluator, g, [2.2, 2.3])
    @test g == [0.0, 2 * 2.3]
    return
end

function test_objective_parameter()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    p = Nonlinear.add_parameter(model, 1.2)
    Nonlinear.set_objective(model, :($p * $x))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.eval_objective(evaluator, [1.3]) == 1.2 * 1.3
    g = [NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.3])
    @test g == [1.2]
    return
end

function test_objective_subexpression()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    input = :($x^2 + 1)
    expr = Nonlinear.add_expression(model, input)
    expr_2 = Nonlinear.add_expression(model, :($expr^2))
    Nonlinear.set_objective(model, :($expr_2^2))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    MOI.initialize(evaluator, [:Grad])
    @test MOI.eval_objective(evaluator, [1.3]) == ((1.3^2 + 1)^2)^2
    g = [NaN]
    MOI.eval_objective_gradient(evaluator, g, [1.3])
    @test g ≈ [2 * (1.3^2 + 1)^2 * (2 * (1.3^2 + 1)) * 2 * 1.3]
    return
end

function test_constraint_quadratic_univariate()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.add_constraint(model, :($x^2), MOI.LessThan(2.0))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    g = [NaN]
    x_val = [1.2]
    MOI.eval_constraint(evaluator, g, x_val)
    @test g == x_val .^ 2
    @test MOI.jacobian_structure(evaluator) == [(1, 1)]
    J = [NaN]
    MOI.eval_constraint_jacobian(evaluator, J, x_val)
    @test J == 2 .* x_val
    @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1)]
    H = [NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])
    @test H == 1.5 .* [2.0]
    return
end

function test_constraint_quadratic_multivariate()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.add_constraint(model, :($x^2 + $x * $y + $y^2), MOI.LessThan(2.0))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    g = [NaN]
    x_val = [1.2, 2.3]
    MOI.eval_constraint(evaluator, g, x_val)
    @test g == [x_val[1]^2 + x_val[1] * x_val[2] + x_val[2]^2]
    @test MOI.jacobian_structure(evaluator) == [(1, 1), (1, 2)]
    J = [NaN, NaN]
    MOI.eval_constraint_jacobian(evaluator, J, x_val)
    @test J == [2 * x_val[1] + x_val[2], x_val[1] + 2 * x_val[2]]
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (2, 1)]
    H = [NaN, NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])
    @test H == 1.5 .* [2.0, 2.0, 1.0]
    return
end

function test_constraint_quadratic_multivariate_subexpressions()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    ex = Nonlinear.add_expression(model, :($x^2))
    ey = Nonlinear.add_expression(model, :($y^2))
    exy = Nonlinear.add_expression(model, :($ex + $x * $y))
    Nonlinear.add_constraint(model, :($exy + $ey), MOI.LessThan(2.0))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    g = [NaN]
    x_val = [1.2, 2.3]
    MOI.eval_constraint(evaluator, g, x_val)
    @test g ≈ [x_val[1]^2 + x_val[1] * x_val[2] + x_val[2]^2]
    # Jacobian
    @test MOI.jacobian_structure(evaluator) == [(1, 1), (1, 2)]
    J = [NaN, NaN]
    MOI.eval_constraint_jacobian(evaluator, J, x_val)
    @test J ≈ [2 * x_val[1] + x_val[2], x_val[1] + 2 * x_val[2]]
    # Jv
    y, w = [NaN], [1.1, 2.2]
    MOI.eval_constraint_jacobian_product(evaluator, y, x_val, w)
    @test y ≈ [(2 * x_val[1] + x_val[2]) (x_val[1] + 2 * x_val[2])] * w
    # v'J
    y, w = [NaN, NaN], [1.1]
    MOI.eval_constraint_jacobian_transpose_product(evaluator, y, x_val, w)
    wJ = w' * [(2 * x_val[1] + x_val[2]) (x_val[1] + 2 * x_val[2])]
    @test y ≈ wJ[:]
    # Hessian-lagrangian
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (2, 1)]
    H = [NaN, NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])
    @test H ≈ 1.5 .* [2.0, 2.0, 1.0]
    return
end

function test_hessian_sparsity_registered_function()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    f(x, y) = x^2 + y^2
    function ∇f(g, x...)
        g[1] = 2x[1]
        g[2] = 2x[2]
        return
    end
    function ∇²f(H, x...)
        H[1, 1] = 2
        H[2, 2] = 2
        return
    end
    model = Nonlinear.Model()
    Nonlinear.register_operator(model, :f, 2, f, ∇f, ∇²f)
    Nonlinear.set_objective(model, :(f($x, $z) + $y^2))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y, z])
    @test :Hess in MOI.features_available(evaluator)
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (3, 3), (3, 1)]
    H = fill(NaN, 4)
    MOI.eval_hessian_lagrangian(evaluator, H, rand(3), 1.5, Float64[])
    @test H == 1.5 .* [2.0, 2.0, 2.0, 0.0]
    return
end

function test_hessian_sparsity_registered_rosenbrock()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    f(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
    function ∇f(g, x...)
        g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2
        g[2] = 200 * (x[2] - x[1]^2)
        return
    end
    function ∇²f(H, x...)
        @assert size(H) == (2, 2)
        H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2
        H[2, 1] = -400 * x[1]
        H[2, 2] = 200.0
        return
    end
    model = Nonlinear.Model()
    Nonlinear.register_operator(model, :rosenbrock, 2, f, ∇f, ∇²f)
    Nonlinear.set_objective(model, :(rosenbrock($x, $y)))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    @test :Hess in MOI.features_available(evaluator)
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    @test MOI.hessian_lagrangian_structure(evaluator) ==
          [(1, 1), (2, 2), (2, 1)]
    H = fill(NaN, 3)
    MOI.eval_hessian_lagrangian(evaluator, H, [1.0, 1.0], 1.5, Float64[])
    @test H == 1.5 .* [802.0, 200.0, -400.0]
    return
end

function test_hessian_registered_error()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    f(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
    function ∇f(g, x...)
        g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2
        g[2] = 200 * (x[2] - x[1]^2)
        return
    end
    function ∇²f(H, x...)
        H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2
        # Wrong index! Should be [2, 1]
        H[1, 2] = -400 * x[1]
        H[2, 2] = 200.0
        return
    end
    model = Nonlinear.Model()
    Nonlinear.register_operator(model, :rosenbrock, 2, f, ∇f, ∇²f)
    Nonlinear.set_objective(model, :(rosenbrock($x, $y)))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    H = fill(NaN, 3)
    @test_throws(
        ErrorException("Unable to access upper-triangular component: (1, 2)"),
        MOI.eval_hessian_lagrangian(evaluator, H, [1.0, 1.0], 1.5, Float64[]),
    )
    return
end

struct _ColoringGraph
    num_vertices::Int
    edges::Vector{Tuple{Int,Int}}
end

function to_adjlist(graph::_ColoringGraph)
    I = [i for (i, _) in graph.edges]
    J = [j for (_, j) in graph.edges]
    return Coloring.UndirectedGraph(I, J, graph.num_vertices)
end

function test_coloring_edge_free_graph()
    graph = _ColoringGraph(10, [])
    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 1
    return
end

function test_coloring_one_edge_graph()
    graph = _ColoringGraph(10, [(2, 4)])
    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 2
    @test color[2] != color[4]
    return
end

function test_coloring_two_edge_graph()
    graph = _ColoringGraph(10, [(2, 4), (2, 3)])
    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 2
    @test color[3] == color[4]
    return
end

function test_coloring_three_edge_graph()
    graph = _ColoringGraph(10, [(2, 4), (2, 3), (3, 4)])
    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 3
    # TODO: What is this testing?
    Coloring.recovery_preprocess(to_adjlist(graph), color, numcolors, Int[])
    return
end

function test_coloring_two_edge_three_vertex_graph()
    graph = _ColoringGraph(3, [(1, 3), (2, 3)])
    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 2
    return
end

function test_coloring_four_edge_four_vertex_graph()
    graph = _ColoringGraph(4, [(1, 2), (2, 3), (3, 4), (4, 1)])
    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
    @test numcolors == 3
    return
end

function test_coloring_topological_sort()
    # graph = _ColoringGraph(6, [(1, 2), (1, 3), (1, 6), (2, 4), (2, 5)])
    vec = [3, 6, 2, 1, 4, 5, 1, 2, 2, 1]
    offset = [1, 4, 7, 8, 9, 10, 11]
    v = Coloring.reverse_topological_sort_by_dfs(vec, offset, 6, zeros(Int, 6))
    @test v[1] == [3, 6, 4, 5, 2, 1]
    @test v[2] == [0, 1, 1, 2, 2, 1]
    return
end

function test_coloring_end_to_end_hessian_coloring_and_recovery()
    I, J, rinfo = Coloring.hessian_color_preprocess(Set([(1, 2)]), 2)
    R = Coloring.seed_matrix(rinfo)
    Coloring.prepare_seed_matrix!(R, rinfo)
    @test I == [1, 2, 2]
    @test J == [1, 2, 1]
    @test R == [1.0 0.0; 0.0 1.0]
    hess = [3.4 2.1; 2.1 1.3]
    matmat = hess * R
    V = zeros(3)
    Coloring.recover_from_matmat!(V, matmat, rinfo, zeros(3))
    @test V == [3.4, 1.3, 2.1]
    return
end

function test_derivatives()
    a = MOI.VariableIndex(1)
    b = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(sin($a^2) + cos($b * 4) / 5 - 2.0))
    evaluator =
        Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [a, b])
    MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
    x = [2.0, 3.0]
    @test MOI.eval_objective(evaluator, x) ==
          sin(x[1]^2) + cos(x[2] * 4) / 5 - 2.0
    g = [NaN, NaN]
    MOI.eval_objective_gradient(evaluator, g, x)
    @test g ≈ [2 * x[1] * cos(x[1]^2), -4 * sin(x[2] * 4) / 5]
    @test MOI.hessian_lagrangian_structure(evaluator) == [(1, 1), (2, 2)]
    H = [NaN, NaN]
    MOI.eval_hessian_lagrangian(evaluator, H, x, 1.5, Float64[])
    H_exact = [
        -4 * x[1]^2 * sin(x[1]^2) + 2 * cos(x[1]^2),
        -4 / 5 * 4 * cos(x[2] * 4),
    ]
    @test H == 1.5 .* H_exact
    return
end

function test_NLPBlockData()
    model = Nonlinear.Model()
    x = MOI.VariableIndex(1)
    Nonlinear.add_constraint(model, :($x - 1), MOI.LessThan(0.0))
    Nonlinear.add_constraint(model, :($x - 2), MOI.GreaterThan(0.0))
    Nonlinear.add_constraint(model, :($x - 3), MOI.EqualTo(0.0))
    Nonlinear.add_constraint(model, :($x), MOI.Interval(4.0, 5.0))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
    block = MOI.NLPBlockData(evaluator)
    @test block.constraint_bounds == [
        MOI.NLPBoundsPair(-Inf, 0.0),
        MOI.NLPBoundsPair(0.0, Inf),
        MOI.NLPBoundsPair(0.0, 0.0),
        MOI.NLPBoundsPair(4.0, 5.0),
    ]
    @test block.has_objective == false
    return
end

function test_linearity()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    variables = Dict(x => 1, y => 2, z => 3)
    function _test_linearity(input, test_value, IJ = [], indices = [])
        model = Nonlinear.Model()
        ex = Nonlinear.add_expression(model, input)
        expr = model[ex]
        adj = Nonlinear.adjacency_matrix(expr.nodes)
        nodes = ReverseAD._replace_moi_variables(expr.nodes, variables)
        ret = ReverseAD._classify_linearity(nodes, adj, ReverseAD.Linearity[])
        @test ret[1] == test_value
        indexed_set = Coloring.IndexedSet(100)
        edge_list = ReverseAD._compute_hessian_sparsity(
            nodes,
            adj,
            ret,
            indexed_set,
            Set{Tuple{Int,Int}}[],
            Vector{Int}[],
        )
        if ret[1] != ReverseAD.NONLINEAR
            @test length(edge_list) == 0
        elseif length(IJ) > 0
            @test IJ == edge_list
        end
        if length(indices) > 0
            empty!(indexed_set)
            ReverseAD._compute_gradient_sparsity!(indexed_set, nodes)
            ix = sort(collect(indexed_set))
            @test indices == ix
        end
        return
    end
    _test_linearity(
        :(sin($x^2) + cos($y * 4) - 2.0),
        ReverseAD.NONLINEAR,
        Set([(2, 2), (1, 1)]),
        [1, 2],
    )
    _test_linearity(:(3 * 4 * ($x + $y)), ReverseAD.LINEAR)
    _test_linearity(
        :($z * $y),
        ReverseAD.NONLINEAR,
        Set([(3, 2), (3, 3), (2, 2)]),
        [2, 3],
    )
    _test_linearity(:(3 + 4), ReverseAD.CONSTANT)
    _test_linearity(:(sin(3) + $x), ReverseAD.LINEAR)
    _test_linearity(
        :(cos($z) * sin(3) + $x),
        ReverseAD.NONLINEAR,
        Set([(3, 3)]),
        [1, 3],
    )
    _test_linearity(:($x - 3 * $y), ReverseAD.LINEAR)
    _test_linearity(:(-$x), ReverseAD.LINEAR)
    _test_linearity(:(+$x), ReverseAD.LINEAR)
    _test_linearity(
        :($x^$y),
        ReverseAD.NONLINEAR,
        Set([(2, 2), (1, 1), (2, 1)]),
        [1, 2],
    )
    _test_linearity(:($x / 3 + $y), ReverseAD.LINEAR)
    _test_linearity(
        :(3 / ($x * $y)),
        ReverseAD.NONLINEAR,
        Set([(2, 2), (1, 1), (2, 1)]),
        [1, 2],
    )
    _test_linearity(:(1 / ($x + 3)), ReverseAD.NONLINEAR, Set([(1, 1)]), [1])
    _test_linearity(
        :(ifelse($x <= 1, $x, $y)),
        ReverseAD.PIECEWISE_LINEAR,
        Set([]),
        [],
    )
    _test_linearity(
        :(ifelse($x <= 1, $x^2, $y)),
        ReverseAD.NONLINEAR,
        Set([(1, 1)]),
        [1, 2],
    )
    _test_linearity(:(ifelse(1 <= 1, 2, 3)), ReverseAD.CONSTANT)
    _test_linearity(
        :(1 / ifelse($x < 1, $x, 0)),
        ReverseAD.NONLINEAR,
        Set([(1, 1)]),
        [1],
    )
    return
end

function test_dual_forward()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    function _test_dual_forward(input, x_input, test_value)
        model = Nonlinear.Model()
        Nonlinear.set_objective(model, input)
        evaluator = Nonlinear.Evaluator(
            model,
            Nonlinear.SparseReverseMode(),
            MOI.VariableIndex[x, y],
        )
        MOI.initialize(evaluator, [:Grad])
        ∇f = fill(NaN, 2)
        MOI.eval_objective_gradient(evaluator, ∇f, x_input)
        @test isapprox(∇f, test_value, atol = 1e-10)
        return
    end
    _test_dual_forward(
        :(sin($x^2) + cos($y * 4) / 5 - 2.0),
        [1.0, 2.0],
        [1.0806046117362795, -0.7914865972987055],
    )
    _test_dual_forward(
        :(sin($x^$y) + cos($y * 4) / 5 - 2.0),
        [1.0, 2.0],
        [1.0806046117362795, -0.7914865972987055],
    )
    _test_dual_forward(
        :(sin($x^3) + cos($x * $y * 4) / 5 - 2.0),
        [1.0, 0.0],
        [1.6209069176044193, 0.0],
    )
    _test_dual_forward(
        :($x * $y),
        [3.427139283036299e-206, 1.0],
        [1.0, 3.427139283036299e-206],
    )
    return
end

function test_gradient_registered_function()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    model = Nonlinear.Model()
    f(x, y) = (1 / 3) * y^3 - 2x^2
    function ∇f(g, x, y)
        g[1] = -4x
        g[2] = y^2
        return
    end
    Nonlinear.register_operator(model, :Φ, 2, f, ∇f)
    Nonlinear.register_operator(model, :c, 1, cos, x -> -sin(x), x -> -cos(x))
    Nonlinear.set_objective(model, :(Φ($y, $x - 1) * c($z)))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x, y, z],
    )
    MOI.initialize(evaluator, [:Grad])
    ∇f = fill(NaN, 3)
    x_input = [2.0, 3.0, 4.0]
    @test MOI.eval_objective(evaluator, x_input) == f(3.0, 2.0 - 1) * cos(4.0)
    MOI.eval_objective_gradient(evaluator, ∇f, [2.0, 3.0, 4.0])
    @test ∇f ≈ [
        cos(x_input[3]) * (x_input[1] - 1)^2,
        -4cos(x_input[3]) * x_input[2],
        -sin(x_input[3]) * f(x_input[2], x_input[1] - 1),
    ]
    return
end

function test_gradient_jump_855()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(
        model,
        :(ifelse($x <= 3.0, ($x - 2.0)^2, 2 * log($x - 2.0) + 1.0)),
    )
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    ∇f = fill(NaN, 1)
    MOI.eval_objective_gradient(evaluator, ∇f, [-1.0])
    @test ∇f ≈ [-6.0]
    MOI.eval_objective_gradient(evaluator, ∇f, [2.0])
    @test ∇f ≈ [0.0]
    return
end

function test_gradient_abs()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(abs($x)))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    ∇f = fill(NaN, 1)
    MOI.eval_objective_gradient(evaluator, ∇f, [2.0])
    @test ∇f ≈ [1.0]
    MOI.eval_objective_gradient(evaluator, ∇f, [-2.0])
    @test ∇f ≈ [-1.0]
    return
end

function test_gradient_trig()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(sin($x^2) + cos($y * 4) / 5 - 2.0))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x, y],
    )
    MOI.initialize(evaluator, [:Grad])
    ∇f = fill(NaN, 2)
    MOI.eval_objective_gradient(evaluator, ∇f, [2.0, 3.0])
    @test ∇f ≈ [2 * 2.0 * cos(2.0^2), -4 * sin(3.0 * 4) / 5]
    return
end

function test_gradient_logical()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :($x > 0.5 && $x < 0.9))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    @test MOI.eval_objective(evaluator, [1.5]) == 0.0
    ∇f = fill(NaN, 1)
    MOI.eval_objective_gradient(evaluator, ∇f, [1.5])
    @test iszero(∇f)
    return
end

function test_gradient_ifelse()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(ifelse($x >= 0.5 || $x < 0.1, $x, 5)))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    @test MOI.eval_objective(evaluator, [1.5]) == 1.5
    ∇f = fill(NaN, 1)
    MOI.eval_objective_gradient(evaluator, ∇f, [1.5])
    @test ∇f ≈ [1.0]
    @test MOI.eval_objective(evaluator, [-0.1]) == -0.1
    MOI.eval_objective_gradient(evaluator, ∇f, [-0.1])
    @test ∇f ≈ [1.0]
    @test MOI.eval_objective(evaluator, [0.2]) == 5
    MOI.eval_objective_gradient(evaluator, ∇f, [0.2])
    @test ∇f ≈ [0.0]
    return
end

function test_gradient_sqrt_nan()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :(sqrt($x)))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    @test isnan(MOI.eval_objective(evaluator, [-1.5]))
    ∇f = fill(Inf, 1)
    MOI.eval_objective_gradient(evaluator, ∇f, [-1.5])
    @test isnan(∇f[1])
    return
end

function test_gradient_variable_power()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    model = Nonlinear.Model()
    Nonlinear.set_objective(model, :((1 / $x)^$y - $z))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x, y, z],
    )
    MOI.initialize(evaluator, [:Grad])
    x_input = [2.5, 3.5, 1.0]
    @test MOI.eval_objective(evaluator, x_input) == (1 / 2.5)^3.5 - 1.0
    ∇f = fill(Inf, 3)
    MOI.eval_objective_gradient(evaluator, ∇f, x_input)
    @test ∇f ≈ [
        -x_input[2] * x_input[1]^(-x_input[2] - 1),
        -((1 / x_input[1])^x_input[2]) * log(x_input[1]),
        -1,
    ]
    return
end

function test_single_parameter()
    x = MOI.VariableIndex(1)
    model = Nonlinear.Model()
    p = Nonlinear.add_parameter(model, 105.2)
    Nonlinear.set_objective(model, :($p))
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x],
    )
    MOI.initialize(evaluator, [:Grad])
    @test MOI.eval_objective(evaluator, [-0.1]) == 105.2
    return
end

function test_gradient_nested_subexpressions()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    model = Nonlinear.Model()
    ex1 = Nonlinear.add_expression(model, :(sin($x^2) + cos($y * 4) / 5 - 2.0))
    ex2 = Nonlinear.add_expression(model, :($ex1))
    Nonlinear.set_objective(model, ex2)
    evaluator = Nonlinear.Evaluator(
        model,
        Nonlinear.SparseReverseMode(),
        MOI.VariableIndex[x, y],
    )
    MOI.initialize(evaluator, [:Grad])
    ∇f = fill(NaN, 2)
    MOI.eval_objective_gradient(evaluator, ∇f, [2.0, 3.0])
    @test ∇f ≈ [2 * 2.0 * cos(2.0^2), -4 * sin(3.0 * 4) / 5]
    return
end

function _dense_hessian(hessian_sparsity, V, n)
    I = [i for (i, _) in hessian_sparsity]
    J = [j for (_, j) in hessian_sparsity]
    raw = SparseArrays.sparse(I, J, V, n, n)
    return Matrix(
        raw + raw' -
        SparseArrays.sparse(LinearAlgebra.diagm(0 => LinearAlgebra.diag(raw))),
    )
end

# This covers the code that computes Hessians in odd chunks of Hess-vec
# products.
function test_odd_chunks_Hessian_products()
    for i in 1:18
        _test_odd_chunks_Hessian_products(i)
    end
    return
end

function _test_odd_chunks_Hessian_products(N)
    model = Nonlinear.Model()
    x = MOI.VariableIndex.(1:N)
    Nonlinear.set_objective(model, Expr(:call, :*, x...))
    evaluator = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), x)
    MOI.initialize(evaluator, [:Hess])
    hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator)
    V = zeros(length(hessian_sparsity))
    values = ones(N)
    MOI.eval_hessian_lagrangian(evaluator, V, values, 1.0, Float64[])
    H = _dense_hessian(hessian_sparsity, V, N)
    @test H ≈ (ones(N, N) - LinearAlgebra.diagm(0 => ones(N)))
    values[1] = 0.5
    MOI.eval_hessian_lagrangian(evaluator, V, values, 1.0, Float64[])
    H = _dense_hessian(hessian_sparsity, V, N)
    H_22 = (ones(N - 1, N - 1) - LinearAlgebra.diagm(0 => ones(N - 1))) / 2
    @test H ≈ [0 ones(N - 1)'; ones(N - 1) H_22]
    return
end

end  # module

TestReverseAD.runtests()