text:
@testset "Genus" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a", cached = false)
  OK = maximal_order(K)
  rlp = real_places(K)
  sig = Dict(rlp[1] => 2)

  p2 = prime_decomposition(OK, 2)[1][1]
  p3 = prime_decomposition(OK, 3)[1][1]
  p5 = prime_decomposition(OK, 5)[1][1]

  @test length(Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 0)) == 8
  @test length(Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 1)) == 8
  @test length(Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 2)) == 16
  @test length(Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 3)) == 24
  @test length(Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 4)) == 32

  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 0)) == 2
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 1)) == 4
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 2)) == 8
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 3)) == 16
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 4)) == 28
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 5)) == 46
  @test length(Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 6)) == 72

  @test length(Hecke.local_genera_quadratic(K, p5, rank = 4, det_val = 2)) == 8
  @test length(Hecke.local_genera_quadratic(K, p5, rank = 4, det_val = 1)) == 4
  @test length(Hecke.local_genera_quadratic(K, p5, rank = 2, det_val = 0)) == 2
  @test length(Hecke.local_genera_quadratic(K, p5, rank = 2, det_val = 1)) == 4

  @test length(Hecke.genera_quadratic(K, rank = 2, signatures = sig, det = 5^2 * 8 * 9 * OK)) == 27
  @test length(Hecke.genera_quadratic(K, rank = 2, signatures = sig, det = 5 * 7 * 8 * 9 * OK)) == 36
  @test length(Hecke.genera_quadratic(K, rank = 2, signatures = sig, det = 2^5 * OK)) == 5
  @test length(Hecke.genera_quadratic(K, rank = 2, signatures = sig, det = 11 * 13^2 * OK)) == 6

  G = Hecke.local_genera_quadratic(K, p2, rank = 2, det_val = 4)
  for i in 1:10
    G1 = rand(G)
    G2 = rand(G)
    L1 = representative(G1)
    L2 = representative(G2)
    G3 = @inferred orthogonal_sum(G1, G2)
    L3, = orthogonal_sum(L1, L2)
    @test G3 == genus(L3, p2)
  end

  G = Hecke.local_genera_quadratic(K, p3, rank = 5, det_val = 5)
  for i in 1:10
    G1 = rand(G)
    G2 = rand(G)
    L1 = representative(G1)
    L2 = representative(G2)
    G3 = @inferred orthogonal_sum(G1, G2)
    L3, = orthogonal_sum(L1, L2)
    @test G3 == genus(L3, p3)
  end

  G = Hecke.local_genera_quadratic(K, p5, rank = 2, det_val = 1)
  for i in 1:10
    G1 = rand(G)
    G2 = rand(G)
    L1 = representative(G1)
    L2 = representative(G2)
    G3 = @inferred orthogonal_sum(G1, G2)
    L3, = orthogonal_sum(L1, L2)
    @test G3 == genus(L3, p5)
  end

  G = Hecke.genera_quadratic(K, rank = 3, signatures = sig, det = 2 * 9 * OK)
  for i in 1:10
    G1 = rand(G)
    G2 = rand(G)
    L1 = representative(G1)
    L2 = representative(G2)
    G3 = @inferred orthogonal_sum(G1, G2)
    L3, = orthogonal_sum(L1, L2)
    @test L3 in G3
  end

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [30, -15, 0, -15, 30, -15, 0, -15, 30]);
  gens = [[1, 0, 0], [1, 0, 0], [0, 1, 0], [0, 1, 0], [5//6, 2//3, 1//6], [5//6, 2//3, 1//6]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  p = prime_decomposition(maximal_order(K), 5)[1][1]
  fl, LL = Hecke.ismaximal_integral(L, p)
  @test !fl
end

ext: jl


text:
import Arpack: eigs

abstract type AbstractEigenSolver end

struct EigArpack <: AbstractEigenSolver end

function eigsolve(A, ::EigArpack ; args...)
    eigs(A, which = :LR ; args...)
end
eigsolve(A; args...) = eigsolve(A,EigArpack() ; args...)


text:
__precompile__()

module SPTK

# A thin Julia wrapper for Speech signal Processing Toolkit (SPTK).
#
# The original SPTK: http://sp-tk.sourceforge.net/
# Note that the SPTK.jl is based on the slightly modified SPTK
# https://github.com/r9y9/SPTK
#
# NOTE1: no function is exported; i.e., you should explicitly import functions
# if you need.
#
# NOTE2: Function interfaces might be different betweeen C and Julia,
# as the following convensions.
#
# Conventions for wrapping C interface
#
# 1. Avoid really short names for variables (e.g. a, b, c, aa, bb, dd) .
#
#    Variable names should be informative. If the C functions have such short
#    names, use self-descriptive names instead for Julia interfaces, unless
#    they have clear meanings in their context.
#
# 2. Avoid too many function arguments.
#
#    Less is better. If the C functions have too many function arguments, use
#    keyword arguments with proper default values for optional ones in Julia.
#
# 3. Handle errors in Julia
#
#    Since C functions might `exit` (unfortunately) inside thier functions for
#    unexpected inputs, it should be check if the inputs are supported or not
#    in Julia before `ccall`.

deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("SPTK not properly installed. Please run Pkg.build(\"SPTK\")")
end

for fname in [
              "common",
              "lib",
              "adaptive",
              "conversions",
              "excite",
              "f0",
              "mfcc",
              "mgcep",
              "synthesis_filters",
              "utils",
              "window",
              "extend"
    ]
    include(string(fname, ".jl"))
end

end # module SPTK


text:
function setup_model()
    T = Float32

    x = Var()
    y = pairwise(x)
    #y = wordembeds(x)
    y = convolution(y, (50,), strides=(10,), pads=(20,))
    y = pooling(:max, y)
    charfun = compile(y)
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
import Distributed: @distributed
import LinearAlgebra: norm

"""Find the Lorenz96 orbit length that is reached starting from X0.
Recursive function that starts testing for orbits over period length l,
but enlarges that testing period if no orbit was found."""
function orbit_length(  X0::Vector{T};          # initial condition
                        lspinup::Int=100_000,   # spinup
                        lchunk::Int=100_000,    # chunk size
                        verbose::Bool=true,     # feedback?
                        nchunkmax::Int=10_000,  # abort after this many chunks
                        ) where T

    tic = time()
    X0 = L96(T,X=X0,n=length(X0),N=lspinup)
    toc = time()
    if verbose
        time_elapsed = readable_secs(toc-tic)
        println("Spin-up completed. $lspinup time steps in $time_elapsed.")
    end

    # use several X for testing
    #   1: initial conditions, static
    # dynamic X that move forward in time (in case 1 isn't part of the orbit yet)
    #   2: sqrt(nchunk)
    #   3: log1.5(nchunk)
    #   4: cbrt(nchunk)
    #   5: log2(nchunk)
    log1p5(x) = log(x)/log(1.5)
    tests = [sqrt,log1p5,cbrt,log2]
    ntests = length(tests)+1

    Xtest = Array{T,2}(undef,length(X0),ntests)
    for i in 1:ntests
        Xtest[:,i] = X0     # initialise all with X0
    end

    # respective time stamps
    ltest = zeros(Int,ntests)

    l = 0               # orbit length (0 = not found yet)
    nchunk = 0          # count the chunks

    while l == 0 && nchunk < nchunkmax
        # compute next chunk of L96 trajectories
        X = L96(T,X=X0,n=length(X0),N=lchunk,output=true)
        X0 = X[:,end]   # new initial condition for next chunk

        j = 2   # index within chunk, skip initial condition
        @inbounds while l == 0 && j < lchunk  # stop when orbit is found (l>0)
            # check for periodicity
            Xj = X[:,j]
            for i in 1:ntests
                l = Xtest[:,i] == Xj ? nchunk*lchunk+j-ltest[i]-1 : 0
            end
            j += 1
        end

        nchunk += 1

        # move Xtest[2:end] slowly forwards in time
        for (i,func) in enumerate(tests)
            if func(nchunk) % 1 == 0.0
                Xtest[:,i+1] = copy(X0)
                ltest[i+1] = nchunk*lchunk
            end
        end
    end

    if verbose
        if l == 0
            print("No orbit found after $nchunk chunks of size $lchunk")
        else
            println("Orbit of length=$l found after $nchunk chunks of size $lchunk.")
        end
    end

    return l,X0    # return orbit length and an X that's on the orbit
end

"""Returns the minimum X on the orbit of length l, starting from X0 on the orbit."""
function orbit_minimum( X0::Vector{T},               # initial condition on the orbit
                        l::Int;                      # the period length of the orbit
                        lchunk::Int=100_000) where T     # chunk size to reduce memory stress

    Xmin = copy(X0)         # start with X0 as minimum
    Xnorm = norm(Xmin)

    ll = 0                  # counter to count up till l for each chunk
    while ll < l
        X = L96(T,X=X0,n=length(X0),N=lchunk,output=true)
        X0 = X[:,end]       # last time step are the new initial conditions

        # find the L2-norm minimum by running through X
        @inbounds for i in 2:lchunk
            normXi = norm(X[:,i])
            if normXi < Xnorm
                Xmin = X[:,i]
                Xnorm = normXi
            end
        end

        ll += lchunk
    end

    return Xmin
end

"""Returns both the length and minimum of the orbit reached from initial condition X0."""
function orbit_length_minimum(  X0::Vector{T};              # initial condition on the orbit
                                lspinup::Int=100_000,       # spin-up length
                                lchunk::Int=100_000,
                                nchunkmax::Int=10_000,
                                verbose::Bool=true) where T    # chunk size
    len,X = orbit_length(X0;lspinup,lchunk,nchunkmax,verbose)
    Xmin = orbit_minimum(X,len;lchunk)
    return Orbit(len,Xmin)
end

"""For a given number format T, number of variables N and n initial conditions, find orbits of
Lorenz96 and test for their uniqueness."""
function find_orbits(   ::Type{T},                  # Number format
                        N::Int,                     # number of variables in L96
                        n::Int=100;                 # number of initial conditions
                        lini::Int=100_000,          # initial spin-up length in Float64
                        lspinup::Int=100_000,       # spinup with format T
                        lchunk::Int=100_000,        # chunk size to reduce memory stress
                        nchunkmax::Int=10_000,      # maximum number of chunks to abort
                        verbose::Bool=true          # report every orbit found?
                        ) where T         # n initial conditions

    # pre-allocate empty array of orbits
    orbits = Orbit[]

    # create initial conditions, discard spinup
    lini = max(lini,10*n)       # increase lini in case of many initial conditions
    spinup = 100
    ini = L96(Float64,n=N,N=lini+spinup,output=true,η=0.005+0.01*rand())
    ini = ini[:,unique(rand(spinup:lini+spinup,10*n))[1:n]]

    tic = time()
    orbits = @distributed (reduce_orbits) for i in 1:n            # for n ICs calculate orbit lengths & x
        X = convert(Vector{T},ini[:,i])
        orbit = orbit_length_minimum(X;lspinup,lchunk,nchunkmax,verbose)
        if verbose
            println(orbit)
        end
        orbit
    end

    # in case n=1 pack into a Vector
    orbits = isa(orbits,Orbit) ? [orbits] : orbits

    sort!(orbits)
    normalise_basins!(orbits)

    toc = time()
    time_elapsed = readable_secs(toc-tic)

    norbits = sum([orbit.length > 0 for orbit in orbits])
    println("Found $norbits orbits in $time_elapsed.")
    return orbits
end


text:
struct TrivialMeasure{X} <: AbstractMeasure{X} end

TrivialMeasure(X) = TrivialMeasure{X}()


text:
using Test, SolveDiffEq

print_head(s::String) = println("\033[1m\033[32mTesting Task\033[0m \033[33m$s\033[0m")
print_body(s::String) = println("             $s")

@testset "Test Bulirsch-Stoer Method" begin
    print_head("y'' - 6y' + 15y = 2sin(3t), y(0) = -1, y'(0) = -4")

    function test_func(t::Real)
        temp = sqrt(6.0) * t
        arg1 = cos(3.0t) + sin(3.0t) / 3.0
        arg2 = 11.0 * exp(3.0t) * cos(temp)
        arg3 = 8.0 * exp(3.0t) * sin(temp) / sqrt(6.0)
        return 0.1 * (arg1 - arg2 - arg3)
    end

    function test_diff!(du::AbstractVector, un::AbstractVector, tn::Real)
        @inbounds du[1] = un[2]
        @inbounds du[2] = -15.0 * un[1] + 6.0 * un[2] + 2.0 * sin(3.0tn)
        return nothing
    end

    Δt = 0.05; ans = test_func(1.0)

    obj = BulirschStoer{7}(Δt, 2)                # stepper
    arr = collect(range(0.0, 1.0; step=Δt))      # t vector
    mat = Matrix{Float64}(undef, 2, length(arr)) # y matrix

    # initial values
    @inbounds mat[1, 1], mat[2, 1] = -1.0, -4.0
    # iterations
    @inbounds for i in 1:length(arr)-1
        apply_step!(view(mat,:,i+1), view(mat,:,i), arr[i], test_diff!, obj)
    end

    print_body("answer = $ans")
    print_body("approx = $(mat[1, end])")
    @test mat[1, end] ≈ ans
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random
using ArraysOfArrays, IntervalSets, ValueShapes

struct apb_test <: BAT.AbstractVarBounds
    varndof::Integer
end

ValueShapes.totalndof(a::apb_test) = a.varndof
BAT.unsafe_intersect(a::apb_test, b::apb_test) = true

@testset "parameter bounds" begin
    @testset "BAT.renormalize_variate_impl" begin
        @test BAT.renormalize_variate_impl(+0.3, -1, 2, BAT.hard_bounds) ≈ +0.3
        @test BAT.renormalize_variate_impl(-0.3, -1, 2, BAT.reflective_bounds) ≈ -0.3
        @test BAT.renormalize_variate_impl(+1.3, -1, 2, BAT.cyclic_bounds) ≈ +1.3

        @test BAT.renormalize_variate_impl(-1.3, -1, 2, BAT.hard_bounds) == -1.3
        @test BAT.renormalize_variate_impl(2.3, -1, 2, BAT.hard_bounds) == 2.3

        @test BAT.renormalize_variate_impl(-1.3, -1, 2, BAT.reflective_bounds) ≈ -0.7
        @test BAT.renormalize_variate_impl(-4.3, -1, 2, BAT.reflective_bounds) ≈ +1.7
        @test BAT.renormalize_variate_impl(-7.3, -1, 2, BAT.reflective_bounds) ≈ -0.7
        @test BAT.renormalize_variate_impl(+2.3, -1, 2, BAT.reflective_bounds) ≈ +1.7
        @test BAT.renormalize_variate_impl(+5.3, -1, 2, BAT.reflective_bounds) ≈ -0.7
        @test BAT.renormalize_variate_impl(+8.3, -1, 2, BAT.reflective_bounds) ≈ +1.7

        @test BAT.renormalize_variate_impl(-1.3, -1, 2, BAT.cyclic_bounds) ≈ +1.7
        @test BAT.renormalize_variate_impl(-4.3, -1, 2, BAT.cyclic_bounds) ≈ +1.7
        @test BAT.renormalize_variate_impl(-7.3, -1, 2, BAT.cyclic_bounds) ≈ +1.7
        @test BAT.renormalize_variate_impl(+2.3, -1, 2, BAT.cyclic_bounds) ≈ -0.7
        @test BAT.renormalize_variate_impl(+5.3, -1, 2, BAT.cyclic_bounds) ≈ -0.7
        @test BAT.renormalize_variate_impl(+8.3, -1, 2, BAT.cyclic_bounds) ≈ -0.7

        @test typeof(@inferred BAT.renormalize_variate_impl(+0.3, -1, 2, BAT.hard_bounds)) == Float64
        @test typeof(@inferred BAT.renormalize_variate_impl(+0.3f0, -1, 2, BAT.hard_bounds)) == Float32
        @test typeof(@inferred BAT.renormalize_variate_impl(+0.3f0, -1, 2.0, BAT.hard_bounds)) == Float64

        @test BAT.renormalize_variate_impl(+5.3, ClosedInterval(-1, 2), BAT.reflective_bounds) ≈ -0.7
    end

    @testset "BAT.NoVarBounds" begin
        n = 2
        @test typeof(@inferred BAT.NoVarBounds(n)) == BAT.NoVarBounds

        v = [-1000., 1000]
        nobounds = BAT.NoVarBounds(n)
        @test @inferred(totalndof(nobounds)) == n
        @test @inferred(v in nobounds)
        @test @inferred(BAT.renormalize_variate!(v, nobounds, copy(v))) === v
    end

    @testset "BAT.HyperRectBounds" begin
        @test typeof(@inferred BAT.HyperRectBounds([-1., 0.5], [2.,1], BAT.BAT.hard_bounds)) <: BAT.VarVolumeBounds{Float64, BAT.HyperRectVolume{Float64}}
        @test_throws ArgumentError BAT.HyperRectBounds([-1.], [2.,1], [BAT.hard_bounds, BAT.reflective_bounds])
        @test_throws ArgumentError BAT.HyperRectBounds([-1., 2.], [2.,1], [BAT.hard_bounds, BAT.reflective_bounds])

        hyperRectBounds = @inferred BAT.HyperRectBounds([-1., -1.], [0.5,1], [BAT.hard_bounds, BAT.hard_bounds])
        @test totalndof(hyperRectBounds) == 2
        @test [0.0, 0.0] in hyperRectBounds
        @test ([0.5, 2] in hyperRectBounds) == false

        hyperRectBounds = @inferred BAT.HyperRectBounds([ClosedInterval(-1.,0.5),
            ClosedInterval(-1.,1.)], [BAT.hard_bounds, BAT.hard_bounds])
        @test totalndof(hyperRectBounds) == 2
        @test [0.0, 0.0] in hyperRectBounds
        @test ([0.5, 2] in hyperRectBounds) == false

        @test @inferred(BAT.renormalize_variate!(
            fill(NaN, 3),
            BAT.HyperRectBounds([-1.,-1,-1], [2.,2,2], [BAT.hard_bounds, BAT.reflective_bounds, BAT.cyclic_bounds]),
            [0.3, -4.3, -7.3]
        )) ≈ [+0.3, 1.7, 1.7]

        @test BAT.renormalize_variate!(fill(NaN, 2), hyperRectBounds, [0, 2.]) == [0,2.]
    end

    @testset "similar" begin
        a = @inferred BAT.HyperRectBounds([-1., 0.5], [2.,1], BAT.BAT.hard_bounds)
        c = @inferred similar(a)
        @test typeof(c) <:BAT.HyperRectBounds{Float64}
        @test typeof(c.vol) <: BAT.HyperRectVolume{Float64}
        @test typeof(c.bt) <: Array{BAT.BoundsType,1}
    end

    @testset "BAT.intersect" begin
        a = @inferred apb_test(3)
        b = @inferred apb_test(4)
        @test_throws ArgumentError intersect(a,b)
        b = @inferred apb_test(3)
        @test intersect(a,b)
        b = @inferred BAT.NoVarBounds(3)
        @test BAT.unsafe_intersect(b, BAT.NoVarBounds(3)) == b
        @test BAT.unsafe_intersect(a, BAT.NoVarBounds(3)) == a
        @test BAT.unsafe_intersect(BAT.NoVarBounds(3), a) == a

        hb = BAT.hard_bounds
        rb = BAT.reflective_bounds
        cb = BAT.cyclic_bounds

        @test intersect(hb, rb) == BAT.hard_bounds
        @test intersect(cb, hb) == BAT.hard_bounds
        @test intersect(rb, cb) == BAT.reflective_bounds
        @test intersect(cb, rb) == BAT.reflective_bounds
        @test intersect(cb, cb) == BAT.cyclic_bounds

        hyperRectBounds_a = @inferred BAT.HyperRectBounds([-1., -1.], [2.,3.],
            [BAT.hard_bounds, BAT.cyclic_bounds])
        hyperRectBounds_b = @inferred BAT.HyperRectBounds([-0.5, -0.5], [1.,2.],
            [BAT.cyclic_bounds, BAT.cyclic_bounds])

        res = @inferred intersect(hyperRectBounds_a, hyperRectBounds_b)
        @test res.vol.lo ≈ hyperRectBounds_b.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_b.vol.hi
        @test res.bt == hyperRectBounds_b.bt

        res = @inferred intersect(hyperRectBounds_b, hyperRectBounds_a)
        @test res.vol.lo ≈ hyperRectBounds_b.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_b.vol.hi
        @test res.bt == hyperRectBounds_b.bt

        hyperRectBounds_c = @inferred BAT.HyperRectBounds([-0.5, -0.5], [4.,4.],
            [BAT.cyclic_bounds, BAT.reflective_bounds])

        res = @inferred intersect(hyperRectBounds_a, hyperRectBounds_c)
        @test res.vol.lo ≈ hyperRectBounds_c.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_a.vol.hi
        @test res.bt == [BAT.hard_bounds, BAT.reflective_bounds]

    end
end


text:
# Include packages
using FFTW: fft, ifft

# MAT package is not mature.
# When it returns a error, comment the following line and
# read data from the txt formatted file.
using MAT: MAT

##############
## functions
##############
function fftfreq(nt::Int,dt)
    # 演習で使用するデータのサンプル数が偶数個なので調整
    # 奇数の場合は正常に動かないかも．未確認
    if isodd(nt); nt2 = Int((nt-1)/2); else; nt2 = Int(nt/2); end
    T =nt*dt;
    freq = zeros(nt);
    for it = 2:nt2
        freq[it] = (it-1)/T;
    end
    for it = nt2+1:nt
        freq[it] = (it-(nt+2))/T;
    end
    return freq
end
##############
function loadmat()
    # define the filepath & filename
    fdir = "./data";
    fname = "crf_wind.mat";
    # file open & get variables
    matfile = MAT.matopen(join([fdir,fname],"/"));
    dt, nt, nz = read(matfile,"dt"), Int(read(matfile,"nt")), Int(read(matfile,"nz"));
    z, u = read(matfile,"z"), read(matfile,"u");
    return dt, nt, nz, z, u
end
##############
function loadtxt()
    # define the filepath & filename
    fdir = "./data"
    fname = "crf_wind.dat"
    # file open & get variables
    f = open(join([fdir,fname],"/"),"r")
    txtorg = readlines(f)
    close(f)
    head = split(txtorg[1],r"\s+")
    nt = parse(Int64, head[1])
    nz = parse(Int64, head[2])
    dt = parse(Float64, txtorg[2])
    z = [parse(Float64, txtorg[3][12(i-1)+1:12i]) for i=1:nz]
    u = [parse(Float64, txtorg[i+3][12(j-1)+1:12j]) for i=1:nz, j=1:nt]
    # return varout
   return dt, nt, nz, z, u
end
##############

####################
## main
####################

# Data source mat or txt
dt, nt, nz, z, u = loadmat() # mat
#dt, nt, nz, z, u = loadtxt() # txt

# Fourier transform
data = u[1,:];
F0 = fft(data,1);
P = abs.(F0/(nt/2)); # Power
freq = fftfreq(nt,dt);

# noise reduction
fc = 1/100dt # cut off　※この値に根拠はありません．
cutoff = abs.(freq) .> fc;
freq0 = iszero.(freq);
F0[cutoff .& .!freq0] .= 0.0
datamod = ifft(F0);

####################

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()
#plotlyjs()

########
## Figure 1
#  Power spectrum density
if isodd(nt); nt2 = Int((nt-1)/2); else; nt2 = Int(nt/2); end
plot(freq[2:nt2], P[2:nt2], lab="Power spectrum", legend=false,
     xscale=:log10, yscale=:log10,
     xlabel="Frequency (Hz)", ylabel="Power (m²/s²)",
     xlims=(2e-4,1e+1), ylims=(1e-4, 2e+0),
     )
savefig(joinpath(figdir,"PSD2.png"))
########

########
## Figure 2
# original and noise reduced data
t = 0:dt:(nt-1)*dt;
plot(t, data, line=(:solid, 1), lab="Raw data", size=(1500, 600),
     xlabel="Time (s)", ylabel="Wind Speed (m/s)", guidefont=14,
     tickfont=12,
     )
plot!(t, real.(datamod), line=(:solid, 2), color=:magenta, lab="Noise reduced",
      legend=:topright, legendfont=14,
      )
savefig(joinpath(figdir,"noise_reduced.png"))
########

####################


text:
function execute(problem::FSIProblem{:oscillator}; kwargs...)

  # Define cylinder motion
  u_y = _get_kwarg(:u_y,kwargs)
  du_y = _get_kwarg(:du_y,kwargs)
  @assert(typeof(u_y)<:Function)
  @assert(typeof(du_y)<:Function)

  # Problem setting (Default FSI-2)
  println("Setting Forced Oscillator cylinder problem parameters")
  Um = _get_kwarg(:Um,kwargs,1.0)
  ⌀ = _get_kwarg(:D,kwargs,0.1)
  # Fluid properties
  ρ_f = _get_kwarg(:rho_f,kwargs,1.0e3)
  Re = _get_kwarg(:Re,kwargs, 100.0)
  μ_f = ρ_f * Um * ⌀ / Re
  γ_f = _get_kwarg(:gamma_f,kwargs,1.0)
  # Mesh properties
  E_m = _get_kwarg(:E_m,kwargs,1.0)
  ν_m = _get_kwarg(:nu_m,kwargs,-0.1)
  α_m = _get_kwarg(:alpha_m,kwargs,1.0e-5)
  weight_strategy = _get_kwarg(:alpha_m_weight,kwargs,"constant")
  # Time stepping
  t0 = _get_kwarg(:t0,kwargs,0.0)
  tf = _get_kwarg(:tf,kwargs,0.1)
  dt = _get_kwarg(:dt,kwargs,0.1)
  θ  = _get_kwarg(:theta,kwargs,0.5)
  # Post-process
  is_vtk = _get_kwarg(:is_vtk,kwargs,false)

  # Mesh strategy
  strategyName = _get_kwarg(:strategy,kwargs,"biharmonic")
  strategy = WeakForms.MeshStrategy{Symbol(strategyName)}()

  # Coupling (boundary conditions) type
  couplingName = _get_kwarg(:coupling,kwargs,"strong")
  coupling = WeakForms.Coupling{Symbol(couplingName)}()

  # Define BC functions
  println("Defining Boundary conditions")
  u_in(x, t) = VectorValue( Um , 0.0)
  u_sym(x, t) = 0.0
  u_cylinder(x,t) = VectorValue(0.0,u_y(t))
  u_in(t::Real) = x -> u_in(x, t)
  u_sym(t::Real) = x -> u_sym(x, t)
  u_cylinder(t::Real) = x -> u_cylinder(x, t)
  ∂tu_in(t) = x -> VectorValue(0.0, 0.0)
  ∂tu_sym(t) = x -> 0.0
  ∂tu_cylinder(t) = x -> VectorValue(0.0, du_y(t))
  ∂tu_in(x, t) = ∂tu_in(t)(x)
  ∂tu_sym(x, t) = ∂tu_sym(t)(x)
  ∂tu_cylinder(x, t) = ∂tu_cylinder(t)(x)
  T_in = typeof(u_in)
  T_sym = typeof(u_sym)
  T_cylinder = typeof(u_cylinder)
  @eval ∂t(::$T_in) = $∂tu_in
  @eval ∂t(::$T_sym) = $∂tu_sym
  @eval ∂t(::$T_cylinder) = $∂tu_cylinder
  bconds = get_boundary_conditions(problem,strategy,coupling,u_in,u_sym,u_cylinder,∂tu_cylinder)

  # Forcing terms
  f(t) = x -> VectorValue(0.0,0.0)

  # Discrete model
  println("Defining discrete model")
  modelName = _get_kwarg(:model,kwargs,"../models/cylinder_NSI.json")
  model = DiscreteModelFromFile(modelName)
  writevtk(model,"cylinder")

  # Triangulations
  println("Defining triangulations")
  trian = Triangulation(model)
  trian_Γc = BoundaryTriangulation(model,"cylinder")
  n_Γc = get_normal_vector(trian_Γc)

  # Compute cell area (auxiliar quantity for mesh motion eq.)
  if( weight_strategy == "volume")
    vol = cell_measure(trian,trian)
    α_Ω = α_m * vol
    α_Γc = α_m * reindex(vol,trian_Γc)
  else
    α_Ω = α_m; α_Γc = α_m
  end

  # Compute interface element size (if weak coupling)
  if ( typeof(coupling) == WeakForms.Coupling{:weak} )
    hΓ = reindex(cell_measure(trian_Γc,trian),trian_Γc)
  else
    hΓ = 0.0
  end

  # Quadratures
  println("Defining quadratures")
  order = _get_kwarg(:order,kwargs,2)
  degree = 2*order
  bdegree = 2*order
  quad = CellQuadrature(trian,degree)
  quad_Γc = CellQuadrature(trian_Γc,bdegree)

  # Test FE Spaces
  println("Defining FE spaces")
  Y_ST, X_ST, Y_NSI, X_NSI = get_FE_spaces(problem,strategy,coupling,model,order,bconds)

  # Stokes problem for initial solution
  println("Defining Stokes operator")
  res_ST(x,y) = WeakForms.stokes_residual(strategy,x,y,μ_f,f(0.0))
  jac_ST(x,dx,y) = WeakForms.stokes_jacobian(strategy,dx,y,μ_f)
  t_ST_Ω = FETerm(res_ST, jac_ST, trian, quad)
  op_ST = FEOperator(X_ST,Y_ST,t_ST_Ω)

  # Setup equation parameters
  nsi_f_params = Dict(
  :μ=>μ_f,
  :ρ=>ρ_f,
  :E=>E_m,
  :ν=>ν_m,
  :α=>α_Ω,
  :fu=>f,
  :fv=>f,
  )
  nsi_Γc_params = Dict(
  :n=>n_Γc,
  :E=>E_m,
  :ν=>ν_m,
  :μ=>μ_f,
  :α=>α_Γc,
  :γ=>γ_f,
  :h=>hΓ,
  :vD=>du_y
  )

   # FSI problem
  println("Defining FSI operator")
  res_NSI_Ω(t,x,xt,y) = WeakForms.fluid_residual_Ω(strategy,t,x,xt,y,nsi_f_params)
  jac_NSI_Ω(t,x,xt,dx,y) = WeakForms.fluid_jacobian_Ω(strategy,x,xt,dx,y,nsi_f_params)
  jac_t_NSI_Ω(t,x,xt,dxt,y) = WeakForms.fluid_jacobian_t_Ω(strategy,x,xt,dxt,y,nsi_f_params)
  t_NSI_Ω = FETerm(res_NSI_Ω, jac_NSI_Ω, jac_t_NSI_Ω, trian, quad)
  if(typeof(coupling)==WeakForms.Coupling{:weak})
    res_NSI_Γc(t,x,xt,y) = WeakForms.fluid_residual_Γ(strategy,t,x,xt,y,nsi_Γc_params)
    jac_NSI_Γc(t,x,xt,dx,y) = WeakForms.fluid_jacobian_Γ(strategy,x,xt,dx,y,nsi_Γc_params)
    jac_t_NSI_Γc(t,x,xt,dxt,y) = WeakForms.fluid_jacobian_t_Γ(strategy,dxt,y,nsi_Γc_params)
    t_NSI_Γc = FETerm(res_NSI_Γc,jac_NSI_Γc,jac_t_NSI_Γc,trian_Γc,quad_Γc)
    op_FSI = TransientFEOperator(X_NSI,Y_NSI,t_NSI_Ω,t_NSI_Γc)
  else
    op_FSI = TransientFEOperator(X_NSI,Y_NSI,t_NSI_Ω)
  end

  # Setup output files
  folderName = "oscillator-results"
  fileName = "fields"
  if !isdir(folderName)
    mkdir(folderName)
  end
  filePath = join([folderName, fileName], "/")

  # Solve Stokes problem
  @timeit "ST problem" begin
    println("Solving Stokes problem")
    xh = solve(op_ST)
    if(is_vtk)
      writePVD(filePath, trian, [(xh, 0.0)])
    end
  end

  # Solve FSI problem
  @timeit "FSI problem" begin
    println("Solving FSI problem")
    xh0  = interpolate(xh,X_NSI(0.0))
    nls = NLSolver(
    show_trace = true,
    method = :newton,
    linesearch = BackTracking(),
    ftol = 1.0e-6,
    iterations = 50
    )
    odes =  ThetaMethod(nls, dt, θ)
    solver = TransientFESolver(odes)
    sol_NSI = solve(solver, op_FSI, xh0, t0, tf)
  end

  # Compute outputs
  out_params = Dict(
  :μ=>μ_f,
  "Um"=>Um,
  "⌀"=>⌀,
  :ρ=>ρ_f,
  "θ"=>θ,
  "model"=>model,
  "bdegree"=>bdegree,
  "trian"=>trian,
  "trian_Γc"=>trian_Γc,
  "quad_Γc"=>quad_Γc,
  "n_Γc"=>n_Γc,
  "xh0"=>xh0,
  "sol"=>sol_NSI,
  "filePath"=>filePath,
  "is_vtk"=>is_vtk,
  "coupling"=>coupling,
  )
  output = computeOutputs(problem,strategy;params=out_params)

end

function get_boundary_conditions(
  problem::FSIProblem{:oscillator},
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:weak},
  u_in,
  u_sym,
  u_cylinder,
  v_cylinder
  )

  u0 = VectorValue(0.0, 0.0)
  u0t(x,t) = u0
  u0t(t::Real) = x -> u0t(x,t)

  boundary_conditions = (
  # Tags
  NSI_Vw_tags = ["inlet", "noslip", "outlet"],
  NSI_Vu_tags = ["inlet", "noslip", "outlet"],
  NSI_Vv_tags = ["inlet", "noslip"],
  NSI_Vv_masks = [(true,true), (false,true)],
  ST_Vu_tags = ["inlet", "noslip", "cylinder","outlet"],
  ST_Vv_tags = ["inlet", "noslip", "cylinder"],
  ST_Vv_masks = [(true,true), (false,true), (true,true)],
  # Values,
  NSI_Vw_values = [u0, u0, u0],
  NSI_Vu_values = [u0t, u0t, u0t],
  NSI_Vv_values = [u_in, u_sym],
  ST_Vu_values = [u0, u0, u0, u0],
  ST_Vv_values = [u_in(0.0), u_sym(0.0), u_cylinder(0.0)],
  )
end

function get_boundary_conditions(
  problem::FSIProblem{:oscillator},
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:strong},
  u_in,
  u_sym,
  u_cylinder,
  v_cylinder
  )

  u0 = VectorValue(0.0, 0.0)
  u0t(x,t) = u0
  u0t(t::Real) = x -> u0t(x,t)

  boundary_conditions = (
  # Tags
  NSI_Vw_tags = ["inlet", "noslip", "cylinder", "outlet"],
  NSI_Vu_tags = ["inlet", "noslip", "cylinder", "outlet"],
  NSI_Vv_tags = ["inlet", "noslip", "cylinder"],
  NSI_Vv_masks = [(true,true), (false,true), (true,true)],
  ST_Vu_tags = ["inlet", "noslip", "cylinder","outlet"],
  ST_Vv_tags = ["inlet", "noslip", "cylinder"],
  ST_Vv_masks = [(true,true), (false,true), (true,true)],
  # Values,
  NSI_Vw_values = [u0, u0, u0, u0],
  NSI_Vu_values = [u0t, u0t, u_cylinder, u0t],
  NSI_Vv_values = [u_in, u_sym, v_cylinder],
  ST_Vu_values = [u0, u0, u_cylinder(0.0), u0],
  ST_Vv_values = [u_in(0.0), u_sym(0.0), v_cylinder(0.0)],
  )
