text:
module RulesBIFMNextTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:BIFM:znext" begin

     @testset "Belief Propagation: (m_out::MultivariateNormalDistributionsFamily, m_in::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, meta::BIFMMeta)" begin

        meta = BIFMMeta([2.0 0; 0 1], # A
                        [3.0 0; 0 2], # B
                        [4.0 0; 0 3], # C
                        [5.0 0; 0 4], # H
                        [6.0 0; 0 5], # BHBt
                        [1.0, 2],     # ξz
                        [7.0 0; 0 6], # Λz
                        [2.0, 3],     # ξztilde
                        [8.0 0; 0 7], # Λztilde
                        [3.0, 4],     # μu
                        [9.0 0; 0 8]  # Σu
                        )

        @test_rules [ with_float_conversions = true, float32_atol = 1e-2, bigfloat_atol = 1e-8 ] BIFM(:znext, Marginalisation) [
            (input = (m_out = MvNormalMeanPrecision([1,2], [2 0; 0 1]),     m_in = MvNormalMeanPrecision([1,2], [1 0; 0 2]), m_zprev = ProdFinal(MvNormalMeanPrecision([1,2], [1 0; 0 2])), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-61.0, -48.0], [6730.0 0.0; 0.0 425.5]))),
            (input = (m_out = MvNormalMeanPrecision([3,4], [6 0; 0 1]),     m_in = MvNormalMeanPrecision([8,2], [2 0; 0 2]), m_zprev = ProdFinal(MvNormalMeanPrecision([1,6], [2 0; 0 2])), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-49, -164], [3368 0.0; 0.0 425.5]))),
            (input = (m_out = MvNormalMeanPrecision([5,6], [2 0; 0 5]),     m_in = MvNormalMeanPrecision([1,9], [1 0; 0 1]), m_zprev = ProdFinal(MvNormalMeanPrecision([6,2], [1 0; 0 1])), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-471.0, -28.0], [6730.0 0.0; 0.0 846.0])))
        ]

    end

end

end


text:
using FunctionWrappers

# when to switch aggregators from tuples for rates/affects to FunctionWrappers
const TUPLE_TO_FWRAPPER_CUTOFF = 10

function get_jump_info_tuples(constant_jumps)
    rates    = ((c.rate for c in constant_jumps)...)
    affects! = ((c.affect! for c in constant_jumps)...)

    return rates, affects!
  end

  function get_jump_info_fwrappers(u, p, t, constant_jumps)
    RateWrapper   = FunctionWrappers.FunctionWrapper{typeof(t),Tuple{typeof(u), typeof(p), typeof(t)}}
    rates         = [RateWrapper(c.rate) for c in constant_jumps]
    AffectWrapper = FunctionWrappers.FunctionWrapper{Void,Tuple{Any}}
    affects!      = [AffectWrapper(x->(c.affect!(x);nothing)) for c in constant_jumps]

    return rates, affects!
  end

ext: jl

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
"Square the sum of the numbers up to the given number" # These descriptions were the wrong way round in the skeleton
# I knew GCSE maths would come in handy one day
# a = 1 + 2 + ... + n
# b = n + ... + 2 + 1
# a+b = n(n+1)
# => sum = n(n+1) ÷ 2 # even * odd = even
function square_of_sum(n::Integer)
#   These two lines are just as fast as each other: it turns out that Julia knows GCSE maths too
#   (see sum(::Abstractrange{<:Real}) in base/range.jl)
#   (n*(n+1) ÷ 2) ^ 2
    sum(Base.OneTo(n))^2
end

# a = 1^2 + 2^2 + ... + n
# b = n^2 + (n-1)^2 + ...
#
# Telescoping series trick:
# sum(k^2 - (k-1)^2) = n^2
# k^2 - (k^2 - 2k +1) = 2k - 1
# => sum(k^2 - (k-1)^2) = sum(2k - 1) = 2*sum(k) - n = n^2 => sum(k) = (n^2 + n) / 2
#
# We want sum([telescoping thing]) ∝ (sum(k^2) + tractable stuff)
# => try sum(k^3 - (k-1)^3) = n^3
# = sum(k^3 - (k-1)(k^2 - 2k +1)) = sum(k^3 - (k^3 -3k^2 +3k -1)) = sum(3k^2 - 3k +1) = 3*sum(k^2) - 3*sum(k) + n = n^3
# => sum(k^2) = 1/3 * (n^3 - n + 3(n(n+1))/2) = 1/3 n^3 + 1/2 n^2 + 1/6 n = n/6 * (2n^2 + 3n + 1) = n/6 * (2n + 1)(n + 1)
#
# Slightly faster way of getting this:
# ```
# using SymPy
# @vars(k,n)
# summation(k^2,(k,(1,n)))
# ```
#
# Can turn it into a function with summation(...).evalf() but it uses Python so it's slow - better off copying it out manually
"Sum the squares of the numbers up to the given number" # These descriptions were the wrong way round in the skeleton
function sum_of_squares(n::Integer)
    n*(2n+1)*(n+1) ÷ 6
end

"Subtract sum of squares from square of sums"
function difference(n::Integer)
    square_of_sum(n) - sum_of_squares(n)
end


text:
using VLKeggSDK
using DataFrames

# setup an array of ec numbers -
ec_number_array = [
    "ec:2.7.1.2"   ;
    "ec:5.3.1.9"   ;
    "ec:2.7.1.11"  ;
    "ec:4.1.2.13"  ;
]

# download these reactions -
reaction_array = get_reactions_for_ec_number(ec_number_array) |> check

# build data table -
data_table = DataFrame(reaction_array)


text:
import Knet: update!

warmup_cosine(x, warmup=0.002) = x < warmup ? x/warmup : 0.5 * (1.0 + cos(π * x))
warmup_constant(x, warmup=0.002) = x < warmup ? x/warmup : 1.0
warmup_linear(x, warmup=0.002) = x < warmup ? x/warmup : 1.0 - x

mutable struct BertAdam
    lr::AbstractFloat
    beta1::AbstractFloat
    beta2::AbstractFloat
    eps::AbstractFloat
    t::Int
    gclip::AbstractFloat
    fstm
    scndm
    w_decay_rate::AbstractFloat
    schedule
    warmup
    t_total
end

BertAdam(; lr=0.001, gclip=1.0, beta1=0.9, beta2=0.999, eps=1e-6, w_decay_rate=0.0, schedule="warmup_linear", warmup=-1, t_total=-1)=BertAdam(lr, beta1, beta2, eps, 0, gclip, nothing, nothing, w_decay_rate, schedule, warmup, t_total)

for T in (Array{Float32},Array{Float64},KnetArray{Float32},KnetArray{Float64}); @eval begin
    function update!(w::$T, g::$T, p::BertAdam)
        Knet.gclip!(g, p.gclip)
        if p.fstm===nothing; p.fstm=zero(w); p.scndm=zero(w); end
        lmul!(p.beta1, p.fstm)
        axpy!(1-p.beta1, g, p.fstm)
        lmul!(p.beta2, p.scndm)
        axpy!(1-p.beta2, g .* g, p.scndm)
        # They don't do bias correction for some reason
        #fstm_corrected = p.fstm / (1 - p.beta1 ^ p.t)
        #scndm_corrected = p.scndm / (1 - p.beta2 ^ p.t)
        if p.t_total !== -1
            schedule_func = eval(Meta.parse(p.schedule))
            lr_scheduled = p.lr * schedule_func(p.t/p.t_total, p.warmup)
        else
            lr_scheduled = p.lr
        end

        if p.w_decay_rate > 0.0
            axpy!(-lr_scheduled, (p.fstm ./ (sqrt.(p.scndm) .+ p.eps)) .+ (p.w_decay_rate * w), w)
        else
            axpy!(-lr_scheduled, (p.fstm ./ (sqrt.(p.scndm) .+ p.eps)), w)
        end

        p.t += 1
    end
