text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# diffusion.jl contains several functions focused on the treatment of infrared spectra along diffusion profiles in minerals
#
#
#############################################################################
"""
The 'IRdataprep' function allows to obtain the frequency, distance and y vectors in a region of interest, ready to input in an optimisation algorithm.

   IRdataprep(data,distance_step,start_sp,stop_sp,low_x,high_x,norm_low_x,norm_high_x)
with
data (array, Float64): an array containing in the first column the x values and in the subsequent columns the correcponding intensities;
distance_step (Float64): the steps between each spectrum along the diffusion profile, in metres;
start_sp and stop_sp (Int): the starting and stopping spectra of interest (warning: number of data columns - 1, as data also contains the x axis such that size(data)[2] = number of spectras + 1);
low_x,high_x (Float64): frequencies between which the signal is of interest, defined as low_x < high_x;
norm_low_x and norm_high_x (Float64): frequencies between which the signal must be integrated for mormalising the final signals.

This function first checks that the x values are increasing, and correct that if this is not the case.

Then, it selects the data in the region of interest, delimited by frequencies defined as low_x < high_x.
Normalisation to the area between the norm_low_x and norm_high_x frequency is also performed (norm_low_x < norm_high_x)

The code returns the x and the y arrays, plus arrays for optimisation: x_input_fit contains the frequency-distance couples and y_input_fit their corresponding y values.
"""

function IRdataprep(data::Array{Float64},distance_step::Float64,start_sp::Int,stop_sp::Int,low_x::Float64,high_x::Float64,norm_low_x::Float64,norm_high_x::Float64,ese_low_x::Float64,ese_high_x::Float64)

   if data[end,1] <= data[1,1]
      data=data[end:-1:1,:]
      data[data[:,:] .== 0]= 1e-20# to avoid pure 0 values
   end

   x::Array{Float64} = data[find(low_x .< data[:,1] .< high_x),1]
   x_sili::Array{Float64} = data[find(norm_low_x .< data[:,1] .< norm_high_x),1]

   if stop_sp .>= size(data)[2]
        y::Array{Float64} = data[find(low_x .< data[:,1] .< high_x),start_sp:end]
        y_sili::Array{Float64} = data[find(norm_low_x .< data[:,1] .< norm_high_x),start_sp:end]
        ese::Array{Float64} =  std(data[find(ese_low_x .< data[:,1] .< ese_high_x),start_sp:end],1) #relative error
        y_ese_r::Array{Float64} = ese[1,:]./y[:,:]
   else
        y = data[find(low_x .< data[:,1] .< high_x),start_sp:stop_sp]
        y_sili = data[find(norm_low_x .< data[:,1] .< norm_high_x),start_sp:stop_sp]
        ese =  std(data[find(ese_low_x .< data[:,1] .< ese_high_x),start_sp:stop_sp],1) #relative error
        y_ese_r = ese[1,:]./y[:,:]

   end

   for i = 1:size(y)[2]
       y[:,i] = y[:,i]/trapz(x_sili[:,1],y_sili[:,i]); # integration of silicate bands for normalisation
   end

   # constructing the good distance vector + x and Y associated values
   steps::Float64 = 0.0
   x_input_fit::Array{Float64} = zeros(length(x),2)
   y_input_fit::Array{Float64} = zeros(length(x),1)
   ese_input_fit::Array{Float64} = zeros(length(x),1)
   for i = 1:size(y)[2]
      x_temp::Array{Float64} = zeros(length(x),2)
      if i == 1
         x_input_fit[:,1] = steps
         x_input_fit[:,2] = x
         y_input_fit[:,1] = y[:,i]
         ese_input_fit[:,1]  = y[:,i] .* y_ese_r[:,i]
      else
         x_temp[:,1] = steps
         x_temp[:,2] = x

         x_input_fit = vcat(x_input_fit,x_temp)
         y_input_fit = vcat(y_input_fit,y[:,i])
         ese_input_fit = vcat(ese_input_fit,y[:,i] .* y_ese_r[:,i])
      end
   steps = steps + distance_step
   end

    return x, y, y_ese_r, x_input_fit, y_input_fit, ese_input_fit
end

"""
The `peak_diffusion` function allows to fit/plot gaussian peaks for which the amplitude is defined by a 1D diffusion law
Call as:
    peak_diffusion(g_c0::Array{Float64,1},g_c1::Array{Float64,1},g_D::Array{Float64,1}, g_freq::Array{Float64,1}, g_hwhm::Array{Float64,1},x::Array{Float64,2},time::Float64)

where:
g_c0 is the vector of concentration imposed, at the diffusive bondary (a.u., float64 number);
g_c1 is the vector of initial concentration in the host (a.u., float64 number);
g_D is the vector of diffusivity coefficient, in log units
g_freq is the vector containing the frequencies fo the peaks
g_hwhm is the vector containing the hwhm of the peaks
x is an array containing the distances and associated frequencies
time is the duration of the run in seconds

The function returns a vector of the values of the peak intensity at the input frequencies.

The current version assumes the peak being a Gaussian peak. Updates will integrate Lorentzian and Pseudo-Voigt peaks.
"""
function peak_diffusion(g_c0::Array{Float64,1},g_c1::Array{Float64,1},g_D::Array{Float64,1}, g_freq::Array{Float64,1}, g_hwhm::Array{Float64,1},x::Array{Float64,2},time::Float64)
    segments = zeros(size(x)[1],size(g_c0)[1])
    for i = 1:size(g_c0)[1]
        segments[:,i] = ((g_c1[i] - g_c0[i]) .* erfc(x[:,1]./(2. .* sqrt( 10.^g_D[i] .* time))) + g_c0[i]) .*exp(-log(2) .* ((x[:,2]-g_freq[i])./g_hwhm[i]).^2)
    end
    return sum(segments,2), segments
end

"""
The global model for using with the Julia.optim package. Not necessarily up to date because I prefer using the JuMP way.
   model(p,x_fit,temps;compo_out = 0)
"""
function model(p::Array{Float64,1},x_fit::Array{Float64,2},temps::Float64;compo_out::Int = 0)
    number_gauss = length(p)/5
    components = zeros(length(x_fit[:,2]),round(Int8,number_gauss))
    compteur = 1
    for i = 1:5:length(p)
        components[:,compteur] = peak_diffusion(x_fit[:,1],temps,x_fit[:,2],p[i],p[i+1],p[i+2],p[i+3],p[i+4])
        compteur = compteur + 1
    end
    if compo_out == 0
        return sum(components,2)
    else
        return components
    end
end


text:
const PRIOR_PROPOSAL_NAME = "prior"
const OPTIMAL_PROPOSAL_NAME = "optimal"