end

function get_FE_spaces(
  problem::FSIProblem{:oscillator},
  strategy::WeakForms.MeshStrategy{:biharmonic},
  coupling::WeakForms.Coupling,
  model,
  order,
  bconds)

  Vw_NSI = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:NSI_Vw_tags]
    )
  Vu_NSI = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:NSI_Vu_tags]
    )
  Vv_NSI = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:NSI_Vv_tags],
    dirichlet_masks=bconds[:NSI_Vv_masks]
    )
  Vw_ST = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:ST_Vu_tags]
    )
  Vu_ST = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:ST_Vu_tags]
    )
  Vv_ST = TestFESpace(
    model=model,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=order,
    conformity =:H1,
    dirichlet_tags=bconds[:ST_Vv_tags],
    dirichlet_masks=bconds[:ST_Vv_masks]
    )
  Q = TestFESpace(
    model=model,
    valuetype=Float64,
    order=order-1,
    reffe=:Lagrangian,
    conformity=:C0
    )

  # Trial FE Spaces
  Uw_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values])
  Uw_NSI = TrialFESpace(Vw_NSI,bconds[:NSI_Vw_values])
  Uu_NSI = TransientTrialFESpace(Vu_NSI,bconds[:NSI_Vu_values])
  Uv_NSI = TransientTrialFESpace(Vv_NSI,bconds[:NSI_Vv_values])
  P = TrialFESpace(Q)

  # Multifield FE Spaces
  fe_spaces = (
    Y_ST = MultiFieldFESpace([Vw_ST,Vu_ST,Vv_ST,Q]),
    X_ST = MultiFieldFESpace([Uw_ST,Uu_ST,Uv_ST,P]),
    Y_NSI = MultiFieldFESpace([Vw_NSI,Vu_NSI,Vv_NSI,Q]),
    X_NSI = TransientMultiFieldFESpace([Uw_NSI,Uu_NSI,Uv_NSI,P])
  )