end;end

max_issues_repo_name: OsmanMutlu/BERT-in-Knet

max_forks_repo_name: OsmanMutlu/BERT-in-Knet


text:
using RHEOSafm
using RHEOS
using Test

include("importJPK.jl")
include("preprocessing.jl")


text:
using QuantumOptics: DenseOperator
using LinearAlgebra: diagm, norm
using Test, IonSim
using Suppressor

Base.:(*)(s::String, i::Int) = i==1 ? s : 0
Base.:(*)(i::Int, s::String) = Base.:(*)(s, i)
Base.:(+)(s::String, i::Int) = i==0 ? s : nothing
Base.:(+)(i::Int, s::String) = Base.:(+)(s, i)

function get_indices(ion_participation, vibrational_mode_dimension, two_modes; cutoff=[Inf,Inf])
    N = vibrational_mode_dimension[1]
    σ1 = [0 0; "a" 0]; σ2 = [0 0; "b" 0]; I = [1 0; 0 1]
    c = ["c$i$j," for i in 1:N, j in 1:N]
    if two_modes
        M = vibrational_mode_dimension[2]
        d = ["d$i$j," for i in 1:M, j in 1:M]
        subarrays = [[d, c, I, σ1], [d, c, σ2, I]]
    else
        subarrays = [[c, I, σ1], [c, σ2, I]]
    end
    if ion_participation == 1
        k = kron(subarrays[1]...)
    elseif ion_participation == 2
        k = kron(subarrays[2]...)
    elseif ion_participation == 12
        k = kron(subarrays[1]...) + kron(subarrays[2]...)
    end
    ridxs = Dict()
    cidxs = Dict()
    nzidxs = []
    for i in 1:size(k, 2), j in 1:size(k, 2)
        if k[i, j] == 0
            continue
        end
        mn = map(x->[parse(Int,x[2]),parse(Int,x[3])], split(k[i,j], ",")[1:end-1])
        if sum([abs(x[1] - x[2]) > cutoff[l] for (l,x) in enumerate(mn)]) > 0
            continue
        end
        if i > j
            if !(k[i, j] in keys(ridxs))
                ridxs[k[i, j]] = [(i,j)]
            else
                push!(ridxs[k[i,j]], (i,j))
            end
        else
            if !(k[i, j] in keys(cidxs))
                parity = sum(map(x->isodd(abs(x[1]-x[2])), mn))
                cidxs[k[i, j]] = [(-1 * isodd(parity),0), (i,j)]
            else
                push!(cidxs[k[i,j]], (i,j))
            end
        end
    end
    collect(values(ridxs)), collect(values(cidxs))
end

@suppress_err begin

@testset "hamiltonians -- misc functions" begin
    # make sure _D and _D_constant_eta return the same thing for constant eta
    L = 10
    Ω = 1e6randn()
    Δ = 10randn()
    η = [abs(randn()) for _ in 1:L]
    ν = [10randn() for _ in 1:L]
    timescale = 10randn()
    n = [[rand(1:20) for _ in 1:L], [rand(1:5) for _ in 1:L]]
    t = randn()
    d1 = IonSim._D(Ω, Δ, η, ν, timescale, n, t, L)
    D = [IonSim._Dnm_cnst_eta(η[i], n[1][i], n[2][i]) for i in 1:L]
    d2 = IonSim._D_cnst_eta(Ω, Δ, ν, timescale, n, D, t, L)
    @test abs(sum(d1 .- d2)) < 1e-8

    # _flattenall (this function should completely flatten nested arrays)
    a = [[[1], [2], [3]], [[4], [5], [6]]]
    @test IonSim._flattenall(a) == collect(1:6)

    # _get_kron_indxs
    a = ["a$i$j," for i in 1:3, j in 1:3]
    b = ["b$i$j," for i in 1:3, j in 1:3]
    c = ["c$i$j," for i in 1:3, j in 1:3]
    abc = kron(a, b, c)
    z = 1:3
    i1 = rand(z); i2 = rand(z)
    j1 = rand(z); j2 = rand(z)
    k1 = rand(z); k2 = rand(z)
    ijk = [(i1, i2), (j1, j2), (k1, k2)]
    (I, J) = IonSim._get_kron_indxs([(i1, i2), (j1, j2), (k1, k2)], [3, 3, 3])
    @test [(parse(Int,i[2]), parse(Int,i[3])) for i in split(abc[I, J], ",")[1:end-1]] == ijk
    # ^this tests that  abc[I, J] = a[i1,i2] * a[j1,j2] * c[k1,k2]

    # _inv_get_kron_indxs
    indxs = IonSim._inv_get_kron_indxs((I, J), (3,3,3))
    indxs = collect(zip(indxs[1], indxs[2]))
    @test indxs == ijk
end