const COL_PROPOSAL_NAME = :proposal_name
const COL_NUM_SWEEPS = :num_sweeps
const COL_NUM_PARTICLES = :num_particles
const COL_AIDE_ESTIMATE = :aide_estimate
const COL_AIDE_STDERR = :aide_stderr
const COL_NUM_CLUSTERS_ESTIMATE = :num_clusters_estimate
const COL_NUM_CLUSTERS_STDERR = :num_clusters_stderr

ext: jl


text:
import Random
using Random: MersenneTwister
using Distributions
using Test
using StaticArrays
using LinearAlgebra
using Unitful

μ = rand()
x = rand()
σ = rand()
Σ = σ*σ'

# Check type conversions
GFloat = Gaussian{Vector{Float64}, Matrix{Float64}}
v = GFloat[Gaussian([1.0], Matrix(1.0I, 2, 2)),
           Gaussian(SVector(1.0), @SMatrix [1.0 0.0; 0.0 1.0])]
@test mean.(v) == [[1.0], [1.0]]

G = Gaussian([1.0], Matrix(1.0I, 2, 2))
m, K = G
@test m === mean(G)
@test K === cov(G)

p = pdf(Normal(μ, √Σ), x)
@test pdf(Gaussian(μ, Σ), x) ≈ p
@test pdf(Gaussian(μ, Σ*I), x) ≈ p
@test pdf(Gaussian([μ], [σ]*[σ]'), [x]) ≈ p

@test pdf(Gaussian((@SVector [μ]), @SMatrix [Σ]), @SVector [x]) ≈ p

for d in 1:3
    μ = rand(d)
    x = rand(d)
    σ = tril(rand(d,d))
    Σ = σ*σ'
    p = pdf(MvNormal(μ, Σ), x)

    @test pdf(Gaussian(μ, Σ), x) ≈ p
    @test pdf(Gaussian(μ, PSD(σ)), x) ≈ p
    @test pdf(Gaussian(SVector{d}(μ), SMatrix{d,d}(Σ)), x) ≈ p
    @test pdf(Gaussian(SVector{d}(μ), PSD(SMatrix{d,d}(σ))), x) ≈ p
end

for d in 1: 3
    local μ = rand(d)
    local x = rand(d)
    local σ = rand()
    local Σ = Matrix(1.0I, d, d).*σ^2
    local p = pdf(MvNormal(μ, Σ), x)

    @test pdf(Gaussian(μ, σ^2*I), x) ≈ p
    @test pdf(Gaussian(SVector{d}(μ), SDiagonal(σ^2*ones(SVector{d}))), x) ≈ p
    @test pdf(Gaussian(SVector{d}(μ), SMatrix{d,d}(Σ)), x) ≈ p
end

@test rand(Random.GLOBAL_RNG, Gaussian(1.0, 0.0)) == 1.0
@test mean(rand(MersenneTwister(1), Gaussian([1., 2], Matrix(1.0I, 2, 2)), 100000)) ≈ 1.5 atol=0.02
@test mean(rand(MersenneTwister(1), Gaussian(1.5u"m", 2.0u"m^2"), 100000)) ≈ 1.5u"m" atol=0.02u"m"
@test std(Gaussian(1.0u"m", 4.0u"m^2")) == 2.0u"m"

@test rand(Gaussian(1.0,0.0)) == 1.0
Random.seed!(5)
x = randn()
Random.seed!(5)
@test rand(Gaussian(1.0,2.0)) == 1.0 .+ sqrt(2)*x

Random.seed!(5)
@test rand(Gaussian(1.0,2.0), (1,)) == [1.0 + sqrt(2)*x]

g = Gaussian([1,2], Matrix(1.0I, 2, 2))
@test mean(g + [10,10]) == [11,12]
@test (g - [10,10]) + [10,10] == g
@test mean(g + [10, 20]) == [11,22]
m = [1.0 2.0; 0.0 2.0]
@test cov(m * g) == m * m'

using GaussianDistributions: ⊕
@test cov(m*g ⊕ g) == cov(m*g) + cov(g)
@test mean(m*g ⊕ g) == mean(m*g) + mean(g)

@test cov(g ⊕ m*mean(g)) == cov(g)
@test mean(m*mean(g) ⊕ g) == mean(m*g) + mean(g)

g2 = Gaussian(SVector(1.000001, 2.0), [1.00001 0.0; 0.0 0.999999])
@test !(g ≈ g2)
@test g ≈ g2 rtol=0.0001


text:
function loadchar(face::FTFont, c::Char)
    err = FT_Load_Char(face, c, FT_LOAD_RENDER)
    check_error(err, "Could not load char to render.")
end

function renderface(face::FTFont, c::Char, pixelsize::Integer)
    set_pixelsize(face, pixelsize)
    loadchar(face, c)
    glyph = unsafe_load(face.glyph)
    @assert glyph.format == FreeType.FT_GLYPH_FORMAT_BITMAP
    return glyphbitmap(glyph.bitmap), FontExtent(glyph.metrics)
end

function glyphbitmap(bitmap::FreeType.FT_Bitmap)
    @assert bitmap.pixel_mode == FreeType.FT_PIXEL_MODE_GRAY
    bmp = Matrix{UInt8}(undef, bitmap.width, bitmap.rows)
    row = bitmap.buffer
    if bitmap.pitch < 0
        row -= bitmap.pitch * (rbmpRec.rows - 1)
    end
    for r in 1:bitmap.rows
        src = unsafe_wrap(Array, row, bitmap.width)
        bmp[:, r] = src
        row += bitmap.pitch
    end
    return bmp
end

one_or_typemax(::Type{T}) where {T<:Union{Real,Colorant}} = T<:Integer ? typemax(T) : oneunit(T)

"""
    renderstring!(img::AbstractMatrix, str::String, face, pixelsize, y0, x0;
    fcolor=one_or_typemax(T), bcolor=zero(T), halign=:hleft, valign=:vbaseline) -> Matrix

Render `str` into `img` using the font `face` of size `pixelsize` at coordinates `y0,x0`.

# Arguments
* `y0,x0`: origin is in upper left with positive `y` going down
* `fcolor`: foreground color; typemax(T) for T<:Integer, otherwise one(T)
* `bcolor`: background color; set to `nothing` for transparent
* `halign`: :hleft, :hcenter, or :hright
* `valign`: :vtop, :vcenter, :vbaseline, or :vbottom
"""
function renderstring!(
        img::AbstractMatrix{T}, str::String, face::FTFont, pixelsize::Union{Int, Tuple{Int, Int}}, y0, x0;
        fcolor::T = one_or_typemax(T), bcolor::Union{T,Nothing} = zero(T),
        halign::Symbol = :hleft, valign::Symbol = :vbaseline
    ) where T<:Union{Real,Colorant}

    if pixelsize isa Tuple
        @warn "using tuple for pixelsize is deprecated, please use one integer"
        pixelsize = pixelsize[1]
    end

    set_pixelsize(face, pixelsize)

    bitmaps = Vector{Matrix{UInt8}}(undef, lastindex(str))
    metrics = Vector{FontExtent{Int}}(undef, lastindex(str))
    # ymin and ymax w.r.t the baseline
    ymin = ymax = sumadvancex = 0

    for (istr, char) = enumerate(str)
        bitmap, metric_float = renderface(face, char, pixelsize)
        metric = round.(Int, metric_float)
        bitmaps[istr] = bitmap
        metrics[istr] = metric

        # scale of glyph (size of bitmap)
        w, h = metric.scale
        # offset between glyph origin and bitmap top left corner
        bx, by = metric.horizontal_bearing
        ymin = min(ymin, by - h)
        ymax = max(ymax, by)
        sumadvancex += metric.advance[1]
    end

    px = x0 - (halign == :hright ? sumadvancex : halign == :hcenter ? sumadvancex >> 1 : 0)
    py = y0 + (
        valign == :vtop ? ymax : valign == :vbottom ? ymin :
        valign == :vcenter ? (ymax - ymin) >> 1 + ymin : 0
    )

    bitmapmax = typemax(eltype(bitmaps[1]))

    imgh, imgw = size(img)
    if bcolor !== nothing
        img[
            clamp(py-ymax+1, 1, imgh) : clamp(py-ymin-1, 1, imgh),
            clamp(px-1, 1, imgw) : clamp(px+sumadvancex-1, 1, imgw)
        ] .= bcolor
    end

    local prev_char::Char
    for (istr, char) = enumerate(str)
        w, h = metrics[istr].scale
        bx, by = metrics[istr].horizontal_bearing

        if istr == 1
            prev_char = char
        else
            kx, ky = map(x-> round(Int, x), kerning(prev_char, char, face))
            px += kx
        end

        # trim parts of glyph images that are outside the destination
        cliprowlo, cliprowhi = max(0, -(py-by)), max(0, py - by + h - imgh)
        clipcollo, clipcolhi = max(0, -bx-px),   max(0, px + bx + w - imgw)

        if bcolor === nothing
            for row = 1+cliprowlo : h-cliprowhi, col = 1+clipcollo : w-clipcolhi
                bitmaps[istr][col,row]==0 && continue
                c1 = bitmaps[istr][col,row] / bitmapmax * fcolor
                img[row+py-by, col+px+bx] = T <: Integer ? round(T, c1) : T(c1)
            end
        else
            for row = 1+cliprowlo : h-cliprowhi, col = 1+clipcollo : w-clipcolhi
                bitmaps[istr][col, row] == 0 && continue
                w1 = bitmaps[istr][col, row] / bitmapmax
                c1 = w1 * fcolor
                c0 = (1.0 - w1) * bcolor
                img[row + py - by, col + px + bx] = T <: Integer ? round(T, c1 + c0) : T(c1 + c0)
            end
        end
        px += metrics[istr].advance[1]
    end
    return img
end


text:
using VertexEliminationOrder
using Test
using Graphs

ENV["JULIA_DEBUG"]=VertexEliminationOrder

G = SimpleGraph{Int64}(12)
# 1 to 5 is a clique
for i in 1:5
    for j in (i+1):5
        add_edge!(G, i, j)
    end
end
# 6 and 7 are midle nodes
add_edge!(G, 5, 6)
add_edge!(G, 6, 7)
add_edge!(G, 6, 8)
# 8 to 12 is a tree
add_edge!(G, 8, 9)
add_edge!(G, 8, 10)
add_edge!(G, 9, 11)
add_edge!(G, 9, 12)

@info "Testing iterative_dissection with custom graph"
res = iterative_dissection(G)
@info res
res_expected = treewidth_by_elimination!(G, res[1])
# 6 should be separator but t seems that it depends on the random s and
@info res_expected
@test abs(res[2] - res_expected) <= 1
@test res[2] >=  4# 4 is the tw but it returns an upper bound

for n in 2:20
    g = square_lattice_graph(n)
    @info "Testing iterative_dissection with square lattice graph n = $n"
    tmp = iterative_dissection(g)
    @info tmp
    tmp_expected = treewidth_by_elimination!(g, tmp[1])
    @info tmp_expected
    @test abs(tmp[2] - tmp_expected) <= 1
end


text:
const program = parse.(Int, split(readchomp("input"), ','))

function run!(program::AbstractVector)
    pos = 1
    while pos < length(program) && program[pos] != 99
        opcode = program[pos]
        if opcode == 1
            op = +
        elseif opcode == 2
            op = *
        else
            error("unexpected opcode $opcode")
        end
        p1, p2, out = program[pos+1:pos+3] .+ 1
        program[min(out, end)] = op(program[p1], program[p2])
        pos += 4
    end
    return program
end

function output(input, noun, verb)
    program = copy(input)
    program[2] = noun
    program[3] = verb
    run!(program)
    return first(program)
end

part1(input) = println("Part 1 answer: ", output(input, 12, 2))

function part2(input, target)
    bound = 99
    noun = verb = -1
    for outer noun = 0:bound, outer verb = 0:bound
        output(input, noun, verb) == target && break
    end
    println("Part 2 answer: ", 100 * noun + verb)
end

part1(program)
part2(program, 19690720)

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
# # [How to port rules generated by SparseADRules to Zygote](@id lowrank_svd)
#
# In this demo we'll show how to insert SparseADRules's sparse gradient implementation to boost Zygote's gradient in lowrank svd.
using Zygote
using SparseArrays, LinearAlgebra
using SparseADRules:low_rank_svd
import FiniteDifferences
using BenchmarkTools
using ChainRulesCore, ChainRulesTestUtils
include(joinpath(@__DIR__, "..", "test", "testutils.jl"))

function svd_loss(A::AbstractSparseMatrix{T}, z::Vector{T}, r::Int) where T
    U, S, Vt = low_rank_svd(A, r)
    residual_mat = U * Diagonal(S) * Vt - U * Diagonal(z) * Vt
    return tr(residual_mat'*residual_mat) # square of Frobenius norm
end

# generate test data
r = 10
A = sprand(100, r, 0.2) * Diagonal(rand(r)) * sprand(r, 50, 0.2)
z = rand(r)

@btime grad_A = Zygote.gradient(A -> svd_loss(A, z, r), $A)
grad_A = Zygote.gradient(A -> svd_loss(A, z, r), A)[1]

grad_fA = FiniteDifferences.grad(FiniteDifferences.central_fdm(5, 1), A -> svd_loss(A, z, r), A)[1]

using Test
@testset "test AD on lowrank svd" begin
    @test grad_A ≈ grad_fA
end


text:
using NBodyIPs, JuLIP, BenchmarkTools, StaticArrays
using Test

using NBodyIPs: BondLengthDesc, BondAngleDesc,
                transform, transform_d, fcut, fcut_d,
                invariants, invariants_d, invariants_ed
using NBodyIPs.Polys: NBPoly

using LinearAlgebra: norm
using Printf

profile = false

if profile
   nbasis = [0, 10, 50, 100, 300]
else
   nbasis = [0, 10, 10, 10, 10]
end

const Ir = @SVector [1,2,3]
const Iθ = @SVector [4,5,6]
all_invs(desc::BondLengthDesc, r) = vcat(invariants(desc, r)...)
jac_all_invs(desc::BondLengthDesc, r) = hcat(vcat(invariants_ed(desc, r)[3:4]...)...)'
all_invs(desc::BondAngleDesc, rθ) = vcat(invariants(desc, (rθ[Ir], rθ[Iθ]))...)
jac_all_invs(desc::BondAngleDesc, rθ) = hcat(vcat(invariants_ed(desc, (rθ[Ir], rθ[Iθ]))[3:4]...)...)'

function fdjac(desc, rθ; h=1e-5)
   v = Vector(rθ)
   f = all_invs(desc, rθ)
   J = zeros(length(f), length(v))
   for i = 1:length(v)
      v[i] += h
      fh = all_invs(desc, SVector(v...))
      J[:, i] = (fh-f) / h
      v[i] -= h
   end
   return J
end

##

for Desc in [BondLengthDesc, BondAngleDesc]
   # Desc = BondAngleDesc
   # Desc = BondLengthDesc

   println("------------------------------------------------")
   println(" Testing Polys based on $(Desc) descriptors")
   println("------------------------------------------------")

   println("Setting up the test systems ...")
   r0 = rnn(:Cu)
   TRANSFORM = ExpTransform(3, r0)
   at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)
   rcut3 = 3.1 * r0
   D3 = Desc(TRANSFORM, CosCut(0.66*rcut3, rcut3) )
   rcut4 = 2.1 * r0
   D4 = Desc(TRANSFORM, CosCut(0.66*rcut4, rcut4) )
   rcut5 = 1.5 * r0
   D5 = Desc(TRANSFORM, CosCut(0.66*rcut5, rcut5) )
   DD = [nothing, D3, D3, D4, D5]

   random_nbody(N, ntup) = (
      (N <= 3) ? NBPoly( [tuple( [rand(1:4, ((N*(N-1))÷2)); 0]... ) for n = 1:ntup],
                        (0.1.+rand(ntup))/factorial(N), DD[N] )
              : NBPoly( [tuple( rand(0:N, ((N*(N-1))÷2+1))... ) for n = 1:ntup],
                      (0.1.+rand(ntup))/factorial(N), DD[N] ) )

   rand_rθ(::BondLengthDesc, N) = SVector( (1.0 .+ rand((N*(N-1))÷2))... )
   rand_rθ(::BondAngleDesc, N) = SVector((1.0.+rand(N-1))...),
                                  SVector( (-0.5.+rand(((N-1)*(N-2))÷2))... )

   println("testing the transform")
   println("---------------------")
   D = D3
   rr = range(r0, stop=r0+1, length=10)
   tdh = (transform.(Ref(D), rr .+ 1e-5) - transform.(Ref(D), rr .- 1e-5)) ./ (2e-5)
   (@test norm(tdh - transform_d.(Ref(D), rr), Inf) < 1e-8) |> println

   println("testing the cutoff")
   println("------------------")
   D = D3
   rr = range(rcut3-0.9, stop=rcut3+0.1, length=100)
   dfh = (fcut.(Ref(D.cutoff), rr .+ 1e-5) - fcut.(Ref(D.cutoff), rr .- 1e-5)) ./ (2e-5)
   (@test norm(dfh - fcut_d.(Ref(D.cutoff), rr), Inf) < 1e-8) |> println

   println(" testing a transformed 4B invariant")
   println("------------------------------------")
   r = SVector( (r0 .+ rand(6))... )
   DI = jac_all_invs(D4, r)
   errs = []
   for p = 2:11
      h = 0.1^p
      DIh = fdjac(D4, r; h = h)
      push!(errs, norm(DIh - DI, Inf))
      @printf(" %.2e | %.2e \n", h, errs[end])
   end
   (@test minimum(errs) <= 1e-3 * maximum(errs)) |> println

   if D isa BondLengthDesc
      println("`NBPoly` gradient-test on simplices")
      println("----------------------------------")
      for N in 2:5, ntup = [1,3]
         VN = random_nbody(N, ntup)
         println("[$N-body, ntup=$ntup, t[1] = $(VN.t[1])]")
         # r = SVector( (r0 + rand((N*(N-1))÷2))... )
         r = rand_rθ(D, N)
         dvN = @D VN(r)
         vN = VN(r)
         rv = Vector(r)
         errs = []
         println("      h   |    err")
         for p = 2:10
            h = 0.1^p
            dvh = zeros(length(dvN))
            for i = 1:length(rv)
               rv[i] += h
               dvh[i] = (VN( SVector(rv...) ) - vN) / h
               rv[i] -=h
            end
            push!(errs, norm(dvN - dvh, Inf))
            @printf(" %.2e | %.2e \n", h, errs[end])
         end
         (@test minimum(errs) <= 1e-3 * maximum(errs)) |> println
      end
   end

   ##
   println("`NBPoly` finite-difference test on configurations")
   println("------------------------------------------------")
   at1 = rattle!(bulk(:Cu, cubic=true) * (1,2,2), 0.02)
   at2 = bulk(:Cu)
   set_constraint!(at2, VariableCell(at2, free = []))
   for at in [at1, at2]
      for N = 3:4, nbas = [1,3]
         println("[$N-body, nbasis = $nbas]")
         VN = random_nbody(N, 1)
         (@test JuLIP.Testing.fdtest(VN, at)) |> println
      end
   end

   ##
   println("Testing Collective Assembly across a Basis Set")
   println("----------------------------------------------")

   BB = [ nothing; [ [random_nbody(N, 1) for _=1:nbasis[N]] for N = 2:5 ] ]
   at = rattle!(bulk(:Cu, cubic=true) * 2, 0.05)

   for N = 2:5
      B = BB[N]
      println("[$N-body]")
      E1 = [energy(b, at) for b in B]
      E2 = energy( B, at )
      print("E: "); (@test E1 ≈ E2) |> println
      F1 = [forces(b, at) for b in B]
      F2 = forces(B, at)
      print("F: "); (@test F1 ≈ F2) |> println
      S1 = [virial(b, at) for b in B]
      S2 = virial(B, at)
      print("S: "); (@test S1 ≈ S2) |> println;
   end

   if profile
      println("Performance Tests")
      println("-----------------")
      println("to get a rough estimate of the cost per atom take the `separate`")
      println("timing and multiply by the factor provided")
      for N = 2:5
         B = BB[N]
         println("[$N-body], length(B) == $(length(B)), factor = $(1/(length(at)*length(B)))")
         print(" E separate:" ); @btime ([energy(b, $at) for b in $B])
         print(" E combined:" ); @btime energy( $B, $at )
         print(" F separate:" ); @btime ([forces(b, $at) for b in $B])
         print(" F combined:" ); @btime forces( $B, $at )
         print(" S separate:" ); @btime ([virial(b, $at) for b in $B])
         print(" S combined:" ); @btime virial( $B, $at )
      end
   end

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE
cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE

cell_to_vtkcell(::Type{Quadrilateral}) = VTKCellTypes.VTK_QUAD
cell_to_vtkcell(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
cell_to_vtkcell(::Type{Triangle}) = VTKCellTypes.VTK_TRIANGLE
cell_to_vtkcell(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE
cell_to_vtkcell(::Type{Cell{2,8,4}}) = VTKCellTypes.VTK_QUADRATIC_QUAD

cell_to_vtkcell(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON
cell_to_vtkcell(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA
cell_to_vtkcell(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA

"""
    vtk_grid(filename::AbstractString, grid::Grid)

Create a unstructured VTK grid from a `Grid`. Return a `DatasetFile`
which data can be appended to, see `vtk_point_data` and `vtk_cell_data`.
"""
function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,N,T}) where {dim,N,T}
    celltype = cell_to_vtkcell(getcelltype(grid))
    cls = MeshCell[]
    for cell in CellIterator(grid)
        push!(cls, MeshCell(celltype, copy(getnodes(cell))))
    end
    coords = reinterpret(T, getnodes(grid), (dim, getnnodes(grid)))
    return vtk_grid(filename, coords, cls)
end

function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, data::Vector{Vec{dim,T}}, name::AbstractString) where {dim,T}
    npoints = length(data)
    data = reinterpret(T, data, (dim, npoints))
    return vtk_point_data(vtk, data, name)
end

function vtk_nodeset(vtk::WriteVTK.DatasetFile, grid::Grid{dim}, nodeset::String) where {dim}
    z = zeros(getnnodes(grid))
    z[collect(getnodeset(grid, nodeset))] = 1.0
    vtk_point_data(vtk, z, nodeset)
end

function vtk_cellset(vtk::WriteVTK.DatasetFile, grid::Grid{dim}, cellset::String) where {dim}
    z = zeros(getncells(grid))
    z[collect(getcellset(grid, cellset))] = 1.0
    vtk_cell_data(vtk, z, cellset)
end


text:
export VaxInt16, VaxInt32

struct VaxInt16 <: VaxInt
    x::UInt16

    VaxInt16(x::UInt16) = new(htol(x))
end

VaxInt16(x::Signed) = VaxInt16(trunc(Int16,x) % UInt16)

Base.convert(::Type{Int16}, x::VaxInt16) = ltoh(x.x) % Int16
function Base.convert(::Type{T}, x::VaxInt16) where T <: Union{Int32,Int64,Int128,BigInt,AbstractFloat}
    return convert(T, convert(Int16, x))
end

struct VaxInt32 <: VaxInt
    x::UInt32

    VaxInt32(x::UInt32) = new(htol(x))
end

VaxInt32(x::Signed) = VaxInt32(trunc(Int32,x) % UInt32)

Base.convert(::Type{Int32}, x::VaxInt32) = ltoh(x.x) % Int32
function Base.convert(::Type{T}, x::VaxInt32) where T <: Union{Int16,Int64,Int128,AbstractFloat}
    return convert(T,convert(Int32,x))
end


text:
"""
ULMFiT - Fine-tuning Language Model

This file contains the funcitons needed to fine-tune a pretrained model.
The novel methods describe the ULMFiT paper are used:

    Discriminative fine-tuning
    Slanted triangular learning rates

"""

"""
    discriminative_step!(layers, ηL::Float64, l, opts::Vector)

Discriminative fine-tuning Step
This function performs the backpropagation step with discriminative fine-tune method,
that is, it uses different learning rates for different layers.

Arguments:

layers      : layers whose weights are going to get updated
ηL          : learning rate of the last layer in `layers`
opts        : `Vector` of optimizers used to update weights for corresponding layers

NOTE: length(opts) == length(layers)
"""
function discriminative_step!(layers, lm::LanguageModel, gen, ηL::Float64, opts::Vector)
    @assert length(opts) == length(layers)
    # Gradient calculation
    grads = Zygote.gradient(() -> loss(lm, gen), get_trainable_params(layers))

    # discriminative step
    ηl = ηL/(2.6^(length(layers)-1))
    for (layer, opt) in zip(layers, opts)
        opt.eta = ηl
        for ps in get_trainable_params([layer])
            Flux.Optimise.update!(opt, ps, grads[ps])
        end
        ηl *= 2.6
    end
    return
end

"""
    fine_tune_lm!(lm::LanguageModel, data_loader::Channel=imdb_fine_tune_data,
            stlr_cut_frac::Float64=0.1, stlr_ratio::Float32=32, stlr_η_max::Float64=4e-3;
            epochs::Integer=1, checkpoint_itvl::Integer=5000)

This function contains main training loops for fine-tuning the language model.
To use this function, an instance of LanguageModel and a data loader is needed.
Read the docs for more info about arguments
"""
function fine_tune_lm!(lm=LanguageModel(), data_loader=imdb_fine_tune_data,
        stlr_cut_frac::Float64=0.1, stlr_ratio::Float32=Float32(32), stlr_η_max::Float64=4e-3;
        epochs::Integer=1, checkpoint_itvl::Integer=5000)

    opts = [ADAM(0.001, (0.7, 0.99)) for i=1:4]

    # Fine-Tuning loops
    for epoch=1:epochs
        println("\nEpoch: $epoch")
        gen = data_loader()
	num_of_iters = take!(gen)
	cut = num_of_iters * epochs * stlr_cut_frac
        T = num_of_iters-Int(floor((num_of_iters*2)/100))
        set_trigger!.(T, lm.layers)
        for i=1:num_of_iters
            # Slanted triangular learning rate step
            t = i + (epoch-1)*num_of_iters
            p_frac = (i < cut) ? i/cut : (1 - ((i-cut)/(cut*(1/stlr_cut_frac-1))))
            ηL = stlr_η_max*((1+p_frac*(stlr_ratio-1))/stlr_ratio)

            # Backprop with discriminative fine-tuning step
            discriminative_step!(lm.layers[[1, 3, 5, 7]], lm, gen, ηL, opts)

            # Resets dropout masks for all the layers with DropOut or DropConnect
            reset_masks!.(lm.layers)
        end
    end
end

"""
    get_vocab(corpus::TokenDocument, lm::LanguageModel)

Returns vocabulary of the given tokenized corpus based on the given language model vocabulary.
"""
function get_vocab(corpus::TokenDocument, lm::LanguageModel)
    u = unique(tokens(corpus))
    new_vocab = intersect(vocab, lm.vocab)
    return new_vocab
end

"""
    set_vocab!(lm::LanguageModel, vocab::Vector)

Sets new vocabulary to the LanguageModel in `vocab` field.
With keeping the same embeddings of previous vocabulary.

Recommended to use this with `get_vocab` function:

Example:
julia> new_vocab = get_vocab(corpus, lm);

julia> set_vocab!(lm, new_vocab)

In addtion to this also incorporate special tokens for UNKNOWN (e.g. "_unk_"), PADDING (e.g. "_pad_") etc.
To add these use `insert!` function after getting vocabulary from `get_vocab` function.
To use `set_vocab!` it is important that the tokens added are same as in the previous vocabulary of Language model,
so that the embeddings for these tokens can be extracted for the previous embeddings of Language model.

# Example:

julia> insert!(vocab, 1, "_unk_")
julia> insert!(vocab, 2, "_pad_")
"""
function set_vocab!(lm::LanguageModel, vocab::Vector)
    idxs = indices(vocab, lm.vocab)
    lm.vocab = vocab
    lm.layers[1].emb = param(lm.layers[1].emb[idxs, :])
    lm.layers[1].mask = gpu(drop_mask((length(vocab),), lm.layers[1].p))
    return
end


text:
module FHCustomSnowballLeaveModule

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/StopCustomSnowballTrigger" CustomSnowballLeaveTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16)

const placements = Ahorn.PlacementDict(
	"Stop Custom Snowball (Frost Helper)" => Ahorn.EntityPlacement(
		CustomSnowballLeaveTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
function gradients!(dynamics::Model, data::ProblemData;
        mode=:nominal)
    x, u, w = trajectories(data,
        mode=mode)
    jx = data.model.jacobian_state
    ju = data.model.jacobian_action
    jacobian!(jx, ju, dynamics, x, u, w)
end

function gradients!(objective::Objective, data::ProblemData;
        mode=:nominal)
    x, u, w = trajectories(data,
        mode=mode)
    gx = data.objective.gradient_state
    gu = data.objective.gradient_action
    gxx = data.objective.hessian_state_state
    guu = data.objective.hessian_action_action
    gux = data.objective.hessian_action_state
    cost_gradient!(gx, gu, objective, x, u, w)
    cost_hessian!(gxx, guu, gux, objective, x, u, w)
end

function gradients!(objective::AugmentedLagrangianCosts, data::ProblemData;
    mode=:nominal)
    # objective
    gx = data.objective.gradient_state
    gu = data.objective.gradient_action
    gxx = data.objective.hessian_state_state
    guu = data.objective.hessian_action_action
    gux = data.objective.hessian_action_state

    # constraints
    constraints = objective.constraint_data.constraints
    c = objective.constraint_data.violations
    cx = objective.constraint_data.jacobian_state
    cu = objective.constraint_data.jacobian_action
    ρ = objective.constraint_penalty
    λ = objective.constraint_dual
    a = objective.active_set
    Iρ = objective.constraint_penalty_matrix
    c_tmp = objective.constraint_tmp
    cx_tmp = objective.constraint_jacobian_state_tmp
    cu_tmp = objective.constraint_jacobian_action_tmp

    # horizon
    H = length(objective)

    # derivatives
    gradients!(objective.costs, data,
        mode=mode)
    gradients!(objective.constraint_data, data,
        mode=mode)

    for t = 1:H
        num_constraint = constraints[t].num_constraint
        for i = 1:num_constraint
            Iρ[t][i, i] = ρ[t][i] * a[t][i]
        end
        c_tmp[t] .= λ[t]

        # gradient_state
        mul!(c_tmp[t], Iρ[t], c[t], 1.0, 1.0)
        mul!(gx[t], transpose(cx[t]), c_tmp[t], 1.0, 1.0)

        # hessian_state_state
        mul!(cx_tmp[t], Iρ[t], cx[t])
        mul!(gxx[t], transpose(cx[t]), cx_tmp[t], 1.0, 1.0)

        t == H && continue

        # gradient_action
        mul!(gu[t], transpose(cu[t]), c_tmp[t], 1.0, 1.0)

        # hessian_action_action
        mul!(cu_tmp[t], Iρ[t], cu[t])
        mul!(guu[t], transpose(cu[t]), cu_tmp[t], 1.0, 1.0)

        # hessian_action_state
        mul!(gux[t], transpose(cu[t]), cx_tmp[t], 1.0, 1.0)
    end
end

function gradients!(constraints_data::ConstraintsData, problem::ProblemData;
    mode=:nominal)
    x, u, w = trajectories(problem,
        mode=mode)
    cx = constraints_data.jacobian_state
    cu = constraints_data.jacobian_action
    jacobian!(cx, cu, constraints_data.constraints, x, u, w)
end

function gradients!(problem::ProblemData;
    mode=:nominal)
    gradients!(problem.model.dynamics, problem,
        mode=mode)
    gradients!(problem.objective.costs, problem,
        mode=mode)
end


text:
using LinearAlgebra, Distances
export Diffusion
export hamming_dist, jaccard_dist, diffusion_dist

"""
Compute the Hamming distance between two adjacency matrices (directed graphs)
"""
function hamming_dist(A1::Matrix{Int}, A2::Matrix{Int})
    return sum(abs(x-y) for (x,y) in zip(A1,A2))
end

(d::Distances.Hamming)(A1::Matrix{Int}, A2::Matrix{Int}) = hamming_dist(A1,A2)

(d::Distances.Hamming)(A1::Matrix{Int}, A2::Nothing) = sum(A1)
(d::Distances.Hamming)(A1::Nothing, A2::Matrix{Int}) = d(A2,A1)

"""
Compute the Jaccard distance between two adjacency matrices (directed graphs)
"""
function jaccard_dist(A1::Matrix{Int}, A2::Matrix{Int})
    return 1 - sum(min(x,y) for (x,y) in zip(A1,A2)) / sum(max(x,y) for (x,y) in zip(A1,A2))
end

(d::Distances.Jaccard)(A1::Matrix{Int}, A2::Matrix{Int}) = jaccard_dist(A1,A2)

struct Diffusion <: Metric end

"""
Compute the diffusion distance between two adjacency matrices with t the time
of diffusion.
"""
function diffusion_dist(A1::Matrix{Int}, A2::Matrix{Int}, t)
    L1 = Diagonal(vec(sum(A1, dims = 2))) - A1
    L2 = Diagonal(vec(sum(A2, dims = 2))) - A2
    function matrixExp(A, t)
        F = eigen(A)
        return F.vectors * Diagonal(exp.(F.values * t)) * transpose(F.vectors)
    end
    return sum((matrixExp(-L1, t) - matrixExp(-L2, t)) .^ 2)
end

(d::Diffusion)(A1::Matrix{Int}, A2::Matrix{Int}) = diffusion_dist(A1,A2)


text:
using OrdinaryDiffEq, DiffEqCallbacks, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear

prob = prob_ode_2Dlinear

integrator = init(prob,Tsit5();dt=1//2^(4))
dt = integrator.dt
solve!(integrator)

u = copy(integrator.sol.u)
t = copy(integrator.sol.t)

reinit!(integrator)
integrator.dt = dt
solve!(integrator)

@test u == integrator.sol.u
@test t == integrator.sol.t

# Implicitly test if dt resets
integrator = init(prob,Tsit5())
dt = integrator.dt
solve!(integrator)

u = copy(integrator.sol.u)
t = copy(integrator.sol.t)

reinit!(integrator)
integrator.dt
solve!(integrator)

@test u == integrator.sol.u
@test t == integrator.sol.t

integrator = init(prob,Tsit5();dt=1//2^(4),tstops = [1//2], saveat = [1//4])
dt = integrator.dt
solve!(integrator)

u = copy(integrator.sol.u)
t = copy(integrator.sol.t)

reinit!(integrator)
integrator.dt = dt
solve!(integrator)

@test u == integrator.sol.u
@test t == integrator.sol.t

#callback test
g(u,p,t) = 2.0*t-2.0
u0=0.0
tspan = (0.0,4.0)
prob = ODEProblem(g,u0,tspan)
saved_values = SavedValues(Float64, Float64)
cb = SavingCallback((u,t,integrator)->u, saved_values)
integrator = init(prob,Tsit5();dt=1//2^(4),callback = cb)
dt = integrator.dt
solve!(integrator)

u = saved_values.saveval
t = saved_values.t
resize!(saved_values.t,0)
resize!(saved_values.saveval,0)
reinit!(integrator)
integrator.dt = dt
solve!(integrator)

@test u == saved_values.saveval
@test t == saved_values.t

@testset "set u0" begin
  prob = prob_ode_2Dlinear
  integrator = init(prob,Tsit5())
  u0 = prob.u0 .+ 1  # just make it different
  @test u0 != prob.u0
  reinit!(integrator, u0)
  @test integrator.u == u0
  @test integrator.sol.u[1] == u0
  @test integrator.sol.interp.timeseries[1] == u0
end

@testset "set u0 with save_idxs" begin
  save_idxs = [1]
  prob = prob_ode_2Dlinear
  integrator = init(prob,Tsit5(); save_idxs=save_idxs)
  u0 = prob.u0 .+ 1  # just make it different
  @test u0 != prob.u0
  reinit!(integrator, u0)
  @test integrator.u == u0
  @test integrator.sol.u[1] == u0[save_idxs]
  @test integrator.sol.interp.timeseries[1] == u0[save_idxs]
end

@testset "set t0" begin
  prob = prob_ode_2Dlinear
  integrator = init(prob,Tsit5())
  t0 = prob.tspan[1] - 1  # just make it different
  @test t0 != prob.tspan[1]
  reinit!(integrator; t0=t0)
  @test integrator.t == t0
  @test integrator.sol.t[1] == t0
  @test integrator.sol.interp.ts[1] == t0
end


text:
using OrthogonalPolynomialsQuasi, FillArrays, BandedMatrices, ContinuumArrays, ArrayLayouts, LazyArrays, Base64, Test
import OrthogonalPolynomialsQuasi: NormalizationConstant, NormalizedBasisLayout, recurrencecoefficients, Normalized, Clenshaw, weighted
import LazyArrays: CachedVector, PaddedLayout
import ContinuumArrays: MappedWeightedBasisLayout

@testset "Normalized" begin
    @testset "Legendre" begin
        P = Legendre()
        Q = Normalized(P)

        @testset "Basic" begin
            @test MemoryLayout(Q) isa NormalizedBasisLayout
            @test @inferred(Q\Q) ≡ Eye(∞)
            @test Q == Q
            @test P ≠ Q
            @test Q ≠ P
            @test Q ≠ P[:,1:end]
            @test P[:,1:end] ≠ Q
        end

        @testset "recurrencecoefficients" begin
            A,B,C = recurrencecoefficients(Q)
            @test A[3:∞][1:10] == A[3:12]
            @test B[3:∞] ≡ Zeros(∞)
        end

        @testset "Evaluation" begin
            M = P'P
            @test Q[0.1,1] == 1/sqrt(2)
            @test Q[0.1,2] ≈ sqrt(1/M[2,2]) * P[0.1,2]
            @test Q[0.1,Base.OneTo(10)] ≈ Q[0.1,1:10] ≈ sqrt.(inv(M)[1:10,1:10]) * P[0.1,Base.OneTo(10)]
            @test (Q'Q)[1:10,1:10] ≈ I
        end

        @testset "Expansion" begin
            f = Q*[1:5; zeros(∞)]
            @test f[0.1] ≈ Q[0.1,1:5]'*(1:5) ≈ f[[0.1]][1]
            x = axes(f,1)
            @test MemoryLayout(Q \ (1 .- x.^2)) isa PaddedLayout
            w = Q * (Q \ (1 .- x.^2));
            @test w[0.1] ≈ (1-0.1^2) ≈ w[[0.1]][1]
        end

        @testset "Conversion" begin
            @test ((P \ Q) * (Q \ P))[1:10,1:10] ≈ (Q \Q)[1:10,1:10] ≈ I
            @test (Jacobi(1,1) \ Q)[1:10,1:10] ≈ ((Jacobi(1,1) \ P) * (P \ Q))[1:10,1:10]
        end

        @testset "Derivatives" begin
            D = Derivative(axes(Q,1))
            f = Q*[1:5; zeros(∞)]
            h = 0.000001
            @test (D*f)[0.1] ≈ (f[0.1+h]-f[0.1])/h atol=1E-4
        end

        @testset "Jacobi" begin
            X = jacobimatrix(Q)
            M = P'P
            @test X[1:10,1:10] ≈ sqrt(M)[1:10,1:10] * jacobimatrix(P)[1:10,1:10] * inv(sqrt(M))[1:10,1:10]
            @test 0.1*Q[0.1,1:10] ≈ (Q*X)[0.1,1:10]
        end

        @testset "Multiplication" begin
            x = axes(Q,1)
            @test Q \ (x .* Q) isa Symmetric

            w = P * (P \ (1 .- x.^2));
            W = Q \ (w .* Q)
            @test W isa Clenshaw
            W̃ = Q' * (w .* Q)
            @test bandwidths(W) == bandwidths(W̃) == (2,2)
            @test W[1:10,1:10] ≈ W[1:10,1:10]' ≈ W̃[1:10,1:10]

            w = @. x + x^2 + 1 # w[x] == x + x^2 + 1
            W = Q \ (w .* Q)
            @test W isa Clenshaw
        end

        @testset "show" begin
            @test stringmime("text/plain", Normalized(Legendre())) == "Normalized(Legendre{Float64})"
            @test summary(Normalized(Legendre()).scaling) == "NormalizationConstant{Float64}"
        end
    end

    @testset "Chebyshev" begin
        T = ChebyshevT()
        w = ChebyshevWeight()
        wT = WeightedChebyshevT()
        Q = Normalized(T)

        @testset "Basic" begin
            @test MemoryLayout(Q) isa NormalizedBasisLayout
            @test @inferred(Q\Q) ≡ Eye(∞)
        end

        @testset "recurrencecoefficients" begin
            A,B,C = recurrencecoefficients(Q)
            @test A[1] ≈ sqrt(2)
            @test A[2:5] ≈ fill(2,4)
            @test C[1:3] ≈ [0,sqrt(2),1]
            @test A[3:∞][1:10] == A[3:12]
            @test B[3:∞] ≡ Zeros(∞)
        end

        @testset "Evaluation" begin
            M = T'wT
            @test Q[0.1,1] == 1/sqrt(π)
            @test Q[0.1,2] ≈ sqrt(1/M[2,2]) * T[0.1,2]
            @test Q[0.1,Base.OneTo(10)] ≈ Q[0.1,1:10] ≈ sqrt.(inv(M)[1:10,1:10]) * T[0.1,Base.OneTo(10)]
            @test (Q'*(w .* Q))[1:10,1:10] ≈ I
        end

        @testset "Expansion" begin
            f = Q*[1:5; zeros(∞)]
            @test f[0.1] ≈ Q[0.1,1:5]'*(1:5) ≈ f[[0.1]][1]
            x = axes(f,1)
            w = Q * (Q \ (1 .- x.^2));
            @test w[0.1] ≈ (1-0.1^2) ≈ w[[0.1]][1]
        end

        @testset "Conversion" begin
            @test ((T \ Q) * (Q \ T))[1:10,1:10] ≈ (Q \Q)[1:10,1:10] ≈ I
            @test (ChebyshevU() \ Q)[1:10,1:10] ≈ ((ChebyshevU() \ T) * (T \ Q))[1:10,1:10]
        end

        @testset "Derivatives" begin
            D = Derivative(axes(Q,1))
            f = Q*[1:5; zeros(∞)]
            h = 0.000001
            @test (D*f)[0.1] ≈ (f[0.1+h]-f[0.1])/h atol=1E-4
        end

        @testset "Multiplication" begin
            x = axes(Q,1)
            @test Q \ (x .* Q) isa Symmetric

            w = T * (T \ (1 .- x.^2));
            W = Q \ (w .* Q)
            @test W isa Clenshaw
            @test bandwidths(W) == (2,2)
            W̃ = Q\  (w .* Q)
            @test W[1:10,1:10] ≈ W[1:10,1:10]' ≈ W̃[1:10,1:10]
        end
    end

    @testset "Jacobi" begin
        Q = Normalized(Jacobi(1/2,0))
        # Emperical from Mathematica
        @test Q[0.1,1:4] ≈ [0.728237657560985,0.41715052371131806,-0.6523500049588019,-0.5607891513201705]
        w = JacobiWeight(1/2,0)
        @test (Q'*(w .* Q))[1:10,1:10] ≈ I
    end

    @testset "Mapped" begin
        P = legendre(0..1)
        x = axes(P,1)
        Q = Normalized(P)

        # Emperical from Mathematica
        @test Q[0.1,1:4] ≈ [1,-1.3856406460551018,1.028591269649903,-0.21166010488516684]

        u = Q[:,1:20] * (Q[:,1:20] \ exp.(x))
        @test u[0.1] ≈ exp(0.1)
        u = Q * (Q \ exp.(x))
        @test u[0.1] ≈ exp(0.1)

        Q = Normalized(jacobi(1/2,0,0..1))
        @testset "Recurrences" begin
            A,B,C = recurrencecoefficients(Q)
            Ã,B̃,C̃ = recurrencecoefficients(Normalized(Jacobi(1/2,0)))
            @test A[1:10] ≈ 2Ã[1:10]
            @test B[1:10] ≈ B̃[1:10] .- Ã[1:10]
            @test C[1:10] ≈ C̃[1:10]
        end
        wQ = weighted(Q)
        x = axes(Q,1)
        @test wQ[0.1,1:10] ≈ Q[0.1,1:10] * sqrt(1-(2*0.1-1))

        u = wQ[:,1:20] * (wQ[:,1:20] \  @.(sqrt(1-x^2)))
        @test u[0.1] ≈ sqrt(1-0.1^2)
        u = wQ * (wQ \ @.(sqrt(1-x^2)))
        @test u[0.1] ≈ sqrt(1-0.1^2)
    end

    @testset "Christoffel–Darboux" begin
        Q = Normalized(Legendre())
        X = Q\ (axes(Q,1) .* Q)
        x,y = 0.1,0.2
        n = 10
        Pn = Diagonal([Ones(n); Zeros(∞)])
        @test (X*Pn - Pn*X)[1:n,1:n] ≈ zeros(n,n)
        @test Pn * Q[y,:] isa CachedVector

        # @test (x-y) * Q[x,1:n]'*Q[y,1:n] ≈ (x-y) * Q[x,:]'*Pn*Q[y,:] ≈ (x-y) * Q[x,:]'*Pn*Q[y,:]
        # Q[x,:]' * ((X*Pn - Pn*X)* Q[y,:])
        @test (x-y) * Q[x,1:n]'*Q[y,1:n] ≈ Q[x,n:n+1]' * (X*Pn - Pn*X)[n:n+1,n:n+1] * Q[y,n:n+1]
    end
end