end
function computeOutputs(problem::FSIProblem{:oscillator},strategy::WeakForms.MeshStrategy;params=Dict())#, sol, xh0)

  # Unpack parameters
  model = params["model"]
  bdegree = params["bdegree"]
  xh0 = params["xh0"]
  sol = params["sol"]
  μ = params[:μ]
  trian = params["trian"]
  trian_Γc = params["trian_Γc"]
  quad_Γc = params["quad_Γc"]
  n_Γc = params["n_Γc"]
  Um = params["Um"]
  ⌀ = params["⌀"]
  ρ = params[:ρ]
  θ = params["θ"]
  filePath = params["filePath"]
  is_vtk = params["is_vtk"]
  coupling = params["coupling"]
  if( typeof(strategy) == WeakForms.MeshStrategy{:biharmonic} )
    uvpindex = [2,3,4]
  else
    uvpindex = [1,2,3]
  end

  #=
  # Aux function
  traction_boundary(n,u,v,p) = n ⋅ WeakForms.Pf_dev(μ,u,v)  + WeakForms.Pf_vol(u,p) * n
  function traction_closure(coupling,n,u,v,p)
    if typeof(coupling) == WeakForms.Coupling{:weak}
      traction_boundary(n,u,v,p).⁺
    else
      traction_boundary(n,u,v,p)
    end
  end
  traction_interface(n,u,v,p) = traction_closure(coupling,n,u,v,p)

  ## Initialize arrays
  tpl = Real[]
  FDpl = Real[]
  FLpl = Real[]
  uhn = xh0[uvpindex[1]]
  vhn = xh0[uvpindex[2]]
  phn = xh0[uvpindex[3]]
=#
  ## Loop over steps
  outfiles = paraview_collection(filePath, append=true
  ) do pvd
  for (i,(xh, t)) in enumerate(sol)
    println("STEP: $i, TIME: $t")
    println("============================")

    #=
    ## Get the solution at n+θ (where velocity and pressure are balanced)
    uh = xh[uvpindex[1]]
    vh = xh[uvpindex[2]]
    ph = xh[uvpindex[3]]
    uh_Γc = restrict(uh, trian_Γc)
    vh_Γc = restrict(vh, trian_Γc)
    ph_Γc = restrict(ph, trian_Γc)
    uhn_Γc = restrict(uhn, trian_Γc)
    vhn_Γc = restrict(vhn, trian_Γc)
    phn_Γc = restrict(phn, trian_Γc)
    uhθ_Γc = θ*uh_Γc + (1.0-θ)*uhn_Γc
    vhθ_Γc = θ*vh_Γc + (1.0-θ)*vhn_Γc
    phθ_Γc = θ*ph_Γc + (1.0-θ)*phn_Γc
    uh_Γi = restrict(uh, trian_Γi)
    vh_Γi = restrict(vh, trian_Γi)
    ph_Γi = restrict(ph, trian_Γi)
    uhn_Γi = restrict(uhn, trian_Γi)
    vhn_Γi = restrict(vhn, trian_Γi)
    phn_Γi = restrict(phn, trian_Γi)
    uhθ_Γi = θ*uh_Γi + (1.0-θ)*uhn_Γi
    vhθ_Γi = θ*vh_Γi + (1.0-θ)*vhn_Γi
    phθ_Γi = θ*ph_Γi + (1.0-θ)*phn_Γi

    # Integrate on the cylinder
    FDc, FLc = sum(integrate(traction_boundary(n_Γc,uhθ_Γc,vhθ_Γc,phθ_Γc), trian_Γc, quad_Γc))

    # Integrate on the interface
    FDi, FLi = sum(integrate(traction_interface(n_Γi,uhθ_Γi,vhθ_Γi,phθ_Γi), trian_Γi, quad_Γi))

    FD = FDc + FDi
    FL = FLc + FLi

    ## Drag and lift Forces
    push!(tpl, t)
    push!(FDpl, -FD)
    push!(FLpl, -FL)

    ## store step n
    uhn = uh
    vhn = vh
    phn = ph
=#
    # Write to PVD
    if(is_vtk)
      uh = xh[uvpindex[1]]
      vh = xh[uvpindex[2]]
      ph = xh[uvpindex[3]]
      pvd[t] = createvtk(
      trian,
      filePath * "_$t.vtu",
      cellfields = ["uh" => uh, "vh" => vh, "ph" => ph]
      )
    end
  end
end

#return (tpl, FDpl, FLpl)
end


text:
function PolynomialFunction(arg0::Vector{jdouble})
    return PolynomialFunction((Vector{jdouble},), arg0)
end

function add(obj::PolynomialFunction, arg0::PolynomialFunction)
    return jcall(obj, "add", PolynomialFunction, (PolynomialFunction,), arg0)
end

function anti_derivative(obj::PolynomialFunction)
    return jcall(obj, "antiDerivative", PolynomialFunction, ())
end

function degree(obj::PolynomialFunction)
    return jcall(obj, "degree", jint, ())
end

function equals(obj::PolynomialFunction, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_coefficients(obj::PolynomialFunction)
    return jcall(obj, "getCoefficients", Vector{jdouble}, ())
end

function hash_code(obj::PolynomialFunction)
    return jcall(obj, "hashCode", jint, ())
end

function integrate(obj::PolynomialFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "integrate", jdouble, (jdouble, jdouble), arg0, arg1)
end

function multiply(obj::PolynomialFunction, arg0::PolynomialFunction)
    return jcall(obj, "multiply", PolynomialFunction, (PolynomialFunction,), arg0)
end

function negate(obj::PolynomialFunction)
    return jcall(obj, "negate", PolynomialFunction, ())
end

function polynomial_derivative(obj::PolynomialFunction)
    return jcall(obj, "polynomialDerivative", PolynomialFunction, ())
end

function subtract(obj::PolynomialFunction, arg0::PolynomialFunction)
    return jcall(obj, "subtract", PolynomialFunction, (PolynomialFunction,), arg0)
end

function to_real_field_univariate_function(obj::FieldUnivariateFunction, arg0::Field)
    return jcall(obj, "toRealFieldUnivariateFunction", RealFieldUnivariateFunction, (Field,), arg0)
end

function to_string(obj::PolynomialFunction)
    return jcall(obj, "toString", JString, ())
end

function value(obj::PolynomialFunction, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::PolynomialFunction, arg0::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement,), arg0)
end