@testset "hamiltonians -- parameter arrays" begin
    # setup system
    C = Ca40()
    C1 = copy(C)
    L1 = Laser(); L1.pointing = [(1, 1.0), (2, 1.0)]
    L2 = copy(L1)

    # _Ωmatrix
    L1.E = 1
    L1.ϕ = 0
    L2.E = 2
    L2.ϕ = 2
    chain = LinearChain(
            ions=[C, C1], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
        )
    T = Trap(configuration=chain, lasers=[L1, L2])
    Ωnmkj = IonSim._Ωmatrix(T, 1)
    t = rand(0:1e-3:100)

    # coupling strength between ion1-laser1 and ion2-laser2 should be identical for all
    # transitions since the ions are identical
    resolve(x, t) = typeof(x)<:Number ? x : x(t)
    @test [resolve(i, t) for i in Ωnmkj[1, 1]] == [resolve(i, t) for i in Ωnmkj[2, 1]]

    # coupling strength between ion1-laser1 and ion1-laser2 should be proportional according
    # to the factor 2exp(-2im) due to the differences we've set in L2.E and L2.ϕ
    @test [2exp(-2im)*resolve(i,t) for i in Ωnmkj[1, 1]] == [resolve(i,t) for i in Ωnmkj[1, 2]]
    @test [2exp(-2im)*resolve(i,t) for i in Ωnmkj[2, 1]] == [resolve(i,t) for i in Ωnmkj[2, 2]]

    # make sure time-dep L.E and L.ϕ propagate appropriately
    L1.E = cos; L1.ϕ = t -> t^2
    Ωnmkj = IonSim._Ωmatrix(T, 1)
    t = 0:1e-3:100
    for Ω in Ωnmkj[1, 1]
        if typeof(Ω)<:Number
            continue
        end
        @test Ω.(t) ≈ @.(Ω(0) * cos(t) * exp(-im * t^2))
    end

    # _Δmatrix
    # zero B, zero laser detuning, zero stark shift should give an array of zeros
    Δ = IonSim._Δmatrix(T, 1)
    for i in 1:2, j in 1:2
        @test Δ[i, j][1] == 0
    end

    # zero B, zero stark shift, non-zero laser detuning should return columns that all have
    # the same value
    L1.Δ = 1; L2.Δ = -1
    Δ = IonSim._Δmatrix(T, 1)
    @test Δ[1, 1][1] ≈ 2π && Δ[2, 1][1] ≈ 2π && Δ[1, 2][1] ≈ -2π && Δ[2, 2][1] ≈ -2π

    # zero B, zero laser detuning, now add stark shift to just one of the ions
    L1.Δ = 0; L2.Δ = 0
    C.stark_shift["S-1/2"] = 1
    Δ = IonSim._Δmatrix(T, 1)
    @test Δ[1, 1][1] ≈ -2π && Δ[1, 2][1] ≈ -2π && Δ[2, 1][1] ≈ 0 && Δ[2, 2][1] ≈ 0

    # lastly let's test when resonant
    C.selected_level_structure = ["S-1/2", "D-1/2"]
    T.B = 1e-4
    f = transition_frequency(T, C, ("S-1/2", "D-1/2"))
    L1.Δ = f
    Δ = IonSim._Δmatrix(T, 1)
    @test Δ[1, 1][1] ≈ 0

    # _ηmatrix
    chain = LinearChain(
                ions=[C, C1], com_frequencies=(x=2e6,y=2e6,z=1e6),
                vibrational_modes=(x=[1], y=[2], z=[1])
            )
    L1.k = (x̂ + ẑ)/√2
    L2.k = (ŷ + ẑ)/√2
    L3 = Laser(); L3.pointing = [(1, 1.0), (2, 1.0)]
    L3.k = ẑ
    T = Trap(configuration=chain, lasers=[L1, L2, L3])
    η = IonSim._ηmatrix(T)

    # η[1,1,1] corresponds laser1, ion1, mode1 (x̂) and η[1,1,3] corresponds laser1, ion1,
    # mode3 (ẑ). They both have the same projection on L1.k, but the frequency of mode1 is
    # twice as large as mode2, so it's L-D factor should be √2 times smaller
    @test abs(η[1, 1, end](1.0)) ≈ abs(η[1, 1, end-2](1.0) / √2)
    # With this setup, the L-D factor should be the same for ion1 and ion2
    @test η[1, 1, end](1.0) ≈ η[2, 1, end](1.0)
    # η[1, 2, 2] is the 1st ion, 2nd laser and y-stretch-mode. The L-D factor should be
    # opposite in sign, equal in magnitude to the 2nd ion, 2nd laser and y-stretch-mode
    @test η[1, 2, end-1](1.0) ≈ -η[2, 2, end-1](1.0)
    # L3, which is in the ẑ direction should only have projection on zmode (mode3)
    @test η[1, 3, end] ≡ 0
    @test η[1, 3, end-1] ≡ 0
    @test η[1, 3, end-2](0.0) != 0
    # test construction of time-dep δν. If δν = 1e6*t, then after 3e-6 seconds (and since
    # ν=1e6), √(ν+δν(t)) = √2 * √(ν) = √2 * √(ν+δν(0))
    chain.vibrational_modes.z[1].δν = t -> 1e6t
    η = IonSim._ηmatrix(T)
    @test η[1, 1, end-2](0.0) ≈ 2 * η[1, 1, end-2](3)
end