function value(obj::PolynomialFunction, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
using Test
using SphericalHarmonicModes

n1, n2 = HelioseismicKernels.Point2D(0, 0), HelioseismicKernels.Point2D(pi/4, 0);
n1′, n2′ = HelioseismicKernels.Point2D(pi/2, 0), HelioseismicKernels.Point2D(pi/2, pi/4);
r_src = HelioseismicKernels.r_src_default;
r_obs = HelioseismicKernels.r_obs_default;
xsrc1 = HelioseismicKernels.Point3D(r_src, n1);
xsrc2 = HelioseismicKernels.Point3D(r_src, n2);
xobs1 = HelioseismicKernels.Point3D(r_obs, n1);
xobs2 = HelioseismicKernels.Point3D(r_obs, n2);

@testset "all tests" begin
	@testset "rotation of Y₀₀" begin
		Y1′2′ = HelioseismicKernels.los_projected_biposh_spheroidal(HelioseismicKernels.computeY₀₀, n1′, n2′, HelioseismicKernels.los_earth(), 1:10)
		Y1′2′_12 = HelioseismicKernels.los_projected_biposh_spheroidal_Y₀₀(n1, n2, n1′, n2′, HelioseismicKernels.los_earth(), 1:10)
		@test Y1′2′ ≈  Y1′2′_12
	end
	@testset "rotated cross covariances" begin
		C = HelioseismicKernels.Cω(nothing, (n1, n1′), (n2, n2′), HelioseismicKernels.los_earth());
		C1 = HelioseismicKernels.Cω(nothing, n1, n2, HelioseismicKernels.los_earth());
		C2 = HelioseismicKernels.Cω(nothing, n1′, n2′, HelioseismicKernels.los_earth());
		@test C[:,1] ≈ C1
		@test C[:,2] ≈ C2
	end
	@testset "rotated biposh flipped" begin
		for los in Any[HelioseismicKernels.los_earth(), HelioseismicKernels.los_radial()]
			@testset "$los" begin
				SHModes = LM(0:10, 0:10)
				jₒjₛ_allmodes = L2L1Triangle(1:4, 3, 1:4)
				Y12, Y21, Y1′2′, Y2′1′ = HelioseismicKernels.los_projected_spheroidal_biposh_flippoints(
					(n1, n1′), (n2, n2′), los, SHModes, jₒjₛ_allmodes);
				Y12_2, Y21_2 = HelioseismicKernels.los_projected_spheroidal_biposh_flippoints(
					n1, n2, los, SHModes, jₒjₛ_allmodes);
				Y1′2′_2, Y2′1′_2 = HelioseismicKernels.los_projected_spheroidal_biposh_flippoints(
					n1′, n2′, los, SHModes, jₒjₛ_allmodes);

				@test all(all(isapprox(Y1[i], Y2[i], atol=1e-14, rtol=1e-8) for i in eachindex(Y1)) for (Y1,Y2) in zip(Y1′2′, Y1′2′_2))
				@test all(all(isapprox(Y1[i], Y2[i], atol=1e-14, rtol=1e-8) for i in eachindex(Y1)) for (Y1,Y2) in zip(Y2′1′, Y2′1′_2))
				@test all(all(isapprox(Y1[i], Y2[i], atol=1e-14, rtol=1e-8) for i in eachindex(Y1)) for (Y1,Y2) in zip(Y12, Y12_2))
				@test all(all(isapprox(Y1[i], Y2[i], atol=1e-14, rtol=1e-8) for i in eachindex(Y1)) for (Y1,Y2) in zip(Y21, Y21_2))
			end
		end
	end
	@testset "kernel components" begin
		@testset "SoundSpeed" begin
			function uniformvalidation(::HelioseismicKernels.SoundSpeed, m::HelioseismicKernels.SeismicMeasurement, xobs1, xobs2,
				los::HelioseismicKernels.los_direction = HelioseismicKernels.los_radial(); SHModes = LM(0:0, 0:0))

				ℓ_range = 1:1

				hω = HelioseismicKernels.hω(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range)

				Kδcₛₜ = HelioseismicKernels.kernel_δcₛₜ(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range, SHModes = SHModes, hω=hω, save = false);

				Kδc₀₀ = HelioseismicKernels.kernel_δc₀₀(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range, hω=hω, save = false);

				@test Kδcₛₜ[:,(0,0)]≈Kδc₀₀
			end

			@testset "TravelTimes" begin
				@testset "radial" begin
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_radial())
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_radial(),
						SHModes = LM(0:2, 0:2))
				end
				@testset "line-of-sight" begin
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_earth())
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_earth(),
						SHModes = LM(0:2, 0:2))
				end
			end
			@testset "Amplitudes" begin
				@testset "radial" begin
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_radial())
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_radial(),
						SHModes = LM(0:2, 0:2))
				end
				@testset "line-of-sight" begin
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_earth())
					uniformvalidation(HelioseismicKernels.SoundSpeed(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_earth(),
						SHModes = LM(0:2, 0:2))
				end
			end
		end

		@testset "flows" begin
			function uniformvalidation(::HelioseismicKernels.Flow, m::HelioseismicKernels.SeismicMeasurement, xobs1, xobs2,
				los::HelioseismicKernels.los_direction = HelioseismicKernels.los_radial();
				SHModes = LM(1:1, 0:0))

				ℓ_range = 1:1

				hω = HelioseismicKernels.hω(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range)

				Krθϕₛ₀ = HelioseismicKernels.kernel_uₛ₀_rθϕ(nothing, m, xobs1, xobs2, los; hω=hω, SHModes = SHModes,ℓ_range = ℓ_range, save = false)
				Krθϕₛ₀_2 = HelioseismicKernels.kernel_uₛ₀_rθϕ_2(nothing, m, xobs1, xobs2, los; hω=hω, SHModes = SHModes,ℓ_range = ℓ_range, save = false)

				Kℑu⁺₁₀ = HelioseismicKernels.kernel_ℑu⁺₁₀(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range, hω=hω, save = false);
				@testset "Ks0 from Kst" begin
					@test @view(Krθϕₛ₀[:,3,(1,0)]) ≈ Kℑu⁺₁₀
				end
				@testset "Ks0 directly" begin
					@test @view(Krθϕₛ₀_2[:,3,(1,0)]) ≈ Kℑu⁺₁₀
				end
			end

			@testset "TravelTimes" begin
				@testset "radial" begin
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_radial())
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_radial(),
						SHModes = LM(1:2, 0:2))
				end
				@testset "line-of-sight" begin
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_earth())
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.TravelTimes(), xobs1, xobs2, HelioseismicKernels.los_earth(),
						SHModes = LM(1:2, 0:2))
				end
			end

			@testset "Amplitudes" begin
				@testset "radial" begin
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_radial())
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_radial(),
						SHModes = LM(1:2, 0:2))
				end
				@testset "line-of-sight" begin
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_earth())
					uniformvalidation(HelioseismicKernels.Flow(), HelioseismicKernels.Amplitudes(), xobs1, xobs2, HelioseismicKernels.los_earth(),
						SHModes = LM(1:2, 0:2))
				end
			end
		end
	end;

	@testset "Ks0 computed directly and from Kst" begin

		function test(m, los; SHModes = LM(1:1, 0:0), ℓ_range = 2:2)
			hω = HelioseismicKernels.hω(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range)

			Krθϕₛ₀ = HelioseismicKernels.kernel_uₛ₀_rθϕ(nothing, m, xobs1, xobs2, los; hω=hω,
						SHModes = SHModes, ℓ_range = ℓ_range, save = false)
			Krθϕₛ₀_2 = HelioseismicKernels.kernel_uₛ₀_rθϕ_2(nothing, m, xobs1, xobs2, los; hω=hω,
						SHModes = SHModes, ℓ_range = ℓ_range, save = false)

			@test Krθϕₛ₀ ≈ Krθϕₛ₀_2
		end

		test(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_radial())
		test(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_radial(), SHModes = LM(0:4, 0:4))
		test(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_earth())
		test(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_earth(), SHModes = LM(0:4, 0:4))
	end;

	@testset "multiple points" begin
		@testset "flows" begin
			function testmultiplepoints(m, los)
				ℓ_range = 1:1
				hω = HelioseismicKernels.hω(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range)
				K = HelioseismicKernels.kernel_ℑu⁺₁₀(nothing, m, n1, n2, los,ℓ_range = ℓ_range, hω = hω,
					ν_ind_range = axes(hω, 1), save = false);
				K2 = HelioseismicKernels.kernel_ℑu⁺₁₀(nothing, m, n1,[n2], los,ℓ_range = ℓ_range, hω = hω,
					ν_ind_range = axes(hω, 1), save = false);
				@test K ≈ K2
			end

			@testset "TravelTimes" begin
				@testset "radial" begin
					testmultiplepoints(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_radial())
				end
				@testset "line-of-sight" begin
					testmultiplepoints(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_earth())
				end
			end
			@testset "Amplitudes" begin
				@testset "radial" begin
					testmultiplepoints(HelioseismicKernels.Amplitudes(), HelioseismicKernels.los_radial())
				end
				@testset "line-of-sight" begin
					testmultiplepoints(HelioseismicKernels.Amplitudes(), HelioseismicKernels.los_earth())
				end
			end
		end
		@testset "SoundSpeed" begin
			function testmultiplepoints(m, los)
				ℓ_range = 1:1
				hω = HelioseismicKernels.hω(nothing, m, xobs1, xobs2, los; ℓ_range = ℓ_range)
				K = HelioseismicKernels.kernel_δc₀₀(nothing, m, n1, n2, los, ℓ_range = ℓ_range, hω = hω,
					ν_ind_range = axes(hω, 1), save = false);
				K2 = HelioseismicKernels.kernel_δc₀₀(nothing, m, n1,[n2], los, ℓ_range = ℓ_range, hω = hω,
					ν_ind_range = axes(hω, 1), save = false);
				@test K ≈ K2
			end

			@testset "TravelTimes" begin
				@testset "radial" begin
					testmultiplepoints(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_radial())
				end
				@testset "line-of-sight" begin
					testmultiplepoints(HelioseismicKernels.TravelTimes(), HelioseismicKernels.los_earth())
				end
			end
			@testset "Amplitudes" begin
				@testset "radial" begin
					testmultiplepoints(HelioseismicKernels.Amplitudes(), HelioseismicKernels.los_radial())
				end
				@testset "line-of-sight" begin
					testmultiplepoints(HelioseismicKernels.Amplitudes(), HelioseismicKernels.los_earth())
				end
			end
		end
	end;
	@testset "rotated kernels" begin
		for los in Any[HelioseismicKernels.los_earth(), HelioseismicKernels.los_radial()]
			@testset "$los" begin
				K = HelioseismicKernels.kernel_uₛₜ(nothing, HelioseismicKernels.TravelTimes(), (n1, n1′), (n2, n2′), los, s_max = 1, t_max = 0, ℓ_range = 2:2, save = false);
				K2 = HelioseismicKernels.kernel_uₛₜ(nothing, HelioseismicKernels.TravelTimes(), n1′, n2′, los, s_max = 1, t_max = 0, ℓ_range = 2:2, save = false);
				@test view(K, :, :, :, 2) ≈ K2
			end
		end
	end;
end;


text:
"""
`ConvergenceTests` implements a few simple convergence tests for verifying
the expected numerical behavior of the spatial discretization
and time-stepping
"""
module ConvergenceTests

export test_rate_of_convergence

using Printf
using Statistics
using Logging

using JLD2
using OffsetArrays

using Oceananigans
using Oceananigans.Fields

using Oceananigans: Face, Center
using Oceananigans.Fields: Face, Center, nodes

import Oceananigans: RegularRectilinearGrid
import Oceananigans.Fields: location

figspath = joinpath(@__DIR__, "..", "figs")
ispath(figspath) || mkpath(figspath)

include("file_wrangling.jl")
include("analysis.jl")

# Exponential decay at a point:
#
# c = exp(-t) | ∂c/∂t = - c
#
# Tests time-stepper convergence.
#
include("PointExponentialDecay.jl")

# Utilities for analyzing and plotting the results
# of 1D convergence tests.
#
include("OneDimensionalUtils.jl")

# Advection and diffusion of a 1D Gaussian:
#
# c = 1/√(4πκt) * exp( -(x - Ut)^2 / 4κt )
#
# Tests x-advection and y-advection.
#
include("OneDimensionalGaussianAdvectionDiffusion.jl")

# Advection and diffusion of a 1D cosine:
#
# c = exp(-κt) * cos(x - Ut)
#
# Tests x-advection and y-advection.
#
include("OneDimensionalCosineAdvectionDiffusion.jl")

# Free decay of a 2D sinusoidal tracer patch
#
# c(x, y, t) = exp(-2t) * cos(x) * cos(y)
#
include("TwoDimensionalDiffusion.jl")

# Free decay of a horizontally-advected Taylor-Green vortex:
#
# u(x, y, t) = U + cos(x - Ut) sin(y)
# v(x, y, t) =   - sin(x - Ut) cos(y)
#
include("DoublyPeriodicTaylorGreen.jl")

include("ForcedFlowFreeSlip.jl")

include("ForcedFlowFixedSlip.jl")

function __init__()
    Logging.global_logger(OceananigansLogger())
end

end


text:
using Pkg
Pkg.activate(".")
pkg"dev .. ../CairoMakie ../GLMakie ../WGLMakie"
Pkg.instantiate()
Pkg.precompile()

Pkg.add(Pkg.PackageSpec(name="Franklin", rev="versions"));

using NodeJS
run(`$(npm_cmd()) install highlight.js`)
run(`$(npm_cmd()) install lunr`)
run(`$(npm_cmd()) install cheerio`)

using Franklin
using Documenter: deploydocs, deploy_folder, GitHubActions

cfg = GitHubActions() # this should pick up all details via GHA environment variables

# repo = "github.com/JuliaPlots/Makie.jl.git"
# push_preview = true
#
# deploydecision = deploy_folder(cfg; repo, push_preview,
#     devbranch = "master",
#     devurl = "dev",
# )

# @info "Setting PREVIEW_FRANKLIN_WEBSITE_URL to $repo"
# ENV["PREVIEW_FRANKLIN_WEBSITE_URL"] = repo
# @info "Setting PREVIEW_FRANKLIN_PREPATH to $(deploydecision.subfolder)"
ENV["PREVIEW_FRANKLIN_PREPATH"] = "Makie.jl"  # XXX since I forked it

serve(single=true, cleanup=false)
lunr()
optimize(minify=false, prerender=false, version="0.15.2")

# deploydocs(;
#     repo,
#     push_preview,
#     target = "__site",
# )


text:
module CbiomesProcessing

using JLD, MAT, YAML, SparseArrays, MeshArrays, Plots, Dates
using Distributed, DistributedArrays, FortranFiles

include("utilities.jl")
include("fileloops.jl")
include("examples.jl")

export StartWorkers, TaskDriver, MetaFileRead, MatrixInterp
export cbioproc_example1, cbioproc_example2

#export AverageYear

end # module

max_issues_repo_name: darwinproject/CBIOMES-Processing-legacy

max_forks_repo_name: darwinproject/CBIOMES-Processing-legacy


text:
#xortraintest-5.jl

#trains an xornet using noisy and unreliable data.  Ten data inputs, top two xor'ed to get the
#correct values, 5% of the time.

function xortrain_4()

  srand(10)

  println("working on unreliable xor data set with backpropagation")

  input_matrix = rand(Bool, 10, 500)
  training_results = Array{Bool,2}(1,500)
  training_results[:] = [input_matrix[1, col] $ input_matrix[2, col] $ (rand() < 0.05) for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{Float64,(10,2,1)}(randn)

  for rounds = 1:25
  for idx = 1:500
    input_column = input_matrix[:, idx]
    training_column = training_results[:, idx]

    GenML.Optimizers.backpropagationoptimize(xornet, input_column, training_column, GenML.CF.crossentropy)
  end
  end

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > 0.5) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for noisy xor: $wrongcount")

  return wrongcount
end

xortrain_4() == 0 || xortrain_4() == 0 || throw(ErrorException("two executions failed."))

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
using NonlinearSolve
using StaticArrays
using BenchmarkTools
using Test

function benchmark_immutable(f, u0)
    probN = NonlinearProblem{false}(f, u0)
    solver = init(probN, NewtonRaphson(), tol = 1e-9)
    sol = solve!(solver)
end

function benchmark_mutable(f, u0)
    probN = NonlinearProblem{false}(f, u0)
    solver = init(probN, NewtonRaphson(), immutable = false, tol = 1e-9)
    sol = (reinit!(solver, probN); solve!(solver))
end

function benchmark_scalar(f, u0)
    probN = NonlinearProblem{false}(f, u0)
    sol = (solve(probN, NewtonRaphson()))
end

function ff(u,p)
    u .* u .- 2
end
const cu0 = @SVector[1.0, 1.0]
function sf(u,p)
    u * u - 2
end
const csu0 = 1.0

sol = benchmark_immutable(ff, cu0)
@test sol.retcode === Symbol(NonlinearSolve.DEFAULT)
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = benchmark_mutable(ff, cu0)
@test sol.retcode === Symbol(NonlinearSolve.DEFAULT)
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = benchmark_scalar(sf, csu0)
@test sol.retcode === Symbol(NonlinearSolve.DEFAULT)
@test sol.u * sol.u - 2 < 1e-9

@test (@ballocated benchmark_immutable(ff, cu0)) == 0
@test (@ballocated benchmark_mutable(ff, cu0)) < 200
@test (@ballocated benchmark_scalar(sf, csu0)) == 0

# AD Tests
using ForwardDiff

# Immutable
f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]

g = function (p)
    probN = NonlinearProblem{false}(f, u0, p)
    sol = solve(probN, NewtonRaphson(), immutable = true, tol = 1e-9)
    return sol.u[end]
end

for p in 1.0:0.1:100.0
    @test g(p) ≈ sqrt(p)
    @test ForwardDiff.derivative(g, p) ≈ 1/(2*sqrt(p))
end

# Scalar
f, u0 = (u, p) -> u * u - p, 1.0

# NewtonRaphson
g = function (p)
    probN = NonlinearProblem{false}(f, oftype(p, u0), p)
    sol = solve(probN, NewtonRaphson())
    return sol.u
end

@test ForwardDiff.derivative(g, 1.0) ≈ 0.5

for p in 1.1:0.1:100.0
    @test g(p) ≈ sqrt(p)
    @test ForwardDiff.derivative(g, p) ≈ 1/(2*sqrt(p))
end

u0 = (1.0, 20.0)
# Falsi
g = function (p)
    probN = NonlinearProblem{false}(f, typeof(p).(u0), p)
    sol = solve(probN, Falsi())
    return sol.left
end

for p in 1.1:0.1:100.0
    @test g(p) ≈ sqrt(p)
    @test ForwardDiff.derivative(g, p) ≈ 1/(2*sqrt(p))
end

f, u0 = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
t = (p) -> [sqrt(p[2] / p[1])]
p = [0.9, 50.0]
for alg in [Bisection(), Falsi()]
    global g, p
    g = function (p)
        probN = NonlinearProblem{false}(f, u0, p)
        sol = solve(probN, Bisection())
        return [sol.left]
    end

    @test g(p) ≈ [sqrt(p[2] / p[1])]
    @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
end

gnewton = function (p)
    probN = NonlinearProblem{false}(f, 0.5, p)
    sol = solve(probN, NewtonRaphson())
    return [sol.u]
end
@test gnewton(p) ≈ [sqrt(p[2] / p[1])]
@test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p)

# Error Checks

f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0]
probN = NonlinearProblem(f, u0)

@test solve(probN, NewtonRaphson()).u[end] ≈ sqrt(2.0)
@test solve(probN, NewtonRaphson(); immutable = false).u[end] ≈ sqrt(2.0)
@test solve(probN, NewtonRaphson(;autodiff=false)).u[end] ≈ sqrt(2.0)
@test solve(probN, NewtonRaphson(;autodiff=false); immutable = false).u[end] ≈ sqrt(2.0)

for u0 in [1.0, [1, 1.0]]
  local f, probN, sol
  f = (u, p) -> u .* u .- 2.0
  probN = NonlinearProblem(f, u0)
  sol = sqrt(2) * u0

  @test solve(probN, NewtonRaphson()).u ≈ sol
  @test solve(probN, NewtonRaphson(); immutable = false).u ≈ sol
  @test solve(probN, NewtonRaphson(;autodiff=false)).u ≈ sol
  @test solve(probN, NewtonRaphson(;autodiff=false); immutable = false).u ≈ sol
end

# Bisection Tests
f, u0 = (u, p) -> u .* u .- 2.0, (1.0, 2.0)
probB = NonlinearProblem(f, u0)

# Falsi
solver = init(probB, Falsi())
sol = solve!(solver)
@test sol.left ≈ sqrt(2.0)

# this should call the fast scalar overload
@test solve(probB, Bisection()).left ≈ sqrt(2.0)

# these should call the iterator version
solver = init(probB, Bisection())
@test solver isa NonlinearSolve.BracketingImmutableSolver
# Question: Do we need BracketingImmutableSolver? We have fast scalar overload and
# Bracketing solvers work only for scalars.

solver = init(probB, Bisection(); immutable = false)
# @test solver isa NonlinearSolve.BracketingSolver
@test solve!(solver).left ≈ sqrt(2.0)

# Garuntee Tests for Bisection
f = function (u, p)
    if u < 2.0
        return u - 2.0
    elseif u > 3.0
        return u - 3.0
    else
        return 0.0
    end
end
probB = NonlinearProblem(f, (0.0, 4.0))

solver = init(probB, Bisection(;exact_left = true); immutable = false)
sol = solve!(solver)
@test f(sol.left, nothing) < 0.0
@test f(nextfloat(sol.left), nothing) >= 0.0

solver = init(probB, Bisection(;exact_right = true); immutable = false)
sol = solve!(solver)
@test f(sol.right, nothing) > 0.0
@test f(prevfloat(sol.right), nothing) <= 0.0

solver = init(probB, Bisection(;exact_left = true, exact_right = true); immutable = false)
sol = solve!(solver)
@test f(sol.left, nothing) < 0.0
@test f(nextfloat(sol.left), nothing) >= 0.0
@test f(sol.right, nothing) > 0.0
@test f(prevfloat(sol.right), nothing) <= 0.0


text:
mutable struct PaoSahModel <: TransistorModel
    mos::MOSStructure
    W::Number
    L::Number
    μ::AbstractMobilityModel
end

function Fpao(ψ,V,mod::PaoSahModel)
   if Na(mod.mos.semiconductor) > Nd(mod.mos.semiconductor)
    return (abs(ψ-V) < 1e-3) ? 0.0 : -q*(n(mod.mos,ψ-V))/(ξ(mod.mos,ψ))

   else
    return (abs(ψ-V) < 1e-3) ? 0.0 : q*(p(mod.mos,ψ))/(ξ(mod.mos,ψ))
   end
end

function FpaoC(ψ,V,mod::PaoSahModel)
    if Na(mod.mos.semiconductor) > Nd(mod.mos.semiconductor)
        #return μ(m,mod.mos.semiconductor,ψ)*ni(mod.mos)^2/abs(C(mod.mos))*exp((ψ-V)/kb/Temperature(mod.mos))/(ξ(mod.mos,ψ))
     return (abs(ψ-V) < 1e-3) ? 0.0 : -μ(mod.μ,mod.mos.semiconductor,ψ,V)* q*(n(mod.mos,ψ-V))/(ξ(mod.mos,ψ))
    else
     return (abs(ψ-V) < 1e-3) ? 0.0 : μ(mod.μ,mod.mos.semiconductor,ψ)*q*(p(mod.mos,ψ-V))/(ξ(mod.mos,ψ))
    end
 end

function id(Vg,Vd,Vs,Vb,mod::PaoSahModel)
    return hcubature(x->Fpao(x[1],x[2],mod), [0.0,Vs-Vb], [ψs(Vg,mod.mos),Vd];maxevals=500,reltol=1e-12)[1]
end

function Id(Vg,Vd,Vs,Vb,mod::PaoSahModel)
    return hcubature(x->FpaoC(x[1],x[2],mod), [0.0,Vs-Vb], [ψs(Vg,mod.mos),Vd];maxevals=500,reltol=1e-12)[1]
end

id(Vg,Vd,Vs,m::PaoSahModel) = id(Vg,Vd,Vs,Vs,m::PaoSahModel)
Id(Vg,Vd,Vs,m::PaoSahModel) = Id(Vg,Vd,Vs,Vs,m::PaoSahModel)
id(Vg,Vd,m::PaoSahModel) = id(Vg,Vd,0.0,m::PaoSahModel)
Id(Vg,Vd,m::PaoSahModel) = Id(Vg,Vd,0.0,m::PaoSahModel)
gm(Vg,Vd,Vs,Vb,m::PaoSahModel;order=5) = central_fdm(order,1)(x->Id(x,Vd,Vs,Vb,m),Vg)
gm(Vg,Vd,Vs,m::PaoSahModel;order=5) = central_fdm(order,1)(x->Id(x,Vd,Vs,m),Vg)
gm(Vg,Vd,m::PaoSahModel;order=5) = central_fdm(order,1)(x->Id(x,Vd,0.0,m),Vg)
gmId(Vg,Vd,m::PaoSahModel;order=5) = gm(Vg,Vd,m;order=order)/Id(Vd,Vg,m)
gds(Vg,Vd,m::PaoSahModel;order=5) = central_fdm(order,1)(x->Id(Vg,x,m),Vd)
gdsId(Vg,Vd,m::PaoSahModel;order=5) = gds(Vg,Vd,m;order=order)/Id(Vd,Vg,m)


text:
using Nabla

# Generate the values required to compute a matrix quadratic form.
N = 5
B = randn(5, 5)
A = B.'B + UniformScaling(1e-6)

# Low-level API computation of derivatives.

# Construct some trackable data.
x_ = randn(N)
x = Leaf(Tape(), x_)

# Compute the forward pass.
y = x.' * (A * x) # Temporary bracketting because we don't support RowVectors yet.
println("Output of the forward pass is:")
println(y)
println()
println("y is $(Nabla.unbox(y)).")
println()

# Get the reverse tape.
ȳ = ∇(y)
println("Output of reverse-pass is")
println(ȳ)
println()

# Index into the reverse tape using x to get the gradient of `y` w.r.t. `x`.
x̄ = ȳ[x]
println("Gradient of y w.r.t. x at $(Nabla.unbox(x)) is $x̄.")
println()

# (Current) High-Level API computation of derivatives. I will probably maintain this
# interface and extend is significantly as it is currently rather limited. It just returns
# a function which returns the gradient, which isn't really what you want.

# Define the function to be differentiated. Parameters w.r.t. which we want gradients must
# be arguments. Parameters that we don't want gradients w.r.t. should be passed in via a
# closure.
@unionise f(x::AbstractVector) = x.' * (A * x)

# Compute a function `∇f` which computes the derivative of `f` w.r.t. the input.
∇f = ∇(f)

# Compute the derivative of `f` w.r.t. `x` at `x_`. Result is currently a 1-Tuple. Might
# introduce a special case for unary functions where it just returns the result.
x̄ = ∇f(x_)

@assert x̄[1] == ȳ[x]


text:
module TestOperators

using JuMP
using LinearAlgebra
using SparseArrays
using Test

const MA = JuMP._MA

include(joinpath(@__DIR__, "utilities.jl"))
include(joinpath(@__DIR__, "JuMPExtension.jl"))

# For "DimensionMismatch when performing vector-matrix multiplication with custom types #988"
import Base: +, *
struct MyType{T}
    a::T
end
struct MySumType{T}
    a::T
end
Base.copy(t::MyType) = t
Base.zero(::Type{MyType{T}}) where {T} = MyType(zero(T))
Base.one(::Type{MyType{T}}) where {T} = MyType(one(T))
Base.zero(::Type{MySumType{T}}) where {T} = MySumType(zero(T))
Base.zero(::MySumType{T}) where {T} = MySumType(zero(T))
Base.transpose(t::MyType) = MyType(t.a)
Base.transpose(t::MySumType) = MySumType(t.a)
LinearAlgebra.adjoint(t::Union{MyType,MySumType}) = t
function +(
    t1::MyT,
    t2::MyS,
) where {MyT<:Union{MyType,MySumType},MyS<:Union{MyType,MySumType}}
    return MySumType(t1.a + t2.a)
end
*(t1::MyType{S}, t2::T) where {S,T} = MyType(t1.a * t2)
*(t1::S, t2::MyType{T}) where {S,T} = MyType(t1 * t2.a)
*(t1::MyType{S}, t2::MyType{T}) where {S,T} = MyType(t1.a * t2.a)