@testset "hamiltonians -- setup functions" begin
    # tests for _setup_global_B_hamiltonian

    # setup system
    C = Ca40()
    L = Laser()
    chain = LinearChain(
        ions=[C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
    )
    T = Trap(configuration=chain, lasers=[L], δB=0)
    global_B_indices, global_B_scales, bfunc = IonSim._setup_global_B_hamiltonian(T, 1)

    # T.δB = 0 -> global_B_indices and global_B_scales should be empty arrays
    @test length(global_B_indices) == 0 && length(global_B_scales) == 0
    # and bfunc should be identically zero like T.δB
    @test sum(bfunc.(0:1e4)) == 0

    # now let's test nontrivial T.δB
    T.δB = sin
    t = 0:0.1:10
    T.configuration.vibrational_modes.z[1].N = 3
    global_B_indices, global_B_scales, bfunc = IonSim._setup_global_B_hamiltonian(T, 1)
    # make sure bfunc is correct
    @test bfunc.(t) == 2π .* sin.(t)
    # make sure all energy levels are being recorded
    @test sort(IonSim._flattenall(global_B_indices)) == collect(1:32)
    # and make sure that the susceptibilites are correct
    zs = [zeeman_shift(1, sls) for sls in values(C.selected_level_structure)]
    @test length(unique([global_B_scales; zs])) == length(zs)

    # test _setup_δν_hamiltonian

    # setup system
    chain = LinearChain(
        ions=[C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
    )
    T = Trap(configuration=chain, lasers=[L], δB=0)
    # should return empty arrays if δν=0
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    @test length(δν_indices) == 0 && length(δν_functions) == 0

    # test output for simple case
    T.configuration.vibrational_modes.z[1].N = 3
    T.configuration.vibrational_modes.z[1].δν = 1
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    indxs = [[8j + i for i in 1:8] for j in 1:3]
    @test δν_indices[1] == indxs

    # add another mode with δν=0 and test output
    chain = LinearChain(
        ions=[C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(y=[1], z=[1])
    )
    T = Trap(configuration=chain, lasers=[L], δB=0)
    T.configuration.vibrational_modes.y[1].N = 3
    T.configuration.vibrational_modes.z[1].N = 3
    T.configuration.vibrational_modes.z[1].δν = 1
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    indxs1 = [[8*4*j + i for i in 1:8*4] for j in 1:3]
    @test δν_indices[1] == indxs1

    # test output when both modes have nonzero δν
    T.configuration.vibrational_modes.z[1].δν = 1
    T.configuration.vibrational_modes.y[1].δν = 1
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    indxs = [[8*j + i for i in 1:8*8] for j in 1:3]
    @test δν_indices[1][1] == [vcat([[8*(4(j-1)+1) + i for i in 1:8] for j in 1:4]...);]
    @test δν_indices[2] == indxs1
    # we have two modes but only one nontrivial δν, so length(δν_functions) should equal 1
    @test length(δν_functions) == 2

    # finally, make sure δν_functions are being constructed appropriately
    T.configuration.vibrational_modes.y[1].δν = cos
    T.configuration.vibrational_modes.z[1].δν = sin
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    t = 0:0.1:10
    @test δν_functions[1].(t) == 2π .* cos.(t)
    @test δν_functions[2].(t) == 2π .* sin.(t)

    # _setup_fluctuation_hamiltonian
    T = Trap(configuration=chain, lasers=[L], δB=1)
    T.configuration.vibrational_modes.y[1].δν = cos
    T.configuration.vibrational_modes.z[1].δν = sin
    δν_indices, δν_functions = IonSim._setup_δν_hamiltonian(T, 1)
    global_B_indices, global_B_scales, bfunc = IonSim._setup_global_B_hamiltonian(T, 1)
    all_unique_indices, gbi, gbs, bfunc1, δνi, δνfuncs = IonSim._setup_fluctuation_hamiltonian(T, 1)
    # make sure information from _setup_global_B_hamiltonian and _setup_δν_hamiltonian
    # is propagated appropriately
    @test any(gbi .== global_B_indices)
    @test any(gbs .== global_B_scales)
    @test bfunc(17.0) == bfunc1(17.0)
    @test any(δνi .== δν_indices)
    @test δνfuncs[end](92.0) == δν_functions[end](92.0)

    # _setup_base_hamiltonian
    # Let's make sure that _setup_base_hamiltonian is recording the appropriate indices.
    C = Ca40(["S-1/2", "D-1/2"])
    L = Laser()
    chain = LinearChain(
            ions=[C, C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
        )
    T = Trap(configuration=chain, B=4e-4, Bhat=(x̂ + ŷ + ẑ)/√3, lasers=[L])
    mode = T.configuration.vibrational_modes.z[1]
    mode.N = rand(1:8)
    N = mode.N + 1
    Efield_from_rabi_frequency!(1e6, T, 1, 1, ("S-1/2", "D-1/2"))

    ## first just shine light on 1st ion
    L.pointing = [(1, 1.0)]
    ridxs, cidxs = get_indices(1, [N], false)
    _, r, c = IonSim._setup_base_hamiltonian(T, 1e-6, 100, Inf, "analytic", true)
    @test length(unique([ridxs; r])) == (N^2 + N) / 2 == length(ridxs)
    @test length(unique([cidxs; c]))-1 == (N^2 - N) / 2 == length(cidxs)

    ## now just shine light on second ion
    L.pointing = [(2, 1.0)]
    ridxs, cidxs = get_indices(2, [N], false)
    _, r, c = IonSim._setup_base_hamiltonian(T, 1e-6, 100, Inf, "analytic", true)
    @test length(unique([ridxs; r])) == (N^2 + N) / 2 == length(ridxs)
    @test length(unique([cidxs; c]))-1 == (N^2 - N) / 2 == length(cidxs)

    ## now shine light on both ions
    L.pointing = [(1, 1.0), (2, 1.0)]
    ridxs, cidxs = get_indices(12, [N], false)
    _, r, c = IonSim._setup_base_hamiltonian(T, 1e-6, 100, Inf, "analytic", true)
    @test length(unique([ridxs; r])) == (N^2 + N) == length(ridxs)
    @test length(unique([cidxs; c]))-1 == (N^2 - N) == length(cidxs)

    ## test for two modes
    C1 = Ca40(["S-1/2", "D-1/2"])
    L1 = Laser()
    chain1 = LinearChain(
            ions=[C1, C1], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1,2])
        )
    T1 = Trap(configuration=chain1, B=4e-4, Bhat=(x̂ + ŷ + ẑ)/√3, lasers=[L1])
    mode1 = T1.configuration.vibrational_modes.z[1]
    mode2 = T1.configuration.vibrational_modes.z[2]
    mode1.N = N-1
    mode2.N = rand(1:8)
    M = mode2.N + 1
    NM = N*M
    Efield_from_rabi_frequency!(1e6, T1, 1, 1, ("S-1/2", "D-1/2"))

    ridxs, cidxs = get_indices(12, [N, M], true)
    _, r, c = IonSim._setup_base_hamiltonian(T, 1e-6, 100, Inf, "analytic", true)
    @test length(unique([ridxs; r])) == (NM^2 + NM) == length(ridxs)
    @test length(unique([cidxs; c]))-1 == (NM^2 - NM) == length(cidxs)

    ## now set non-trivial lamb_dicke_order
    lamb_dicke_order = [rand(1:N), rand(1:M)]
    ridxs, cidxs = get_indices(12, [N, M], true, cutoff=reverse(lamb_dicke_order))
    _, r, c = IonSim._setup_base_hamiltonian(T1, 1e-6, lamb_dicke_order, Inf, "analytic", true)
    @test length(unique([ridxs; r])) == length(ridxs)
    @test length(unique([cidxs; c]))-1 == length(cidxs)

    ## when laser tuned to carrier, setting an rwa_cutoff below the vibrationl_mode frequency
    ## should have the same effect
    L.Δ = transition_frequency(T1, C1, ("S-1/2", "D-1/2"))
    _, repeated_indices, conj_repeated_indices = IonSim._setup_base_hamiltonian(T1, 1, 100, 1e5, "analytic", true)
    @test length(unique([ridxs; r])) == length(ridxs)
    @test length(unique([cidxs; c]))-1 == length(cidxs)
end

@testset "hamiltonian" begin
    C1 = Ca40(["S-1/2", "D-1/2"])
    C2 = Ca40(["S-1/2", "D-1/2"])
    L = Laser()
    L.pointing = [(1, 1.0), (2, 1.0)]
    chain = LinearChain(
                ions=[C1, C2], com_frequencies=(x=3e6,y=3e6,z=1e6),
                vibrational_modes=(;z=[1,2])
            )
    T = Trap(configuration=chain, B=4e-4, Bhat=(x̂ + ŷ + ẑ)/√3, lasers=[L])
    mode1 = T.configuration.vibrational_modes.z[1]
    mode2 = T.configuration.vibrational_modes.z[2]
    Δ = round(randn(), digits=5) * 1e5  # TODO: this begins to fail at below 1 Hz!
    L.Δ = transition_frequency(T, C1, ("S-1/2", "D-1/2")) + Δ
    ϕ = randn()
    L.ϕ = ϕ
    mode1.N = 10
    mode2.N = 9
    Ω = randn()
    Efield_from_rabi_frequency!(Ω * 1e6, T, 1, 1, ("S-1/2", "D-1/2"))

    # full hamiltonian (w conj_repeated_indices)
    timescale = 1e-6
    ion_op(t) = Ω * π * exp(-im * (2π * Δ*t*timescale + ϕ)) * C1["D-1/2"] ⊗ C1["S-1/2"]'
    η11 = get_η(mode1, L, C1)
    η12 = get_η(mode2, L, C1)
    η21 = get_η(mode1, L, C2)
    η22 = get_η(mode2, L, C2)
    mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="truncated")
    mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="truncated")
    Hp(t) = (ion_op(t) ⊗ one(C2) ⊗ mode_op1(t, η=η11) ⊗ mode_op2(t, η=η12)
             + one(C1) ⊗ ion_op(t) ⊗ mode_op1(t, η=η21) ⊗ mode_op2(t, η=η22))
    qoH(t) = Hp(t) + dagger(Hp(t))
    tp = abs(51randn())

    H = hamiltonian(T, lamb_dicke_order=101)
    H1 = hamiltonian(T, lamb_dicke_order=101, time_dependent_eta=true)
    @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4
    @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4

    mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic")
    mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="analytic")
    H1 = hamiltonian(T, lamb_dicke_order=101, displacement="analytic", time_dependent_eta=false)
    H = hamiltonian(T, lamb_dicke_order=101, displacement="analytic", time_dependent_eta=true)
    @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4
    @test H1(tp, 0).data ≈ H(tp, 0).data

    # full hamiltonian (w/o conj_repeated_indices)
    H = hamiltonian(T, lamb_dicke_order=101, rwa_cutoff=1e10)
    H1 = hamiltonian(T, lamb_dicke_order=101, time_dependent_eta=true, rwa_cutoff=1e10)
    mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="truncated")
    mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="truncated")
    @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4
    @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4

    mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic")
    mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="analytic")
    H1 = hamiltonian(T, lamb_dicke_order=101, displacement="analytic", time_dependent_eta=false, rwa_cutoff=1e10)
    H = hamiltonian(T, lamb_dicke_order=101, displacement="analytic", time_dependent_eta=true, rwa_cutoff=1e10)
    @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4
    @test H1(tp, 0).data ≈ H(tp, 0).data

    # # Lamb-Dicke
    mode_op11 = DenseOperator(mode1, diagm(0 => [1 - η11^2 * i for i in 0:mode1.N]))
    mode_op12 = DenseOperator(mode2, diagm(0 => [1 - η12^2 * i for i in 0:mode2.N]))
    mode_op21 = DenseOperator(mode1, diagm(0 => [1 - η21^2 * i for i in 0:mode1.N]))
    mode_op22 = DenseOperator(mode2, diagm(0 => [1 - η22^2 * i for i in 0:mode2.N]))
    mode_op1 = mode_op11 ⊗ mode_op12
    mode_op2 = mode_op21 ⊗ mode_op22
    Hp(t) = ion_op(t) ⊗ one(C2) ⊗ mode_op1 + one(C1) ⊗ ion_op(t) ⊗ mode_op2
    qoH(t) = Hp(t) + dagger(Hp(t))
    H = hamiltonian(T, lamb_dicke_order=0, rwa_cutoff=Inf)
    H1 = hamiltonian(T, lamb_dicke_order=0, rwa_cutoff=Inf, time_dependent_eta=true)
    H2 = hamiltonian(T, lamb_dicke_order=0, rwa_cutoff=Inf, displacement="analytic")
    H3 = hamiltonian(T, lamb_dicke_order=0, rwa_cutoff=Inf, displacement="analytic", time_dependent_eta=true)
    # only considering first order corrections to carrier (propto η^2) so this won't be perfect
    @test norm((qoH(tp) - H(tp, 0)).data) < 2
    @test norm((qoH(tp) - H1(tp, 0)).data) < 2
    @test norm((qoH(tp) - H2(tp, 0)).data) < 2
    @test norm((qoH(tp) - H3(tp, 0)).data) < 2

    # RWA
    H = hamiltonian(T, lamb_dicke_order=30, rwa_cutoff=3e5)
    H1 = hamiltonian(T, lamb_dicke_order=30, rwa_cutoff=3e5, time_dependent_eta=true)
    H2 = hamiltonian(T, lamb_dicke_order=30, rwa_cutoff=3e5, displacement="analytic")
    H3 = hamiltonian(T, lamb_dicke_order=30, rwa_cutoff=3e5, displacement="analytic", time_dependent_eta=true)
    # only considering first order corrections to carrier (propto η^2) so this won't be perfect
    @test norm((qoH(tp) - H(tp, 0)).data) < 2
    @test norm((qoH(tp) - H1(tp, 0)).data) < 2
    @test norm((qoH(tp) - H2(tp, 0)).data) < 2
    @test norm((qoH(tp) - H3(tp, 0)).data) < 2

    @test_throws AssertionError hamiltonian(T, lamb_dicke_order=[1, 2, 3], rwa_cutoff=3e5)
end
end  # end suppress


text:
using Compat.Test

@testset "Site Tag" begin
    include("Site.jl")
end

@testset "Space Tag" begin
    include("Space.jl")
end

@testset "Lattice Tag" begin
    include("Lattice.jl")
end

@testset "General Tag" begin
    include("General.jl")
end


text:
@testset "optimize angle power" begin
    λ0 = 1 #doesn't matter since everything is normalized to λ0
    k0 = 2π/λ0
    kin = 3k0
    ui = PlaneWave(0.0)
    M = 10

    sfun(N) = rounded_star(0.35λ0, 0.1λ0, 4, N)
    # N, errN = minimumN(k0, kin, sfun; tol=1e-6, N_start = 200, N_min = 100,
                # N_max = 400)
    shapes = [sfun(202)]
    P = 12#P, errP = minimumP(k0, kin, shapes[1]; P_max = 30, tol = 1e-7)

    centers = rect_grid(2, div(M,2), λ0, λ0) #2xM/2 grid with distance λ0
    ids = ones(Int, M)
    φs0 = zeros(M)
    sp = ScatteringProblem(shapes, ids, centers, φs0)

    points = [2λ0*ones(15) range(minimum(centers[:,2]), stop=maximum(centers[:,2]), length=15)]
    nhat = [ones(15) zeros(15)]
    len = norm(points[1,:] - points[end,:])
    @assert verify_min_distance(sp, points)

    fmm_options = FMMoptions(true, acc = 6, dx = 2λ0)
    groups,boxSize = divideSpace(sp.centers, fmm_options)
    P2,Q = FMMtruncation(fmm_options.acc, boxSize, k0)
    optim_options = Optim.Options(f_tol = 1e-5, iterations = 50,
                        store_trace = true, show_trace = false)

    # Allocate buffers
    opb = [OptimProblemBuffer(k0, kin, sp.centers, ui, P)]
    power_buffer = [optMatrixPwr(points, sp.centers, M, P, k0, ui, 1, nhat, len)]
    mFMM = [FMMbuildMatrices(k0, P, P2, Q, groups, sp.centers, boxSize, tri = true)]
    buf = FMMbuffer(M, P, Q, length(groups))
    scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids)
    scatteringMatrices = [scatteringMatrices]
    scatteringLU = [[lu(scatteringMatrices[1][iid]) for iid = 1:length(shapes)]]

    function fobj_test(sv::Array{PowerBuffer})
        res = -sv[1].pow
    end

    function gobj_test!(sv::Array{PowerBuffer}, opb::Array{OptimProblemBuffer})
        #calculate -(∂f/∂β)ᵀ for adjoint method
        opb[1].rhs_grad[:] = sv[1].∂pow
    end

    #initial, last must be different before first iteration
    last_φs = similar(φs0); last_φs[1] = NaN; @assert(last_φs != φs0)
    df = Optim.OnceDifferentiable(
        φs -> optimize_pwr_φ_f(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, M, buf, fmm_options, mFMM, fobj_test),
        (grad_stor, φs) -> optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, scatteringMatrices, scatteringLU, P, M, buf, fmm_options, mFMM, gobj_test!),
        φs0)

    res = Optim.optimize(df, φs0, Optim.LBFGS(), optim_options)

    power_after = calc_power([k0], [kin], P, ScatteringProblem(shapes, ids, centers,
                    res.minimizer), points, nhat, ui)[1]*len
    @test res.f_converged
    @test power_after ≈ 0.00094388385
end


text:
#############################################################################
# JuMP
# An algebraic modelling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# speed2.jl
#
# Runs some JuMP benchmarks to test for speed-related regressions.
# Test time to build model in memory.
#############################################################################

include("../solvers.jl")

using JuMP
using Random

function pMedian(solver, numFacility::Int, numCustomer::Int, numLocation::Int)
    Random.seed!(10)
    customerLocations = [rand(1:numLocation) for a in 1:numCustomer]

    buildTime = @elasped begin
        m = Model(solver = solver)

        # Facility locations
        @variable(m, 0 <= s[1:numLocation] <= 1)

        # Aux. variable: x_a,i = 1 iff the closest facility to a is at i
        @variable(m, 0 <= x[1:numLocation, 1:numCustomer] <= 1)

        # Objective: min distance
        @objective(
            m,
            Max,
            sum(
                abs(customerLocations[a] - i) * x[i, a] for
                a in 1:numCustomer, i in 1:numLocation
            )
        )

        # Constraints
        for a in 1:numCustomer
            # Subject to linking x with s
            for i in 1:numLocation
                @constraint(m, x[i, a] - s[i] <= 0)
            end
            # Subject to one of x must be 1
            @constraint(m, sum(x[i, a] for i in 1:numLocation) == 1)
        end

        # Subject to must allocate all facilities
        @constraint(m, sum(s[i] for i in 1:numLocation) == numFacility)
    end

    writeTime = @elapsed JuMP.build(m)

    return buildTime, writeTime
end