function JuMP.isequal_canonical(t::MySumType, s::MySumType)
    return JuMP.isequal_canonical(t.a, s.a)
end
function JuMP.isequal_canonical(
    x::AbstractArray{<:MySumType},
    y::AbstractArray{<:MySumType},
)
    return size(x) == size(y) && all(JuMP.isequal_canonical.(x, y))
end

function test_promotion(ModelType, VariableRefType)
    m = ModelType()
    I = Int
    V = VariableRefType
    A = JuMP.GenericAffExpr{Float64,VariableRefType}
    Q = JuMP.GenericQuadExpr{Float64,VariableRefType}
    @test promote_type(V, I) == A
    @test promote_type(I, V) == A
    @test promote_type(A, I) == A
    @test promote_type(I, A) == A
    @test promote_type(A, V) == A
    @test promote_type(V, A) == A
    @test promote_type(Q, I) == Q
    @test promote_type(I, Q) == Q
    @test promote_type(Q, A) == Q
    @test promote_type(A, Q) == Q
    @test promote_type(Q, V) == Q
    @test promote_type(V, Q) == Q
end

function test_broadcast_division_error(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2, 1:2])
    A = [1 2; 3 4]
    B = sparse(A)
    y = SparseMatrixCSC(2, 2, copy(B.colptr), copy(B.rowval), vec(x))
    @test_throws ErrorException A ./ x
    @test_throws ErrorException B ./ x
    @test_throws ErrorException A ./ y
    @test_throws ErrorException B ./ y

    # TODO: Refactor to avoid calling the internal JuMP function
    # `_densify_with_jump_eltype`.
    #z = JuMP._densify_with_jump_eltype((2 .* y) ./ 3)
    #@test JuMP.isequal_canonical((2 .* x) ./ 3, z)
    #z = JuMP._densify_with_jump_eltype(2 * (y ./ 3))
    #@test JuMP.isequal_canonical(2 .* (x ./ 3), z)
    #z = JuMP._densify_with_jump_eltype((x[1,1],) .* B)
    #@test JuMP.isequal_canonical((x[1,1],) .* A, z)
end

function test_vectorized_comparisons(ModelType, ::Any)
    m = ModelType()
    @variable(m, x[1:3])
    A = [
        1 2 3
        0 4 5
        6 0 7
    ]
    B = sparse(A)
    # force vector output
    cref1 = @constraint(m, reshape(x, (1, 3)) * A * x .>= 1)
    c1 = JuMP.constraint_object.(cref1)
    f1 = map(c -> c.func, c1)
    @test JuMP.isequal_canonical(
        f1,
        [
            x[1] * x[1] +
            2x[1] * x[2] +
            4x[2] * x[2] +
            9x[1] * x[3] +
            5x[2] * x[3] +
            7x[3] * x[3],
        ],
    )
    @test all(c -> c.set.lower == 1, c1)

    cref2 = @constraint(m, x' * A * x >= 1)
    c2 = JuMP.constraint_object.(cref2)
    @test JuMP.isequal_canonical(f1[1], c2.func)

    mat = [
        3x[1] + 12x[3] + 4x[2]
        2x[1] + 12x[2] + 10x[3]
        15x[1] + 5x[2] + 21x[3]
    ]

    cref3 = @constraint(m, (x'A)' + 2A * x .<= 1)
    c3 = JuMP.constraint_object.(cref3)
    f3 = map(c -> c.func, c3)
    @test JuMP.isequal_canonical(f3, mat)
    @test all(c -> c.set.upper == 1, c3)
    @test JuMP.isequal_canonical((x'A)' + 2A * x, (x'A)' + 2B * x)
    @test JuMP.isequal_canonical((x'A)' + 2A * x, (x'B)' + 2A * x)
    @test JuMP.isequal_canonical((x'A)' + 2A * x, (x'B)' + 2B * x)
    @test JuMP.isequal_canonical((x'A)' + 2A * x, MA.@rewrite((x'A)' + 2A * x))
    @test JuMP.isequal_canonical((x'A)' + 2A * x, MA.@rewrite((x'B)' + 2A * x))
    @test JuMP.isequal_canonical((x'A)' + 2A * x, MA.@rewrite((x'A)' + 2B * x))
    @test JuMP.isequal_canonical((x'A)' + 2A * x, MA.@rewrite((x'B)' + 2B * x))

    cref4 = @constraint(m, -1 .<= (x'A)' + 2A * x .<= 1)
    c4 = JuMP.constraint_object.(cref4)
    f4 = map(c -> c.func, c4)
    @test JuMP.isequal_canonical(f4, mat)
    @test all(c -> c.set.lower == -1, c4)
    @test all(c -> c.set.upper == 1, c4)

    cref5 = @constraint(m, -[1:3;] .<= (x'A)' + 2A * x .<= 1)
    c5 = JuMP.constraint_object.(cref5)
    f5 = map(c -> c.func, c5)
    @test JuMP.isequal_canonical(f5, mat)
    @test map(c -> c.set.lower, c5) == -[1:3;]
    @test all(c -> c.set.upper == 1, c4)

    cref6 = @constraint(m, -[1:3;] .<= (x'A)' + 2A * x .<= [3:-1:1;])
    c6 = JuMP.constraint_object.(cref6)
    f6 = map(c -> c.func, c6)
    @test JuMP.isequal_canonical(f6, mat)
    @test map(c -> c.set.lower, c6) == -[1:3;]
    @test map(c -> c.set.upper, c6) == [3:-1:1;]

    cref7 = @constraint(m, -[1:3;] .<= (x'A)' + 2A * x .<= 3)
    c7 = JuMP.constraint_object.(cref7)
    f7 = map(c -> c.func, c7)
    @test JuMP.isequal_canonical(f7, mat)
    @test map(c -> c.set.lower, c7) == -[1:3;]
    @test all(c -> c.set.upper == 3, c7)
end

function test_custom_dimension_mismatch(ModelType, VariableRefType)
    ElemT = MySumType{JuMP.GenericAffExpr{Float64,VariableRefType}}
    model = ModelType()
    @variable(model, Q[1:3, 1:3], PSD)
    x = [MyType(1), MyType(2), MyType(3)]
    y = Q * x
    z = x' * Q
    @test y isa Vector{ElemT}
    @test size(y) == (3,)
    @test z isa Adjoint{ElemT,<:Any}
    @test size(z) == (1, 3)
    for i in 1:3
        # Q is symmetric
        a = zero(JuMP.GenericAffExpr{Float64,VariableRefType})
        a += Q[1, i]
        a += 2Q[2, i]
        a += 3Q[3, i]
        # Q[1,i] + 2Q[2,i] + 3Q[3,i] is rearranged as 2 Q[2,3] + Q[1,3] + 3 Q[3,3]
        @test z[i].a == y[i].a == a
    end
end

function test_matrix_multiplication(ModelType, VariableRefType)
    ElemT = MySumType{JuMP.GenericAffExpr{Float64,VariableRefType}}
    model = ModelType()
    @variable(model, Q[1:3, 1:3], PSD)
    X = MyType.((1:3)' .+ (1:3))
    @test_expression Q * X
    Y = Q * X
    @test_expression X' * Q
    Z = X' * Q
    @test Y isa Matrix{ElemT}
    @test size(Y) == (3, 3)
    @test Z isa Matrix{ElemT}
    @test size(Z) == (3, 3)
    @test JuMP.isequal_canonical(Z', Y)
    @test_expression Q * X'
    Y = Q * X'
    @test_expression X * Q
    Z = X * Q
    @test Y isa Matrix{ElemT}
    @test size(Y) == (3, 3)
    @test Z isa Matrix{ElemT}
    @test size(Z) == (3, 3)
    @test JuMP.isequal_canonical(Z', Y)
end

function test_operator_warn(ModelType, ::Any)
    model = ModelType()
    @variable model x[1:51]
    # JuMPExtension does not have the `operator_counter` field
    if ModelType <: Model
        @test model.operator_counter == 0
    end
    # Triggers the increment of operator_counter since sum(x) has more than 50 terms
    @test_expression(sum(x) + 2x[1])
    if ModelType <: Model
        # The following check verifies that this test covers the code incrementing `operator_counter`
        @test model.operator_counter == 1
    end
end

function test_uniform_scaling(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    @test_expression_with_string x + 2I "x + 2"
    @test_expression_with_string (x + 1) + I "x + 2"
    @test_expression_with_string x - 2I "x - 2"
    @test_expression_with_string (x - 1) - I "x - 2"
    @test_expression_with_string 2I + x "x + 2"
    @test_expression_with_string I + (x + 1) "x + 2"
    @test_expression_with_string 2I - x "-x + 2"
    @test_expression_with_string I - (x - 1) "-x + 2"
    @test_expression_with_string I * x "x"
    @test_expression_with_string I * (x + 1) "x + 1"
    @test_expression_with_string (x + 1) * I "x + 1"
end

function test_basic_operators(ModelType, ::Any)
    model = ModelType()
    @variable(model, w)
    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    aff = @inferred 7.1 * x + 2.5
    @test_expression_with_string 7.1 * x + 2.5 "7.1 x + 2.5"
    aff2 = @inferred 1.2 * y + 1.2
    @test_expression_with_string 1.2 * y + 1.2 "1.2 y + 1.2"
    q = @inferred 2.5 * y * z + aff
    @test_expression_with_string 2.5 * y * z + aff "2.5 y*z + 7.1 x + 2.5"
    q2 = @inferred 8 * x * z + aff2
    @test_expression_with_string 8 * x * z + aff2 "8 x*z + 1.2 y + 1.2"
    @test_expression_with_string 2 * x * x + 1 * y * y + z + 3 "2 x² + y² + z + 3"
end

function test_basic_operators_number(ModelType, ::Any)
    model = ModelType()
    @variable(model, w)
    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    aff = @inferred 7.1 * x + 2.5
    q = @inferred 2.5 * y * z + aff

    # 1-1 Number--Number
    # ???
    # 1-2 Number--Variable
    @test_expression_with_string 4.13 + w "w + 4.13"
    @test_expression_with_string 3.16 - w "-w + 3.16"
    @test_expression_with_string 5.23 * w "5.23 w"
    @test_throws ErrorException 2.94 / w
    # 1-3 Number--AffExpr
    @test_expression_with_string 1.5 + aff "7.1 x + 4"
    @test_expression_with_string 1.5 - aff "-7.1 x - 1"
    @test_expression_with_string 2 * aff "14.2 x + 5"
    @test_throws ErrorException 2 / aff
    # 1-4 Number--QuadExpr
    @test_expression_with_string 1.5 + q "2.5 y*z + 7.1 x + 4"
    @test_expression_with_string 1.5 - q "-2.5 y*z - 7.1 x - 1"
    @test_expression_with_string 2 * q "5 y*z + 14.2 x + 5"
    @test_throws ErrorException 2 / q
end

function test_basic_operators_variable(ModelType, ::Any)
    model = ModelType()
    @variable(model, w)
    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    aff = @inferred 7.1 * x + 2.5
    q = @inferred 2.5 * y * z + aff

    # 2-0 Variable unary
    @test (+x) === x
    @test_expression_with_string -x "-x"
    # 2-1 Variable--Number
    @test_expression_with_string w + 4.13 "w + 4.13"
    @test_expression_with_string w - 4.13 "w - 4.13"
    @test_expression_with_string w * 4.13 "4.13 w"
    @test_expression_with_string w / 2.00 "0.5 w"
    @test w == w
    @test_expression_with_string x * y - 1 "x*y - 1"
    @test_expression_with_string x^2 "x²"
    @test_expression_with_string x^1 "x"
    @test_expression_with_string x^0 "1"
    @test_throws ErrorException x^3
    @test_throws ErrorException x^1.5
    # 2-2 Variable--Variable
    @test_expression_with_string w + x "w + x"
    @test_expression_with_string w - x "w - x"
    @test_expression_with_string w * x "w*x"
    @test_expression_with_string x - x "0"
    @test_throws ErrorException w / x
    @test_expression_with_string y * z - x "y*z - x"
    # 2-3 Variable--AffExpr
    @test_expression_with_string z + aff "z + 7.1 x + 2.5"
    @test_expression_with_string z - aff "z - 7.1 x - 2.5"
    @test_expression_with_string z * aff "7.1 z*x + 2.5 z"
    @test_throws ErrorException z / aff
    @test_throws MethodError z ≤ aff
    @test_expression_with_string 7.1 * x - aff "0 x - 2.5"
    # 2-4 Variable--QuadExpr
    @test_expression_with_string w + q "2.5 y*z + w + 7.1 x + 2.5"
    @test_expression_with_string w - q "-2.5 y*z + w - 7.1 x - 2.5"
    @test_throws ErrorException w * q
    @test_throws ErrorException w / q
    @test transpose(x) === x
    @test conj(x) === x
end

function test_basic_operators_affexpr(ModelType, ::Any)
    model = ModelType()
    @variable(model, w)
    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    aff = @inferred 7.1 * x + 2.5
    aff2 = @inferred 1.2 * y + 1.2
    q = @inferred 2.5 * y * z + aff
    # 3-0 AffExpr unary
    @test_expression_with_string +aff "7.1 x + 2.5"
    @test_expression_with_string -aff "-7.1 x - 2.5"
    # 3-1 AffExpr--Number
    @test_expression_with_string aff + 1.5 "7.1 x + 4"
    @test_expression_with_string aff - 1.5 "7.1 x + 1"
    @test_expression_with_string aff * 2 "14.2 x + 5"
    @test_expression_with_string aff / 2 "3.55 x + 1.25"
    @test_throws MethodError aff ≤ 1
    @test aff == aff
    @test_throws MethodError aff ≥ 1
    @test_expression_with_string aff - 1 "7.1 x + 1.5"
    @test_expression_with_string aff^2 "50.41 x² + 35.5 x + 6.25"
    @test_expression_with_string (7.1 * x + 2.5)^2 "50.41 x² + 35.5 x + 6.25"
    @test_expression_with_string aff^1 "7.1 x + 2.5"
    @test_expression_with_string (7.1 * x + 2.5)^1 "7.1 x + 2.5"
    @test_expression_with_string aff^0 "1"
    @test_expression_with_string (7.1 * x + 2.5)^0 "1"
    @test_throws ErrorException aff^3
    @test_throws ErrorException (7.1 * x + 2.5)^3
    @test_throws ErrorException aff^1.5
    @test_throws ErrorException (7.1 * x + 2.5)^1.5
    # 3-2 AffExpr--Variable
    @test_expression_with_string aff + z "7.1 x + z + 2.5"
    @test_expression_with_string aff - z "7.1 x - z + 2.5"
    @test_expression_with_string aff * z "7.1 x*z + 2.5 z"
    @test_throws ErrorException aff / z
    @test_expression_with_string aff - 7.1 * x "0 x + 2.5"
    # 3-3 AffExpr--AffExpr
    @test_expression_with_string aff + aff2 "7.1 x + 1.2 y + 3.7"
    @test_expression_with_string aff - aff2 "7.1 x - 1.2 y + 1.3"
    @test_expression_with_string aff * aff2 "8.52 x*y + 3 y + 8.52 x + 3"
    @test string((x + x) * (x + 3)) == string((x + 3) * (x + x))  # Issue #288
    @test_throws ErrorException aff / aff2
    @test_expression_with_string aff - aff "0 x"
    # 4-4 AffExpr--QuadExpr
    @test_expression_with_string aff2 + q "2.5 y*z + 1.2 y + 7.1 x + 3.7"
    @test_expression_with_string aff2 - q "-2.5 y*z + 1.2 y - 7.1 x - 1.3"
    @test_throws ErrorException aff2 * q
    @test_throws ErrorException aff2 / q
    @test transpose(aff) === aff
    @test conj(aff) === aff
end

function test_basic_operators_quadexpr(ModelType, ::Any)
    model = ModelType()
    @variable(model, w)
    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    aff = @inferred 7.1 * x + 2.5
    aff2 = @inferred 1.2 * y + 1.2
    q = @inferred 2.5 * y * z + aff
    q2 = @inferred 8 * x * z + aff2
    # 4-0 QuadExpr unary
    @test_expression_with_string +q "2.5 y*z + 7.1 x + 2.5"
    @test_expression_with_string -q "-2.5 y*z - 7.1 x - 2.5"
    # 4-1 QuadExpr--Number
    @test_expression_with_string q + 1.5 "2.5 y*z + 7.1 x + 4"
    @test_expression_with_string q - 1.5 "2.5 y*z + 7.1 x + 1"
    @test_expression_with_string q * 2 "5 y*z + 14.2 x + 5"
    @test_expression_with_string q / 2 "1.25 y*z + 3.55 x + 1.25"
    @test q == q
    @test_expression_with_string aff2 - q "-2.5 y*z + 1.2 y - 7.1 x - 1.3"
    # 4-2 QuadExpr--Variable
    @test_expression_with_string q + w "2.5 y*z + 7.1 x + w + 2.5"
    @test_expression_with_string q - w "2.5 y*z + 7.1 x - w + 2.5"
    @test_throws ErrorException q * w
    @test_throws ErrorException q / w
    # 4-3 QuadExpr--AffExpr
    @test_expression_with_string q + aff2 "2.5 y*z + 7.1 x + 1.2 y + 3.7"
    @test_expression_with_string q - aff2 "2.5 y*z + 7.1 x - 1.2 y + 1.3"
    @test_throws ErrorException q * aff2
    @test_throws ErrorException q / aff2
    # 4-4 QuadExpr--QuadExpr
    @test_expression_with_string q + q2 "2.5 y*z + 8 x*z + 7.1 x + 1.2 y + 3.7"
    @test_expression_with_string q - q2 "2.5 y*z - 8 x*z + 7.1 x - 1.2 y + 1.3"
    @test_throws ErrorException q * q2
    @test_throws ErrorException q / q2
    @test transpose(q) === q
    @test conj(q) === q
end

function test_dot(ModelType, ::Any)
    model = ModelType()
    @variable(model, 0 ≤ x[1:3] ≤ 1)

    @test_expression_with_string dot(x[1], x[1]) "x[1]²"
    @test_expression_with_string dot(2, x[1]) "2 x[1]"
    @test_expression_with_string dot(x[1], 2) "2 x[1]"

    c = vcat(1:3)
    @test_expression_with_string dot(c, x) "x[1] + 2 x[2] + 3 x[3]"
    @test_expression_with_string dot(x, c) "x[1] + 2 x[2] + 3 x[3]"

    A = [1 3; 2 4]
    @variable(model, 1 ≤ y[1:2, 1:2] ≤ 1)
    @test_expression_with_string dot(A, y) "y[1,1] + 2 y[2,1] + 3 y[1,2] + 4 y[2,2]"
    @test_expression_with_string dot(y, A) "y[1,1] + 2 y[2,1] + 3 y[1,2] + 4 y[2,2]"

    B = ones(2, 2, 2)
    @variable(model, 0 ≤ z[1:2, 1:2, 1:2] ≤ 1)
    @test_expression_with_string dot(B, z) "z[1,1,1] + z[2,1,1] + z[1,2,1] + z[2,2,1] + z[1,1,2] + z[2,1,2] + z[1,2,2] + z[2,2,2]"
    @test_expression_with_string dot(z, B) "z[1,1,1] + z[2,1,1] + z[1,2,1] + z[2,2,1] + z[1,1,2] + z[2,1,2] + z[1,2,2] + z[2,2,2]"

    @objective(model, Max, dot(x, ones(3)) - dot(y, ones(2, 2)))
    for i in 1:3
        JuMP.set_start_value(x[i], 1)
    end
    for i in 1:2, j in 1:2
        JuMP.set_start_value(y[i, j], 1)
    end
    for i in 1:2, j in 1:2, k in 1:2
        JuMP.set_start_value(z[i, j, k], 1)
    end
    @test dot(c, JuMP.start_value.(x)) ≈ 6
    @test dot(A, JuMP.start_value.(y)) ≈ 10
    @test dot(B, JuMP.start_value.(z)) ≈ 8

    if ModelType <: Model
        # Only `Model` is guaranteed to have `operator_counter`, so
        # only test for that case.
        @testset "dot doesn't trigger operator_counter" begin
            # Check that dot is not falling back to default, inefficient
            # addition (JuMP PR #943).
            model = ModelType()
            @test model.operator_counter == 0
            @variable(model, x[1:100])
            JuMP.set_start_value.(x, 1:100)
            @expression(model, test_sum, sum(x[i] * i for i in 1:100))
            @expression(model, test_dot1, dot(x, 1:100))
            @expression(model, test_dot2, dot(1:100, x))
            @test model.operator_counter == 0
            test_add = test_dot1 + test_dot2
            @test model.operator_counter == 1  # Check triggerable.
            test_sum_value = JuMP.value(test_sum, JuMP.start_value)
            @test test_sum_value ≈ JuMP.value(test_dot1, JuMP.start_value)
            @test test_sum_value ≈ JuMP.value(test_dot2, JuMP.start_value)
        end
    end
end

function test_higher_level(ModelType, ::Any)
    model = ModelType()
    @variable(model, 0 ≤ matrix[1:3, 1:3] ≤ 1, start = 1)
    @testset "sum(j::DenseAxisArray{Variable})" begin
        @test_expression_with_string sum(matrix) "matrix[1,1] + matrix[2,1] + matrix[3,1] + matrix[1,2] + matrix[2,2] + matrix[3,2] + matrix[1,3] + matrix[2,3] + matrix[3,3]"
    end

    @testset "sum(j::DenseAxisArray{T}) where T<:Real" begin
        @test sum(JuMP.start_value.(matrix)) ≈ 9
    end
    @testset "sum(j::Array{VariableRef})" begin
        @test string(sum(matrix[1:3, 1:3])) == string(sum(matrix))
    end
    @testset "sum(affs::Array{AffExpr})" begin
        @test_expression_with_string sum([
            2 * matrix[i, j] for i in 1:3, j in 1:3
        ]) "2 matrix[1,1] + 2 matrix[2,1] + 2 matrix[3,1] + 2 matrix[1,2] + 2 matrix[2,2] + 2 matrix[3,2] + 2 matrix[1,3] + 2 matrix[2,3] + 2 matrix[3,3]"
    end
    @testset "sum(quads::Array{QuadExpr})" begin
        @test_expression_with_string sum([
            2 * matrix[i, j]^2 for i in 1:3, j in 1:3
        ]) "2 matrix[1,1]² + 2 matrix[2,1]² + 2 matrix[3,1]² + 2 matrix[1,2]² + 2 matrix[2,2]² + 2 matrix[3,2]² + 2 matrix[1,3]² + 2 matrix[2,3]² + 2 matrix[3,3]²"
    end
    S = [1, 3]
    @variable(model, x[S], start = 1)
    @testset "sum(j::JuMPDict{VariableRef})" begin
        @test_expression sum(x)
        @test length(string(sum(x))) == 11 # order depends on hashing
        @test occursin("x[1]", string(sum(x)))
        @test occursin("x[3]", string(sum(x)))
    end
    @testset "sum(j::JuMPDict{T}) where T<:Real" begin
        @test sum(JuMP.start_value.(x)) == 2
    end
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        f = getfield(@__MODULE__, name)
        @testset "$(name)" begin
            f(Model, VariableRef)
        end
        @testset "$(name)-JuMPExtension" begin
            f(JuMPExtension.MyModel, JuMPExtension.MyVariableRef)
        end
    end
end

end

TestOperators.runtests()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct Rr
    V::Array{Float64}
    TtDY::Array{Float64}
    sv::Vector{Float64}
    lb::Float64
    xmeans::Vector{Float64}
    ymeans::Vector{Float64}
    weights::Vector{Float64}
end

"""
    rr(X, Y, weights = ones(size(X, 1)); lb = .01)
Ridge regression (RR) implemented by SVD factorization.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `lb` : A value of the regularization parameter "lambda".

`X` and `Y` are internally centered. The model is computed with an intercept.

## References

Cule, E., De Iorio, M., 2012. A semi-automatic method to guide the choice
of ridge parameter in ridge regression. arXiv:1205.0686.

Hastie, T., Tibshirani, R., 2004. Efficient quadratic regularization
for expression arrays. Biostatistics 5, 329-340. https://doi.org/10.1093/biostatistics/kxh010

Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning: data mining,
inference, and prediction, 2nd ed. Springer, New York.

Hoerl, A.E., Kennard, R.W., 1970. Ridge Regression: Biased Estimation for Nonorthogonal Problems.
Technometrics 12, 55-67. https://doi.org/10.1080/00401706.1970.10488634
"""
function rr(X, Y, weights = ones(size(X, 1)); lb = .01)
    rr!(copy(X), copy(Y), weights; lb = lb)
end

function rr!(X, Y, weights = ones(size(X, 1)); lb = .01)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    sqrtw = sqrt.(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    sqrtD = Diagonal(sqrtw)
    res = LinearAlgebra.svd!(sqrtD * X)
    sv = res.S
    TtDY = Diagonal(sv) * res.U' * (sqrtD * Y)
    Rr(res.V, TtDY, sv, lb, xmeans, ymeans, weights)
end

"""
    rrchol(X, Y, weights = ones(size(X, 1)); lb = .01)
Ridge regression (RR) using the Normal equations and a Cholesky factorization.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `lb` : A value of the regularization parameter "lambda".

`X` and `Y` are internally centered. The model is computed with an intercept.

## References

Cule, E., De Iorio, M., 2012. A semi-automatic method to guide the choice
of ridge parameter in ridge regression. arXiv:1205.0686.

Hastie, T., Tibshirani, R., 2004. Efficient quadratic regularization
for expression arrays. Biostatistics 5, 329-340. https://doi.org/10.1093/biostatistics/kxh010

Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning: data mining,
inference, and prediction, 2nd ed. Springer, New York.

Hoerl, A.E., Kennard, R.W., 1970. Ridge Regression: Biased Estimation for Nonorthogonal Problems.
Technometrics 12, 55-67. https://doi.org/10.1080/00401706.1970.10488634
"""
function rrchol(X, Y, weights = ones(size(X, 1)); lb = .01)
    rrchol!(copy(X), copy(Y), weights; lb = lb)
end

function rrchol!(X, Y, weights = ones(size(X, 1)); lb = .01)
    @assert size(X, 2) > 1 "Method only working for X with > 1 column."
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    p = size(X, 2)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    XtD = X' * Diagonal(weights)
    B = cholesky!(Hermitian(XtD * X + lb^2 * Diagonal(ones(p)))) \ (XtD * Y)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    coef(object::Rr; lb = nothing)
Compute the b-coefficients of a fitted model.
* `object` : The fitted model.
* `lb` : A value of the regularization parameter "lambda".
    If nothing, it is the parameter stored in the fitted model.
"""
function coef(object::Rr; lb = nothing)
    isnothing(lb) ? lb = object.lb : nothing
    eig = object.sv.^2
    z = 1 ./ (eig .+ lb^2)
    beta = Diagonal(z) * object.TtDY
    B = object.V * beta
    int = object.ymeans' .- object.xmeans' * B
    tr = sum(eig .* z)
    (B = B, int = int, df = 1 + tr)
end

"""
    predict(object::Rr, X; lb = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
* `lb` : Regularization parameter, or collection of regularization parameters,
    "lambda" to consider. If nothing, it is the parameter stored in the
    fitted model.
"""
function predict(object::Rr, X; lb = nothing)
    isnothing(lb) ? lb = object.lb : nothing
    le_lb = length(lb)
    pred = list(le_lb, Matrix{Float64})
    @inbounds for i = 1:le_lb
        z = coef(object; lb = lb[i])
        pred[i] = z.int .+ X * z.B
    end
    le_lb == 1 ? pred = pred[1] : nothing
    (pred = pred,)
end


text:
const Code = Union{AbstractVector{<:Integer}, CodeUnits}
const Sequence = Union{AbstractString, Code}

"""
    NGramIterator{T}

Iterates over ngram codes of collection of integers `s` using [`Mill.string_start_code()`](@ref) and [`Mill.string_end_code()`](@ref) for padding. NGram codes are computed as in positional number systems, where items of `s` are digits, `b` is the base, and `m` is modulo.

In order to reduce collisions when mixing ngrams of different order one should avoid zeros and negative integers in `s` and should set base `b` to the expected number of unique tokens in `s`.

See also: [`NGramMatrix`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref).
"""
struct NGramIterator{T}
    s::T
    n::Int
    b::Int
    m::Int

    NGramIterator(s::T, n::Int=3, b::Int=256, m::Int=typemax(Int)) where T <: Maybe{Code} = new{T}(s, n, b, m)
end

"""
    NGramIterator(s, n=3, b=256, m=typemax(Int))

Construct an [`NGramIterator`](@ref). If `s` is an `AbstractString` it is first converted
to integers with `Base.codeunits`.

# Examples
```jldoctest
julia> NGramIterator("deadbeef", 3, 256, 17) |> collect
10-element Vector{Int64}:
  2
 16
  9
  9
  6
 10
 11
 15
  2
  6

julia> NGramIterator(collect(1:9), 3, 10, 1009) |> collect
11-element Vector{Int64}:
 221
 212
 123
 234
 345
 456
 567
 678
 789
 893
 933

julia> Mill.string_start_code()
0x02

julia> Mill.string_end_code()
0x03
```

See also: [`NGramMatrix`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref).
"""
NGramIterator(s::AbstractString, args...) = NGramIterator(codeunits(s), args...)

Base.eltype(::NGramIterator) = Int

@inline _len(s::Sequence, n) = length(s) + n - 1
Base.length(it::NGramIterator) = _len(it.s, it.n)
Base.length(it::NGramIterator{Missing}) = 0

_next_ngram(z, i, it::NGramIterator) = _next_ngram(z, i, it.s, it.n, it.b)
function _next_ngram(z, i, s::Code, n, b)
    z *= b
    z += i ≤ length(s) ? s[i] : string_end_code()
    z -= b^n * (i > n ? s[i - n] : string_start_code())
end

Base.iterate(it::NGramIterator{Missing}) = nothing
function Base.iterate(it::NGramIterator, (z, i) = (_init_z(it), 1))
    if i ≤ length(it)
        z = _next_ngram(z, i, it)
        z % it.m, (z, i+1)
    end
end

_init_z(it) = _init_z(it.n, it.b)
function _init_z(n::T, b::T) where T
    z = zero(T)
    for _ in 1:n
        z *= b
        z += one(T)
    end
    string_start_code() * z
end

Base.hash(it::NGramIterator{T}, h::UInt) where {T} = hash((string(T), it.s, it.n, it.b, it.m), h)
(it1::NGramIterator{T} == it2::NGramIterator{T}) where {T} = it1.s == it2.s &&
    it1.n === it2.n && it1.b === it2.b && it1.m === it2.m

"""
    ngrams!(o, x, n=3, b=256)

Store codes of `n` grams of `x` using base `b` to `o`.

# Examples
```jldoctest
julia> o = zeros(Int, 5)
5-element Vector{Int64}:
 0
 0
 0
 0
 0

julia> ngrams!(o, "foo", 3, 256)
5-element Vector{Int64}:
  131686
  157295
 6713199
 7302915
 7275267
```

See also: [`ngrams`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref), [`NGramMatrix`](@ref), [`NGramIterator`](@ref).
"""
function ngrams!(o, x::Sequence, args...)
    for (i, idx) in enumerate(NGramIterator(x, args...))
        o[i] = idx
    end
    o
end

"""
    ngrams(o, x, n=3, b=256)

Return codes of `n` grams of `x` using base `b`.

# Examples
```jldoctest
julia> ngrams("foo", 3, 256)
5-element Vector{Int64}:
  131686
  157295
 6713199
 7302915
 7275267
```

See also: [`ngrams!`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref), [`NGramMatrix`](@ref), [`NGramIterator`](@ref).
"""
ngrams(x::Sequence, args...) = collect(NGramIterator(x, args...))

"""
    countngrams!(o, x, n, b, m=length(o))

Count the number of of `n` grams of `x` using base `b` and modulo `m` and store the result to `o`.

# Examples
```jldoctest
julia> o = zeros(Int, 5)
5-element Vector{Int64}:
 0
 0
 0
 0
 0

julia> countngrams!(o, "foo", 3, 256)
5-element Vector{Int64}:
 2
 1
 1
 0
 1
```

See also: [`countngrams`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref),
    [`NGramMatrix`](@ref), [`NGramIterator`](@ref).
"""
countngrams!(o, x::Sequence, n::Int, b::Int) = countngrams!(o, x, n, b, length(o))
function countngrams!(o, x::Sequence, args...)
    for idx in NGramIterator(x, args...)
        o[idx + 1] += 1
    end
    o
end

"""
    countngrams(o, x, n, b, m)

Count the number of of `n` grams of `x` using base `b` and modulo `m` into a vector of length `m`
in case `x` is a single sequence or into a matrix with `m` rows if `x` is an iterable of sequences.

# Examples
```jldoctest
julia> countngrams("foo", 3, 256, 5)
5-element Vector{Int64}:
 2
 1
 1
 0
 1

julia> countngrams(["foo", "bar"], 3, 256, 5)
5×2 Matrix{Int64}:
 2  1
 1  0
 1  2
 0  0
 1  2
```

See also: [`countngrams!`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref),
    [`NGramMatrix`](@ref), [`NGramIterator`](@ref).
"""
countngrams(x, n, b, m) = countngrams!(zeros(Int, m), x, n, b)
function countngrams(x::AbstractVector{<:Sequence}, n, b, m)
    o = zeros(Int, m, length(x))
    for (i, s) in enumerate(x)
        countngrams!(view(o, :, i), x[i], n, b, m)
    end
    o
end

"""
    NGramMatrix{T, U, V} <: AbstractMatrix{U}

A matrix-like structure for lazily representing sequences like strings as ngrams of
cardinality `n` using `b` as a base for calculations and `m` as the modulo. Therefore, the matrix
has `m` rows and one column for representing each sequence. Missing sequences are supported.

See also: [`NGramIterator`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref).
"""
struct NGramMatrix{T, U, V} <: AbstractMatrix{V}
    S::U
    n::Int
    b::Int
    m::Int

    function NGramMatrix(S::AbstractVector{T}, n::Int=3, b::Int=256, m::Int=2053) where T <: Sequence
        new{T, typeof(S), Int}(S, n, b, m)
    end

    function NGramMatrix(S::AbstractVector{Missing}, n::Int=3, b::Int=256, m::Int=2053)
        new{Missing, typeof(S), Missing}(S, n, b, m)
    end

    function NGramMatrix(S::AbstractVector{T}, n::Int=3, b::Int=256, m::Int=2053) where T <: Maybe{Sequence}
        new{T, typeof(S), Maybe{Int}}(S, n, b, m)
    end

    function NGramMatrix{T, U, V}(S, n::Int=3, b::Int=256, m::Int=2053) where
            {T <: Maybe{Sequence}, U <: AbstractVector{T}, V <: Maybe{Int}}
        new{T, U, V}(convert(U, S), n, b, m)
    end
end

"""
    NGramMatrix(s, n=3, b=256, m=2053)

Construct an [`NGramMatrix`](@ref). `s` can either be a single sequence or any `AbstractVector`.

# Examples
```jldoctest
julia> NGramMatrix([1,2,3])
2053×1 NGramMatrix{Vector{Int64}, Vector{Vector{Int64}}, Int64}:
 [1, 2, 3]

julia> NGramMatrix(["a", missing, "c"], 2, 128)
2053×3 NGramMatrix{Union{Missing, String}, Vector{Union{Missing, String}}, Union{Missing, Int64}}:
 "a"
 missing
 "c"
```

See also: [`NGramIterator`](@ref), [`ngrams`](@ref), [`ngrams!`](@ref), [`countngrams`](@ref),
    [`countngrams!`](@ref).
"""
NGramMatrix(s::Maybe{Sequence}, args...) = NGramMatrix([s], args...)

NGramIterator(A::NGramMatrix, i::Integer) = NGramIterator(A.S[i], A.n, A.b, A.m)

Base.length(A::NGramMatrix) = A.m * length(A.S)
Base.size(A::NGramMatrix) = (A.m, length(A.S))
Base.size(A::NGramMatrix, d) = (d == 1) ? A.m : length(A.S)

Base.getindex(X::NGramMatrix, idcs...) = (@boundscheck checkbounds(X, idcs...); _getindex(X, idcs...))
_getindex(X::NGramMatrix, ::Colon, i::Integer) = NGramMatrix(X.S[[i]], X.n, X.b, X.m)
_getindex(X::NGramMatrix, ::Colon, i::AbstractArray) = NGramMatrix(X.S[i], X.n, X.b, X.m)
_getindex(X::NGramMatrix, ::Colon, ::Colon) = NGramMatrix(X.S[:], X.n, X.b, X.m)

subset(a::NGramMatrix, i) = NGramMatrix(a.S[i], a.n, a.b, a.m)

function Base.convert(::Type{<:NGramMatrix{T, U}}, A::NGramMatrix) where {T, U <: AbstractVector{T}}
    NGramMatrix(convert(U, A.S), A.n, A.b, A.m)
end

function Base.promote_rule(::Type{NGramMatrix{T, U, V}}, ::Type{NGramMatrix{A, B, C}}) where
    {T, A, U, B, V, C}
    NGramMatrix{promote_type(T, A), promote_type(U, B), promote_type(V, C)}
end
function Base.promote_rule(::Type{<:NGramMatrix{Missing, U, V}}, ::Type{<:NGramMatrix{A, B, C}}) where
    {A <: Sequence, U, B, V, C}
    X = Union{Missing, A}
    NGramMatrix{X, promote_type(U, B){X}, promote_type(V, C)}
end
Base.promote_rule(t1::Type{<:NGramMatrix{<:Sequence}}, t2::Type{<:NGramMatrix{Missing}}) = promote_rule(t2, t1)

function _check_nbm(As::AbstractVecOrTuple{NGramMatrix})
    n, b, m = As[1].n, As[1].b, As[1].m
    if any(!isequal(n), (A.n for A in As)) ||
        any(!isequal(b), (A.b for A in As)) ||
        any(!isequal(m), (A.m for A in As))
        DimensionMismatch(
                          "Matrices do not have the same n, b, or m."
                         ) |> throw
    end
    n, b, m
end

Base.hcat(As::T...) where T <: NGramMatrix = _typed_hcat(T, As)
Base.hcat(As::NGramMatrix...) = _typed_hcat(_promote_types(As...), As)
function _typed_hcat(::Type{T}, As::Tuple{Vararg{NGramMatrix}}) where T <: NGramMatrix
    T(vcat([A.S for A in As]...), _check_nbm(As)...)
end

Base.reduce(::typeof(hcat), As::Vector{<:NGramMatrix}) = _typed_hcat(mapreduce(typeof, promote_type, As), As)
function _typed_hcat(::Type{T}, As::AbstractVector{<:NGramMatrix}) where T <: NGramMatrix
    T(reduce(vcat, [A.S for A in As]), _check_nbm(As)...)
end

SparseArrays.SparseMatrixCSC(x::NGramMatrix) = SparseArrays.SparseMatrixCSC{Int64, UInt}(x)
function SparseArrays.SparseMatrixCSC{Tv, Ti}(x::NGramMatrix) where {Tv, Ti <: Integer}
    size(x, 2) == 0 && return sparse(Ti[],Ti[],Tv[], size(x,1), size(x,2))
    l = sum(map(i -> length(NGramIterator(x, i)), 1:size(x, 2)))
    I = zeros(Ti, l)
    J = zeros(Ti, l)
    V = ones(Tv, l)
    vid = 1
    for j in 1:size(x, 2)
        for i in NGramIterator(x, j)
            I[vid] = i + 1
            J[vid] = j
            vid += 1
        end
    end
    sparse(I, J, V, size(x,1), size(x,2))
end

A::AbstractMatrix * B::NGramMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::AbstractMatrix * B::NGramMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))

_mul(A::AbstractMatrix, B::NGramMatrix{Missing}) = fill(missing, size(A, 1), size(B, 2))

# TODO rewrite this to less parameters once Zygote allows for composite grads
_mul(A::AbstractMatrix, B::NGramMatrix{T}) where T <: Maybe{Sequence} = _mul(A, B.S, B.n, B.b, B.m)

function _mul(A::AbstractMatrix, S::AbstractVector{T}, n, b, m) where T <: Maybe{Sequence}
    T_res = Missing <: T ? Union{eltype(A), Missing} : eltype(A)
    C = zeros(T_res, size(A, 1), length(S))
    z = _init_z(n, b)
    for (k, s) in enumerate(S)
        _mul_vec!(C, k, A, z, s, n, b, m)
    end
    C
end

Zygote.@adjoint function _mul(A::AbstractMatrix, S::AbstractVector{<:Maybe{Sequence}}, n, b, m)
    return _mul(A, S, n, b, m), Δ -> (_mul_∇A(Δ, A, S, n, b, m), nothing, nothing, nothing, nothing)
end

function _mul_∇A(Δ, A, S, n, b, m)
    ∇A = zero(A)
    z = _init_z(n, b)
    for (k, s) in enumerate(S)
        _∇A_mul_vec!(Δ, k, ∇A, z, s, n, b, m)
    end
    ∇A
end
function ChainRulesCore.rrule(::typeof(_mul_∇A), Δ, A, S, n, b, m)
    y = _mul_∇A(Δ, A, S, n, b, m)
    function _mul_∇A_pullback(Δ₂)
        return (NoTangent(), _mul_∇₂A(Δ₂, Δ, A, S, n, b, m)...)
    end
    y, _mul_∇A_pullback
end

function _mul_∇₂A(Δ₂, Δ, A, S, n, b, m)
    _mul(Δ₂, S, n, b, m), nothing, nothing, nothing, nothing, nothing
end

_mul_vec!(C, k, A, z, s::AbstractString, args...) = _mul_vec!(C, k, A, z, codeunits(s), args...)
function _mul_vec!(C, k, A, z, s, n, b, m, ψ=nothing)
    for l in 1:_len(s, n)
        z = _next_ngram(z, l, s, n, b)
        zm = z%m + 1
        for i in 1:size(C, 1)
            @inbounds C[i, k] += A[i, zm]
        end
    end
end
_mul_vec!(C, k, A, z, s::Missing, n, b, m, ψ=missing) = @inbounds C[:, k] .= ψ

_∇A_mul_vec!(Δ, k, ∇A, z, s::AbstractString, args...) = _∇A_mul_vec!(Δ, k, ∇A, z, codeunits(s), args...)
function _∇A_mul_vec!(Δ, k, ∇A, z, s, n, b, m)
    for l in 1:_len(s, n)
        z = _next_ngram(z, l, s, n, b)
        zm = z%m + 1
        for i in 1:size(∇A, 1)
            @inbounds ∇A[i, zm] += Δ[i, k]
        end
    end
end
_∇A_mul_vec!(Δ, k, ∇A, z, s::Missing, n, b, m) = return

Base.hash(M::NGramMatrix, h::UInt) = hash((M.S, M.n, M.b, M.m), h)
(M1::NGramMatrix == M2::NGramMatrix) = isequal(M1.S == M2.S, true) && M1.n == M2.n && M1.b == M2.b && M1.m == M2.m
Base.isequal(M1::NGramMatrix, M2::NGramMatrix) = isequal(M1.S, M2.S) && M1.n == M2.n && M1.b == M2.b && M1.m == M2.m