function cont5(solver, n)
    m = n
    n1 = n - 1
    m1 = m - 1
    dx = 1 / n
    T = 1.58
    dt = T / m
    h2 = dx^2
    a = 0.001
    yt = [0.5 * (1 - (j * dx)^2) for j in 0:n]

    buildTime = @elapsed begin
        mod = Model(solver = solver)
        @variable(mod, 0 <= y[0:m, 0:n] <= 1)
        @variable(mod, -1 <= u[1:m] <= 1)
        @objective(
            mod,
            Min,
            0.25 *
            dx *
            (
                (y[m, 0] - yt[1])^2 +
                2 * sum((y[m, j] - yt[j+1])^2 for j in 1:n1) +
                (y[m, n] - yt[n+1])^2
            ) + 0.25 * a * dt * (2 * sum(u[i]^2 for i in 1:m1) + u[m]^2)
        )

        # PDE
        for i in 0:m1
            for j in 1:n1
                @constraint(
                    mod,
                    h2 * (y[i+1, j] - y[i, j]) ==
                    0.5 *
                    dt *
                    (
                        y[i, j-1] - 2 * y[i, j] + y[i, j+1] + y[i+1, j-1] -
                        2 * y[i+1, j] + y[i+1, j+1]
                    )
                )
            end
        end

        # IC
        for j in 0:n
            @constraint(mod, y[0, j] == 0)
        end

        # BC
        for i in 1:m
            @constraint(mod, y[i, 2] - 4 * y[i, 1] + 3 * y[i, 0] == 0)
            @constraint(
                mod,
                y[i, n-2] - 4 * y[i, n1] + 3 * y[i, n] ==
                (2 * dx) * (u[i] - y[i, n])
            )
        end
    end

    writeTime = @elapsed JuMP.build(mod)

    return buildTime, writeTime
end

function RunTests()
    # Pmedian
    if !isempty(lp_solvers)
        pmedian_build = Float64[]
        pmedian_write = Float64[]
        for runs in 1:9
            bt, wt = pMedian(first(lp_solvers), 100, 100, 5000)
            push!(pmedian_build, bt)
            push!(pmedian_write, wt)
        end
        sort!(pmedian_build)
        sort!(pmedian_write)
        print(
            "PMEDIAN BUILD MIN=",
            minimum(pmedian_build),
            "  MED=",
            pmedian_build[5],
            "\n",
        )
        print(
            "PMEDIAN INTRN MIN=",
            minimum(pmedian_write),
            "  MED=",
            pmedian_write[5],
            "\n",
        )
    else
        @warn("PMEDIAN NOT RUN!")
    end

    # Cont5
    if !isempty(quad_solvers)
        cont5_build = Float64[]
        cont5_write = Float64[]
        for runs in 1:9
            bt, wt = cont5(first(quad_solvers), 500)
            push!(cont5_build, bt)
            push!(cont5_write, wt)
        end
        sort!(cont5_build)
        sort!(cont5_write)
        print(
            "CONT5 BUILD   MIN=",
            minimum(cont5_build),
            "  MED=",
            cont5_build[5],
            "\n",
        )
        print(
            "CONT5 INTRN   MIN=",
            minimum(cont5_write),
            "  MED=",
            cont5_write[5],
            "\n",
        )
    else
        @warn("CONT5 NOT RUN!")
    end
end

RunTests()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# General code to execute the Monte Carlo simulations for the paper.
# Copyright 2015 Gray Calhoun

outputfile = ARGS[1]
codezfile = ARGS[2]
configfile = ARGS[3]

# The next lines define the parameters of the monte carlo exercise,
# including `ncore`, the number of processors to use for the
# calculations, start the additional processors, and then define the
# parameters again on each of the individual processors. (Everything
# in `configfile` is wrapped in an `@everywhere` macro.)
include(configfile)
addprocs(ncore - 1)
include(configfile) # define variables on other processors
include(codezfile) # define estimators, etc.

# These variables are defined in `configfile`
@everywhere eachmc(x) = allmcs(x, nboot, RP, α)
mcres = mean(pmap(eachmc, fill(integer(nsims / nprocs()), nprocs())))

stats = ["naive bootstrap", "our bootstrap"]
f = open(outputfile, "w")
write(f, "T,P,Method,Size\n")
for r in 1:size(RP)[1], m in 1:length(stats)
    write(f, "$(RP[r,1] + RP[r,2]),$(RP[r,2]),$(stats[m]),$(mcres[r,m])\n")
    write(f, "$(RP[r,1] + RP[r,2]),$(RP[r,2]),nominal size,$(α)\n")
end
close(f)

max_stars_repo_name: grayclhn-econ/oosbootstrap

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: grayclhn-econ/oosbootstrap

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: grayclhn-econ/oosbootstrap

max_forks_repo_licenses: MIT Unlicense


text:
using VectorizedStatistics
using Documenter

DocMeta.setdocmeta!(VectorizedStatistics, :DocTestSetup, :(using VectorizedStatistics); recursive=true)

makedocs(;
    modules=[VectorizedStatistics],
    authors="C. Brenhin Keller",
    repo="https://github.com/brenhinkeller/VectorizedStatistics.jl/blob/{commit}{path}#{line}",
    sitename="VectorizedStatistics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://brenhinkeller.github.io/VectorizedStatistics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/brenhinkeller/VectorizedStatistics.jl",
    devbranch = "main",
)


text:
import SpecialFunctions

function flux_quad_GxI(G, nx, ny, u1, u2, rho, pr)
    # G = Array{Float64,1}(undef, 0)
    tx = ny
    ty = -nx
    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny
    beta = 0.5*rho/pr
    S1 = ut*sqrt(beta)
    S2 = un*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(Float64(pi)*beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(Float64(pi)*beta)
    A1neg = 0.5*(1.0 - SpecialFunctions.erf(S1))
    A2neg = 0.5*(1.0 - SpecialFunctions.erf(S2))
    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    G[1] = (rho*A2neg*(ut*A1neg - B1))

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg-ut*B1
    G[2] = (rho*A2neg*temp2)

    temp1 = ut*A1neg - B1
    temp2 = un*A2neg - B2
    G[3] = (rho*temp1*temp2)

    temp1 = (7.0 *pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg

    temp1 = (6.0 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1neg - B1
    temp4 = 0.5*rho*un*B2*temp1

    G[4] = (rho*A2neg*(temp2 - temp3) - temp4)
    return nothing
end

function flux_quad_GxII(G, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx
    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    # if flag == 0
    #     println(IOContext(stdout, :compact => false), rho)
    #     println(IOContext(stdout, :compact => false), pr)
    # end
    S1 = ut*sqrt(beta)
    S2 = un*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(pi*beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(pi*beta)
    A1pos = 0.5*(1.0 + SpecialFunctions.erf(S1))
    A2neg = 0.5*(1.0 - SpecialFunctions.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut^2 + un^2
    G[1] = rho * A2neg* (ut*A1pos + B1)

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    G[2] = rho*A2neg*temp2

    temp1 = ut*A1pos + B1
    temp2 = un*A2neg - B2
    G[3] = rho*temp1*temp2

    # if flag == 0
    #     println("===%%%%===")
    #     println(IOContext(stdout, :compact => false), beta)
    #     println(IOContext(stdout, :compact => false), ut)
    #     println(IOContext(stdout, :compact => false), S1)
    #     println(IOContext(stdout, :compact => false), B1)
    #     println(IOContext(stdout, :compact => false), "####")
    #     println(IOContext(stdout, :compact => false), 0.5* exp(-S1*S1))
    #     println(IOContext(stdout, :compact => false), sqrt(pi*beta))
    #     println(IOContext(stdout, :compact => false), "####")
    #     # println(IOContext(stdout, :compact => false), B2)
    #     println(IOContext(stdout, :compact => false), A1pos)
    #     # println(IOContext(stdout, :compact => false), A2neg)
    #     # println(IOContext(stdout, :compact => false), u_sqr)
    #     println(IOContext(stdout, :compact => false), temp1)
    #     println(IOContext(stdout, :compact => false), temp2)
    #     # println(IOContext(stdout, :compact => false), temp3)
    #     # println(IOContext(stdout, :compact => false), temp4)
    #     println(IOContext(stdout, :compact => false), rho)
    #     println(IOContext(stdout, :compact => false), rho*temp1*temp2)
    #     println()
    # end

    temp1 = (7 *pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos

    temp1 = (6 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1pos + B1
    temp4 = 0.5*rho*un*B2*temp1
    G[4] = rho*A2neg*(temp2 + temp3) - temp4

    return nothing
end

function flux_quad_GxIII(G, nx, ny, u1, u2, rho, pr)
    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*sqrt(beta)
    S2 = un*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(pi*beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(pi*beta)
    A1pos = 0.5*(1.0 + SpecialFunctions.erf(S1))
    A2pos = 0.5*(1.0 + SpecialFunctions.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    G[1] = (rho*A2pos*(ut*A1pos + B1))

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    G[2] = (rho*A2pos*temp2)

    temp1 = ut*A1pos + B1
    temp2 = un*A2pos + B2
    G[3] = (rho*temp1*temp2)

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos

    temp1 = (6*pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1pos + B1
    temp4 = 0.5*rho*un*B2*temp1

    G[4] = (rho*A2pos*(temp2 + temp3) + temp4)

    return nothing
end

function flux_quad_GxIV(G, nx, ny, u1, u2, rho, pr)
    # G = Array{Float64,1}(undef, 0)

    tx = ny
    ty = -nx

    ut = u1*tx + u2*ty
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*sqrt(beta)
    S2 = un*sqrt(beta)
    B1 = 0.5*exp(-S1*S1)/sqrt(pi*beta)
    B2 = 0.5*exp(-S2*S2)/sqrt(pi*beta)
    A1neg = 0.5*(1.0 - SpecialFunctions.erf(S1))
    A2pos = 0.5*(1.0 + SpecialFunctions.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    G[1]  = (rho*A2pos*(ut*A1neg - B1))

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg - ut*B1
    G[2] =  (rho*A2pos*temp2)

    temp1 = ut*A1neg - B1
    temp2 = un*A2pos + B2
    G[3] =  (rho*temp1*temp2)

    temp1 = (7.0*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg

    temp1 = (6.0 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1neg - B1
    temp4 = 0.5*rho*un*B2*temp1
    G[4] = (rho*A2pos*(temp2 - temp3) + temp4)
    return nothing
end

# print(flux_quad_GxI(1,2,3,4,5,5))

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
using Pkg

Pkg.activate()
Pkg.instantiate()
Pkg.precompile()

max_stars_repo_name: BirdClef/birdclef

max_issues_repo_name: BirdClef/birdclef

max_forks_repo_name: BirdClef/birdclef


text:
## OVERLOADING TRAIT DEFAULTS RELEVANT TO MODELS

# unexported aliases:
const Detector = Union{SupervisedDetector, UnsupervisedDetector}
const ProbabilisticDetector = Union{
    ProbabilisticSupervisedDetector,
    ProbabilisticUnsupervisedDetector
}
const DeterministicDetector = Union{
    DeterministicSupervisedDetector,
    DeterministicUnsupervisedDetector
}

const StatTraits = StatisticalTraits

StatTraits.docstring(M::Type{<:MLJType}) = name(M)

function StatTraits.docstring(M::Type{<:Model})
    return "$(name(M)) from $(package_name(M)).jl.\n" *
        "[Documentation]($(package_url(M)))."
end

StatTraits.is_supervised(::Type{<:Supervised}) = true
StatTraits.is_supervised(::Type{<:SupervisedAnnotator}) = true

StatTraits.prediction_type(::Type{<:Deterministic}) = :deterministic
StatTraits.prediction_type(::Type{<:Probabilistic}) = :probabilistic
StatTraits.prediction_type(::Type{<:Interval}) = :interval
StatTraits.prediction_type(::Type{<:ProbabilisticDetector}) = :probabilistic
StatTraits.prediction_type(::Type{<:DeterministicDetector}) = :deterministic

function StatTraits.target_scitype(::Type{<:ProbabilisticDetector})
    return AbstractVector{<:Union{Missing,OrderedFactor{2}}}
end

function StatTraits.target_scitype(::Type{<:DeterministicDetector})
    return AbstractVector{<:Union{Missing, OrderedFactor{2}}}
end

# implementation is deferred as it requires methodswith which depends upon
# InteractiveUtils which we don't want to bring here as a dependency
# (even if it's stdlib).
implemented_methods(M::Type) = implemented_methods(get_interface_mode(), M)
implemented_methods(model) = implemented_methods(typeof(model))
implemented_methods(::LightInterface, M) = errlight("implemented_methods")

for M in ABSTRACT_MODEL_SUBTYPES
    @eval(StatTraits.abstract_type(::Type{<:$M}) = $M)
end

# helper to determine the scitype of supervised models
function supervised_fit_data_scitype(M)
    I = input_scitype(M)
    T = target_scitype(M)
    ret = Tuple{I, T}
    if supports_weights(M)
        W = AbstractVector{Union{Continuous, Count}} # weight scitype
        return Union{ret, Tuple{I, T, W}}
    elseif supports_class_weights(M)
        W = AbstractDict{Finite, Union{Continuous, Count}}
        return Union{ret, Tuple{I, T, W}}
    end
    return ret
end

StatTraits.fit_data_scitype(M::Type{<:Unsupervised}) = Tuple{input_scitype(M)}
StatTraits.fit_data_scitype(::Type{<:Static}) = Tuple{}
StatTraits.fit_data_scitype(M::Type{<:Supervised}) = supervised_fit_data_scitype(M)

# In special case of `UnsupervisedAnnotator`, we allow the target
# as an optional argument to `fit` (that is ignored) so that the
# `machine` constructor will accept it as a valid argument, which
# then enables *evaluation* of the detector with labeled data:
function StatTraits.fit_data_scitype(M::Type{<:UnsupervisedAnnotator})
    return Union{Tuple{input_scitype(M)}, supervised_fit_data_scitype(M)}
end

function StatTraits.fit_data_scitype(M::Type{<:SupervisedAnnotator})
    return supervised_fit_data_scitype(M)
end

StatTraits.transform_scitype(M::Type{<:Unsupervised}) = output_scitype(M)
StatTraits.inverse_transform_scitype(M::Type{<:Unsupervised}) = input_scitype(M)

function StatTraits.predict_scitype(
    M::Type{<:Union{Deterministic, DeterministicDetector}}
)
    return target_scitype(M)
end

## FALLBACKS FOR `predict_scitype` FOR `Probabilistic` and
## `ProbabilisticDetector` MODELS

# This seems less than ideal but should reduce the number of `Unknown`
# in `prediction_type` for models which, historically, have not
# implemented the trait.

function StatTraits.predict_scitype(
    M::Type{<:Union{Probabilistic, ProbabilisticDetector}}
)
    return _density(target_scitype(M))
end

_density(::Any) = Unknown

for T in [:Continuous, :Count, :Textual]
    eval(
        quote
            function _density(::Type{AbstractArray{$T, D}}) where D
                return AbstractArray{Density{$T}, D}
            end
        end
    )
end

for T in [:Finite, :Multiclass, :OrderedFactor, :Infinite, :Continuous, :Count, :Textual]
    eval(
        quote
            function _density(::Type{AbstractArray{<:$T, D}}) where D
                return AbstractArray{Density{<:$T}, D}
            end

            _density(::Type{Table($T)}) = Table(Density{$T})
        end
    )
end

for T in [:Finite, :Multiclass, :OrderedFactor]
    eval(
        quote
            function _density(::Type{AbstractArray{<:$T{N}, D}}) where {N, D}
                return AbstractArray{Density{<:$T{N}}, D}
            end

            function _density(::Type{AbstractArray{$T{N}, D}}) where {N, D}
                return AbstractArray{Density{$T{N}}, D}
            end

            _density(::Type{Table($T{N})}) where N = Table(Density{$T{N}})
        end
    )
end


text:
mutable struct recordwriter
    io :: IO
    schtyp :: SchemaType
    sync :: NTuple{16, UInt8}
    comp
end

"""
recordwriter(schema::Tables.Schema, pth:: String; compress::Union{Nothing, Symbol}=nothing)

Opens file handle at `pth` and writes `schema` as header.
To write records as separate partitions use `writerecord(arr::recordwriter, record)`
and call close(arr::recordwriter) to close file handle.

Compression is supported via the `compress` keyword argument, and can
currently be one of `:zstd`, `:deflate`, `:bzip2`, or `:xz`.

```julia
mktempdir() do p
pth = joinpath(p,"tmp.avro")
sch = Tables.Schema([:a, :b, :c, :d],
    [Union{Missing, Int64}, Float64, Union{Missing, Int64}, Union{Int64, String}])
rw = Avro.recordwriter(sch, pth)
Avro.writerecord(rw, (a=missing, b=2., c=3, d=4))
Avro.writerecord(rw, (a=1, b=4.0, c=missing, d=5))
Avro.writerecord(rw, (a=3, b=3.7, c=3, d="5"))
close(rw)
end
```
"""
function recordwriter(schema::Tables.Schema, pth:: String; compress::Union{Nothing, Symbol}=nothing)
    io = open(pth,"w")
    schtyp = schematype(schema)
    sync = _cast(NTuple{16, UInt8}, rand(UInt128))
    comp = get(COMPRESSORS, compress, nothing)
    meta = Dict("avro.schema" => JSON3.write(schtyp))
    if comp !== nothing
        meta["avro.codec"] = String(compress)
    end
    sync = _cast(NTuple{16, UInt8}, rand(UInt128))
    buf = write((magic=MAGIC, meta=meta, sync=sync); schema=FileHeaderRecordType)
    Base.write(io, buf)
    return recordwriter(io, schtyp, sync, comp)
end

"""
recordwriter(f::Function, schema::Tables.Schema, pth:: String; compress::Union{Nothing, Symbol}=nothing)

Writes `schema` as header to `pth` and ensures file handle is closed.
To write records as separate partitions use `writerecord(arr::recordwriter, record)`.

Compression is supported via the `compress` keyword argument, and can
currently be one of `:zstd`, `:deflate`, `:bzip2`, or `:xz`.

```julia
mktempdir() do p
    pth = joinpath(p,"tmp.avro")
    sch = Tables.Schema([:a, :b, :c, :d],
        [Union{Missing, Int64}, Float64, Union{Missing, Int64}, Union{Int64, String}])
    Avro.recordwriter(sch, pth) do rw
        Avro.writerecord(rw, (a=missing, b=2., c=3, d=4))
        Avro.writerecord(rw, (a=1, b=4.0, c=missing, d=5))
        Avro.writerecord(rw, (a=3, b=3.7, c=3, d="5"))
    end
end
```
"""
function recordwriter(f::Function, schema::Tables.Schema, pth:: String; kw...)
    writer = recordwriter(schema, pth; kw...)
    try
        f(writer)
    finally
        close(writer)
    end
end

function writerecord(arr::recordwriter, record)
    nb = nbytes(arr.schtyp, record)
    bytes = Vector{UInt8}(undef, nb)
    posrecord = writevalue(Binary(), arr.schtyp, record, bytes, 1, nb, Dict())
    if arr.comp !== nothing
        bytes = transcode(arr.comp[Threads.threadid()], unsafe_wrap(Base.Array, pointer(bytes), posrecord - 1))
    end
    block = Block(1, view(bytes, 1:length(bytes)), arr.sync)
    Base.write(arr.io, write(block; schema=BlockType))
end

function close(arr::recordwriter)
    close(arr.io)
end

mutable struct recordreader
    buf:: Vector{UInt8}
    pos:: Int
    firstrecordpos:: Int
    T :: DataType
    schtype :: RecordType
    len :: Int
    io :: IO
    comp
end

"""
recordreader(pth:: String)

Open avro file and read header. Use `iterate` to access one partition at a time.
Close filehandle with `close(arr::recordreader)`

Any compression will be detected and decompressed automatically when reading.

```julia
mktempdir() do p
    pth = joinpath(p,"tmp.avro")
    sch = Tables.Schema([:a, :b], [Int, Int])
    Avro.recordwriter(sch, pth) do rw
        Avro.writerecord(rw, (a=1, b=2))
        Avro.writerecord(rw, (a=1, b=4))
    end
    rr = Avro.recordreader(pth)
    @assert NamedTuple(first(rr)) == (a=1, b=2)
    close(rr)
end
```
"""
function recordreader(pth:: String)
    fh = open(pth)
    buf = Mmap.mmap(fh)
    header, pos = readvalue(Binary(), FileHeaderRecordType, FileHeader, buf, 1, length(buf), nothing)
    sch = JSON3.read(header.meta["avro.schema"], Schema)
    comp = get(header.meta, "avro.codec", nothing)
    comp = get(DECOMPRESSORS, Symbol(comp), nothing)
    T = juliatype(sch)
    len = length(buf)
    schtype = schematype(T)
    return recordreader(buf, pos, pos, T, schtype, len, fh, comp)
end

"""
    recordreader(f::Function, pth:: String)

Open avro file and read header. Use `iterate` to access one partition at a time.

Any compression will be detected and decompressed automatically when reading.

```julia
mktempdir() do p
    pth = joinpath(p,"tmp.avro")
    sch = Tables.Schema([:a, :b], [Int, Int])
    Avro.recordwriter(sch, pth) do rw
        Avro.writerecord(rw, (a=1, b=2))
        Avro.writerecord(rw, (a=1, b=4))
    end
    Avro.recordreader(pth) do rr
        @assert NamedTuple(first(rr)) == (a=1, b=2)
    end
end
```
"""
function recordreader(f::Function, pth:: String)
    reader = recordreader(pth)
    try
        f(reader)
    finally
        close(reader)
    end
end

function iterate(r::recordreader, state=nothing)
    if isnothing(state)
        r.pos = r.firstrecordpos
    end
    try
        read(r), r
    catch e
        if isa(e, EOFError)
            return nothing
        else
            throw(e)
        end
    end
end

function read(reader::recordreader)
    reader.pos >= reader.len && throw(EOFError())
    block, pos = readvalue(Binary(), BlockType, Block, reader.buf, reader.pos, reader.len, nothing)
    reader.pos = pos
    @debug 1 "reading block: count = $(block.count), size = $(length(block.bytes))"
    bytes = block.bytes
    # uncompress
    if reader.comp !== nothing
        bytes = transcode(comp, Vector{UInt8}(bytes))
    end
    v, _ = readvalue(Binary(), reader.schtype, reader.T, bytes, 1, length(bytes), nothing)
    return v
end

function close(arr::recordreader)
    close(arr.io)
end