text:
function solve(obj::BaseSecantSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::AllowedSolution)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, AllowedSolution), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::BaseSecantSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::BaseSecantSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::AllowedSolution)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, jdouble, AllowedSolution), arg0, arg1, arg2, arg3, arg4, arg5)
end

function solve_interval(obj::BaseSecantSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "solveInterval", BracketedUnivariateSolver_Interval, (jint, UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function solve_interval(obj::BracketedUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "solveInterval", BracketedUnivariateSolver_Interval, (jint, UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2, arg3)
end


text:
using Primes

function solve(A, M)
    L = 50 # ≥ log₂ M
    if M == 1
        return L
    end

    φₘ = totient(M)
    x = solve(A, φₘ)
    l = lcm(φₘ, M)
    res = powermod(A, x, l)

    while res < L
        res += l
    end

    return res
end

function main()
    Q = parse(Int, readline())

    for _ in 1:Q
        A, M = parse.(Int, split(readline()))
        println(solve(A, M))
    end
end

main()

ext: jl


text:
function bacs ()
  bulls = cows = turns = 0
  a = {}
  srand(int(time()))
  while length(unique(a)) < 4
    push!(a,rand('0':'9'))
  end
  my_guess = unique(a)
  println("Bulls and Cows")
  while (bulls != 4)
    print("Guess? ")
	if eof(STDIN)
      s = "q"
    else
      s = chomp(readline(STDIN))
    end
    if (s == "q")
      print("My guess was "); [print(my_guess[i]) for i=1:4]
      return
    end
    guess = collect(s)
    if !(length(unique(guess)) == length(guess) == 4 && all(isdigit,guess))
      print("\nEnter four distinct digits or q to quit: ")
      continue
    end
    bulls = sum(map(==, guess, my_guess))
    cows = length(intersect(guess,my_guess)) - bulls
    println("$bulls bulls and $cows cows!")
    turns += 1
  end
  println("You guessed my number in $turns turns.")
end

max_stars_repo_path: Module 3/Chapter02/bacs.jl

max_issues_repo_path: Module 3/Chapter02/bacs.jl

max_forks_repo_path: Module 3/Chapter02/bacs.jl


text:
MOIU.@model(MIPInnerModel,
(),
(MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
(),
(),
(MOI.SingleVariable,),
(MOI.ScalarAffineFunction,),
(),
()
)

const Model{T} = MOIU.UniversalFallback{MIPInnerModel{T}}

"""
Model()
Create an empty instance of MatrixOptInterface.Model.
"""
function Model{T}() where T
    return MOIU.UniversalFallback(MIPInnerModel{T}())
end

function Base.show(io::IO, ::Model)
    print(io, "A MatrixOptInterface Model")
    return
end

abstract type AbstractLPForm{T} end
struct LPMatrixOptInterfaceForm{T} <: AbstractLPForm{T}#, V<:AbstractVector{T}, M<:AbstractMatrix{T}}
    direction::MOI.OptimizationSense
    c::Vector{T}
    A::Matrix{T}
    c_lb::Vector{T}
    c_ub::Vector{T}
    v_lb::Vector{T}
    v_ub::Vector{T}
end
struct LPStandardForm{T} <: AbstractLPForm{T}
    direction::MOI.OptimizationSense
    c::Vector{T}
    A::Matrix{T}
    b::Vector{T}
end
struct LPCannonicalForm{T} <: AbstractLPForm{T}
    direction::MOI.OptimizationSense
    c::Vector{T}
    A::Matrix{T}
    b::Vector{T}
end
struct LPSolverForm{T} <: AbstractLPForm{T}
    direction::MOI.OptimizationSense
    c::Vector{T}
    A::Matrix{T}
    b::Vector{T}
    senses::Vector{ConstraintSense}
    v_lb::Vector{T}
    v_ub::Vector{T}
end
struct MILP{T}
    lp
    variable_type
end

function change_form(::Type{F}, lp::F) where {F <: AbstractLPForm{T}} where T
    return lp
end

function change_form(::Type{LPMatrixOptInterfaceForm{T}}, lp::LPStandardForm{T}) where T
    return LPMatrixOptInterfaceForm{T}(
        lp.direction,
        lp.c,
        lp.A,
        lp.b,
        lp.b,
        fill(zero(T), length(lp.c)),
        fill(typemax(T), length(lp.c)),
    )
end
function change_form(::Type{LPMatrixOptInterfaceForm{T}}, lp::LPCannonicalForm{T}) where T
    return LPMatrixOptInterfaceForm{T}(
        lp.direction,
        lp.c,
        lp.A,
        fill(typemin(T), length(lp.b)),
        lp.b,
        fill(typemin(T), length(lp.c)),
        fill(typemax(T), length(lp.c)),
    )
end
function change_form(::Type{LPMatrixOptInterfaceForm{T}}, lp::LPSolverForm{T}) where T
    c_lb = fill(typemin(T), length(lp.b))
    c_ub = fill(typemax(T), length(lp.b))
    for i in eachindex(lp.b)
        if lp.senses[i] == LESS_THAN
            c_ub[i] = lp.b[i]
        elseif lp.senses[i] == GREATER_THAN
            c_lb[i] = lp.b[i]
        elseif lp.senses[i] == EQUAL_TO
            c_lb[i] = lp.b[i]
            c_ub[i] = lp.b[i]
        else
            error("invalid sign $(lp.senses[i])")
        end
    end
    return LPMatrixOptInterfaceForm{T}(
        lp.direction,
        lp.c,
        lp.A,
        c_lb,
        c_ub,
        lp.v_lb,
        lp.v_ub,
    )
end

function change_form(::Type{LPCannonicalForm{T}}, lp::LPMatrixOptInterfaceForm{T}) where T
    has_c_upper = Int[]
    has_c_lower = Int[]
    sizehint!(has_c_upper, length(lp.c_ub))
    sizehint!(has_c_lower, length(lp.c_ub))
    for i in eachindex(lp.c_ub)
        if lp.c_ub < Inf
            push!(has_c_upper, i)
        end
        if lp.c_lb > -Inf
            push!(has_c_lower, i)
        end
    end
    has_v_upper = Int[]
    has_v_lower = Int[]
    sizehint!(has_v_upper, length(lp.v_ub))
    sizehint!(has_v_lower, length(lp.v_ub))
    for i in eachindex(lp.v_ub)
        if lp.v_ub < Inf
            push!(has_v_upper, i)
        end
        if lp.v_lb > -Inf
            push!(has_v_lower, i)
        end
    end
    Id = Matrix{T}(I, length(lp.c), length(lp.c))
    new_A = vcat(
        lp.A[has_c_upper,:],
        -lp.A[has_c_lower,:],
        Id[has_v_upper,:],
        -Id[has_v_lower,:],
                )
    new_b = vcat(
        lp.c_ub[has_c_upper],
        -lp.c_lb[has_c_lower],
        lp.v_ub[has_v_upper],
        -lp.v_lb[has_v_lower],
    )
    return LPCannonicalForm{T}(
        lp.direction,
        lp.c,
        new_A,
        new_b
    )
end
function change_form(::Type{LPCannonicalForm{T}}, lp::F) where {F <: AbstractLPForm{T}} where T
    temp_lp = change_form(LPMatrixOptInterfaceForm{T}, lp)
    change_form(LPCannonicalForm{T}, lp)
end

function change_form(::Type{LPStandardForm{T}}, lp::LPCannonicalForm{T}) where T
    new_A = hcat(
        lp.A,
        -lp.A,
        Matrix{T}(I, length(lp.b), length(lp.b))
    )
    new_c = vcat(
        lp.c,
        -lp.c,
        fill(0.0, length(lp.b))
    )
    return LPStandardForm{T}(
        lp.direction,
        new_c,
        new_A,
        copy(lp.b)
    )
end
function change_form(::Type{LPStandardForm{T}}, lp::F) where {F <: AbstractLPForm{T}} where T
    temp_lp = change_form(LPMatrixOptInterfaceForm{T}, lp)
    new_lp = change_form(LPCannonicalForm{T}, temp_lp)
    change_form(LPStandardForm{T}, new_lp)
end

function change_form(::Type{LPSolverForm{T}}, lp::LPMatrixOptInterfaceForm{T}) where T
    new_A = copy(lp.A)
    senses = fill(LESS_THAN, length(lp.c_lb))
    new_b = fill(NaN, length(lp.c_lb))
    for i in eachindex(lp.c_lb)
        if lp.c_lb[i] == lp.c_ub[i]
            senses[i] = EQUAL_TO
            new_b[i] = lp.c_lb[i]
        elseif lp.c_lb[i] > -Inf && lp.c_ub[i] < Inf
            senses[i] = GREATER_THAN
            new_b[i] = lp.c_lb[i]
            push!(new_b, lp.c_ub[i])
            push!(sense, LESS_THAN)
            new_A = vcat(new_A, lp.A[i,:])
        elseif lp.c_lb[i] > -Inf
            senses[i] = GREATER_THAN
            new_b[i] = lp.c_lb[i]
        elseif lp.c_ub[i] < Inf
            senses[i] = LESS_THAN
            new_b[i] = lp.c_ub[i]
        end
    end
    return LPMatrixOptInterfaceForm{T}(
        lp.direction,
        lp.c,
        new_A,
        new_b,
        senses,
        lp.v_lb,
        lp.v_ub,
    )
end
function change_form(::Type{LPSolverForm{T}}, lp::F) where {F <: AbstractLPForm{T}} where T
    temp_lp = change_form(LPMatrixOptInterfaceForm{T}, lp)
    change_form(LPSolverForm{T}, temp_lp)
end

"""
Possible forms:

A) LP standard form:

opt <c, x>
s.t.
Ax == b
 x >= 0

B) LP cannonical form:

opt <c, x>
s.t.
Ax <= b

C) MBP (our standard)
opt <c, x>
s.t.
c_lb <= Ax <= c_ub
v_lb <=  x <= v_ub

D) Solver
opt <c, x>
s.t.
Ax sense c_ub
v_lb <=  x <= v_ub

Extra
vartype = {'C','I','B'}
sense = {'<','>','='}

"""
function MatrixOptimizer(raw_lp::F) where {F <: AbstractLPForm{T}} where T
    lp = change_form(LPMatrixOptInterfaceForm{T}, raw_lp)
    num_variables = length(lp.c)
    num_constraints = length(lp.c_lb)

    # use caching
    optimizer = Model{T}()

    x = MOI.add_variables(optimizer, num_variables)

    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(lp.c, x), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            objective_function)
    MOI.set(optimizer, MOI.ObjectiveSense(), lp.direction)

    # Add constraints
    for i in 1:num_constraints
        add_constraint(optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(lp.A[i,:], x), 0.0),
            lp.c_lb[i], lp.c_ub[i])
    end

    # Add bounds
    for i in 1:num_variables
        add_constraint(optimizer,
            MOI.SingleVariable(x[i]), lp.v_lb[i], lp.v_ub[i])
    end

    return optimizer
end

function add_constraint(optimizer, func, lb, ub)
    if lb == ub > -Inf
        MOI.add_constraint(optimizer, func, MOI.EqualTo(lb))
    else
        if lb > -Inf && ub < Inf
            MOI.add_constraint(optimizer, func, MOI.Interval(lb, ub))
        elseif lb > -Inf
            MOI.add_constraint(optimizer, func, MOI.GreaterThan(lb))
        elseif ub < Inf
            MOI.add_constraint(optimizer, func, MOI.LessThan(ub))
        end
    end
end


text:
"""
Store the singular value decomposition of a matrix:

    A = UΣV'

"""
struct LowRankMatrix{T} <: AbstractMatrix{T}
    U::Matrix{T}
    Σ::Matrix{T}
    V::Matrix{T}
end

size(L::LowRankMatrix) = size(L.U, 1), size(L.V, 1)
rank(L::LowRankMatrix) = size(L.Σ,1)
norm(L::LowRankMatrix) = first(L.Σ)
cond(L::LowRankMatrix) = ifelse(rank(L) < minimum(size(L)), Inf, first(L.Σ)/last(L.Σ))

istriu(L::LowRankMatrix) = false
istril(L::LowRankMatrix) = false
issymmetric(L::LowRankMatrix) = false
ishermitian(L::LowRankMatrix) = false

function getindex(L::LowRankMatrix{T},i::Integer,j::Integer) where T
    ret = zero(T)
    U, Σ, V, r = L.U, L.Σ, L.V, rank(L)
    # for k = r:-1:1
    #     ret += U[i,k]*Σ[k,k]*V[j,k]
    # end
    U[i,:]'*Σ*V[j,:]

    # ret
end


text:
# Julia wrapper for header: /usr/include/libavutil/frame.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_frame_get_best_effort_timestamp,
    av_frame_set_best_effort_timestamp,
    av_frame_get_pkt_duration,
    av_frame_set_pkt_duration,
    av_frame_get_pkt_pos,
    av_frame_set_pkt_pos,
    av_frame_get_channel_layout,
    av_frame_set_channel_layout,
    av_frame_get_channels,
    av_frame_set_channels,
    av_frame_get_sample_rate,
    av_frame_set_sample_rate,
    av_frame_get_metadata,
    av_frame_set_metadata,
    av_frame_get_decode_error_flags,
    av_frame_set_decode_error_flags,
    av_frame_get_pkt_size,
    av_frame_set_pkt_size,
    av_frame_get_qp_table,
    av_frame_set_qp_table,
    av_frame_get_colorspace,
    av_frame_set_colorspace,
    av_frame_get_color_range,
    av_frame_set_color_range,
    av_get_colorspace_name,
    av_frame_alloc,
    av_frame_free,
    av_frame_ref,
    av_frame_clone,
    av_frame_unref,
    av_frame_move_ref,
    av_frame_get_buffer,
    av_frame_is_writable,
    av_frame_make_writable,
    av_frame_copy,
    av_frame_copy_props,
    av_frame_get_plane_buffer,
    av_frame_new_side_data,
    av_frame_new_side_data_from_buf,
    av_frame_get_side_data,
    av_frame_remove_side_data,
    av_frame_apply_cropping,
    av_frame_side_data_name

function av_frame_get_best_effort_timestamp(frame)
    ccall((:av_frame_get_best_effort_timestamp, libavutil), Int64, (Ptr{AVFrame},), frame)
end

function av_frame_set_best_effort_timestamp(frame, val::Int64)
    ccall((:av_frame_set_best_effort_timestamp, libavutil), Cvoid, (Ptr{AVFrame}, Int64), frame, val)
end

function av_frame_get_pkt_duration(frame)
    ccall((:av_frame_get_pkt_duration, libavutil), Int64, (Ptr{AVFrame},), frame)
end

function av_frame_set_pkt_duration(frame, val::Int64)
    ccall((:av_frame_set_pkt_duration, libavutil), Cvoid, (Ptr{AVFrame}, Int64), frame, val)
end

function av_frame_get_pkt_pos(frame)
    ccall((:av_frame_get_pkt_pos, libavutil), Int64, (Ptr{AVFrame},), frame)
end

function av_frame_set_pkt_pos(frame, val::Int64)
    ccall((:av_frame_set_pkt_pos, libavutil), Cvoid, (Ptr{AVFrame}, Int64), frame, val)
end

function av_frame_get_channel_layout(frame)
    ccall((:av_frame_get_channel_layout, libavutil), Int64, (Ptr{AVFrame},), frame)
end

function av_frame_set_channel_layout(frame, val::Int64)
    ccall((:av_frame_set_channel_layout, libavutil), Cvoid, (Ptr{AVFrame}, Int64), frame, val)
end

function av_frame_get_channels(frame)
    ccall((:av_frame_get_channels, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_set_channels(frame, val::Integer)
    ccall((:av_frame_set_channels, libavutil), Cvoid, (Ptr{AVFrame}, Cint), frame, val)
end

function av_frame_get_sample_rate(frame)
    ccall((:av_frame_get_sample_rate, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_set_sample_rate(frame, val::Integer)
    ccall((:av_frame_set_sample_rate, libavutil), Cvoid, (Ptr{AVFrame}, Cint), frame, val)
end

function av_frame_get_metadata(frame)
    ccall((:av_frame_get_metadata, libavutil), Ptr{AVDictionary}, (Ptr{AVFrame},), frame)
end

function av_frame_set_metadata(frame, val)
    ccall((:av_frame_set_metadata, libavutil), Cvoid, (Ptr{AVFrame}, Ptr{AVDictionary}), frame, val)
end

function av_frame_get_decode_error_flags(frame)
    ccall((:av_frame_get_decode_error_flags, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_set_decode_error_flags(frame, val::Integer)
    ccall((:av_frame_set_decode_error_flags, libavutil), Cvoid, (Ptr{AVFrame}, Cint), frame, val)
end

function av_frame_get_pkt_size(frame)
    ccall((:av_frame_get_pkt_size, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_set_pkt_size(frame, val::Integer)
    ccall((:av_frame_set_pkt_size, libavutil), Cvoid, (Ptr{AVFrame}, Cint), frame, val)
end

function av_frame_get_qp_table(f, stride, _type)
    ccall((:av_frame_get_qp_table, libavutil), Ptr{Int8}, (Ptr{AVFrame}, Ptr{Cint}, Ptr{Cint}), f, stride, _type)
end

function av_frame_set_qp_table(f, buf, stride::Integer, _type::Integer)
    ccall((:av_frame_set_qp_table, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVBufferRef}, Cint, Cint), f, buf, stride, _type)
end

function av_frame_get_colorspace(frame)
    ccall((:av_frame_get_colorspace, libavutil), AVColorSpace, (Ptr{AVFrame},), frame)
end

function av_frame_set_colorspace(frame, val::AVColorSpace)
    ccall((:av_frame_set_colorspace, libavutil), Cvoid, (Ptr{AVFrame}, AVColorSpace), frame, val)
end

function av_frame_get_color_range(frame)
    ccall((:av_frame_get_color_range, libavutil), AVColorRange, (Ptr{AVFrame},), frame)
end

function av_frame_set_color_range(frame, val::AVColorRange)
    ccall((:av_frame_set_color_range, libavutil), Cvoid, (Ptr{AVFrame}, AVColorRange), frame, val)
end

function av_get_colorspace_name(val::AVColorSpace)
    ccall((:av_get_colorspace_name, libavutil), Cstring, (AVColorSpace,), val)
end

function av_frame_alloc()
    ccall((:av_frame_alloc, libavutil), Ptr{AVFrame}, ())
end

function av_frame_free(frame)
    ccall((:av_frame_free, libavutil), Cvoid, (Ptr{Ptr{AVFrame}},), frame)
end

function av_frame_ref(dst, src)
    ccall((:av_frame_ref, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVFrame}), dst, src)
end

function av_frame_clone(src)
    ccall((:av_frame_clone, libavutil), Ptr{AVFrame}, (Ptr{AVFrame},), src)
end

function av_frame_unref(frame)
    ccall((:av_frame_unref, libavutil), Cvoid, (Ptr{AVFrame},), frame)
end

function av_frame_move_ref(dst, src)
    ccall((:av_frame_move_ref, libavutil), Cvoid, (Ptr{AVFrame}, Ptr{AVFrame}), dst, src)
end

function av_frame_get_buffer(frame, align::Integer)
    ccall((:av_frame_get_buffer, libavutil), Cint, (Ptr{AVFrame}, Cint), frame, align)
end

function av_frame_is_writable(frame)
    ccall((:av_frame_is_writable, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_make_writable(frame)
    ccall((:av_frame_make_writable, libavutil), Cint, (Ptr{AVFrame},), frame)
end

function av_frame_copy(dst, src)
    ccall((:av_frame_copy, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVFrame}), dst, src)
end

function av_frame_copy_props(dst, src)
    ccall((:av_frame_copy_props, libavutil), Cint, (Ptr{AVFrame}, Ptr{AVFrame}), dst, src)
end

function av_frame_get_plane_buffer(frame, plane::Integer)
    ccall((:av_frame_get_plane_buffer, libavutil), Ptr{AVBufferRef}, (Ptr{AVFrame}, Cint), frame, plane)
end

function av_frame_new_side_data(frame, _type::AVFrameSideDataType, size::Integer)
    ccall((:av_frame_new_side_data, libavutil), Ptr{AVFrameSideData}, (Ptr{AVFrame}, AVFrameSideDataType, Cint), frame, _type, size)
end

function av_frame_new_side_data_from_buf(frame, _type::AVFrameSideDataType, buf)
    ccall((:av_frame_new_side_data_from_buf, libavutil), Ptr{AVFrameSideData}, (Ptr{AVFrame}, AVFrameSideDataType, Ptr{AVBufferRef}), frame, _type, buf)
end

function av_frame_get_side_data(frame, _type::AVFrameSideDataType)
    ccall((:av_frame_get_side_data, libavutil), Ptr{AVFrameSideData}, (Ptr{AVFrame}, AVFrameSideDataType), frame, _type)
end

function av_frame_remove_side_data(frame, _type::AVFrameSideDataType)
    ccall((:av_frame_remove_side_data, libavutil), Cvoid, (Ptr{AVFrame}, AVFrameSideDataType), frame, _type)
end

function av_frame_apply_cropping(frame, flags::Integer)
    ccall((:av_frame_apply_cropping, libavutil), Cint, (Ptr{AVFrame}, Cint), frame, flags)
end

function av_frame_side_data_name(_type::AVFrameSideDataType)
    ccall((:av_frame_side_data_name, libavutil), Cstring, (AVFrameSideDataType,), _type)
end


text:
mutable struct KernelMachineRegression{M, O, C}
    machine::M
    output::O
    cost::C
    result::Union{OptimizationResults, Nothing}
end

"""
    KernelMachineRegression(X::AbstractArray, Y::AbstractArray;
        kernel=additivegaussiankernel, cost=0, dims)

Return a `KernelMachineRegression` object. `X` and `Y` are training input and output.
`Y` can be either a vector or a matrix. `kernel` is required to have a method
`kernel(u, v, k=nothing)`,where `u` and `v` are matrices (datapoints are columns),
and `k` is the kernel evaluated in the previous hidden space.
See [`additivegaussiankernel`](@ref) for an example.
`cost` is used to regularize. Finally, `dims` is a tuple containing the dimensionality
of the spaces that compose the machine space.

!!! note
    `dims` here is different than in [`KernelMachine`](@ref), as it excludes
    the dimensionality of the output. For instance, setting `dims=(3, 3, 3)`
    on a `KernelMachineRegression` with scalar output (i.e., `y::AbstractVector`)
    would correspond to a machine with `dims=(3, 3, 3, 1)`.

Use `fit!(kr::KernelMachineRegression)` to fit the `KernelMachineRegression` object, and
`predict(kr::KernelMachineRegression, input)` to get the model prediction on a dataset.
If no `input` is given, it will default to the training data.
"""
function KernelMachineRegression(X::AbstractArray, Y::AbstractArray;
    kernel=additivegaussiankernel, cost=0, dims)

    machine = KernelMachine(kernel, permutedims(X); dims=(dims..., size(Y, 2)))
    output = permutedims(Y)
    return KernelMachineRegression(machine, output, cost, nothing)
end

Base.show(io::IO, d::KernelMachineRegression) = print(io, "KernelMachineRegression {...}")

loss(k::KernelMachineRegression) = loss(k.machine, k.output, k.cost)

function loss(km::KernelMachine, output, cost)
    pred, sq_norm = km(nothing)
    mean(abs2, pred - output) + cost * sq_norm
end

function updatemachine!(k::KernelMachineRegression, w)
    km = k.machine
    aug, cs = km.augmenter, km.cs
    n_aug, n_cs = length(aug), length(cs)
    copyto!(aug, 1, w, 1, n_aug)
    copyto!(cs, 1, w, n_aug + 1, n_cs)
    return
end

function compute_fg(k::KernelMachineRegression)
    function fg!(_, G, w)
        updatemachine!(k, w)
        if isnothing(G)
            l = loss(k)
        else
            l, back = pullback(loss, k.machine, k.output, k.cost)
            gs, = back(one(l))
            copyto!(G, gs.augmenter)
            copyto!(G, length(gs.augmenter) + 1, gs.cs)
        end
        return l
    end
end

const default_method = ConjugateGradient()

function fit!(k::KernelMachineRegression,
    method::AbstractOptimizer=default_method,
    options::Options=Options(; default_options(method)...))

    fg! = compute_fg(k)
    aug, cs = k.machine.augmenter, k.machine.cs
    init = vcat(vec(aug), vec(cs))
    res = optimize(only_fg!(fg!), init, method, options)
    updatemachine!(k, minimizer(res))
    k.result = res
    return k
end

function predict(k::KernelMachineRegression, input=nothing)
    km = k.machine
    if isnothing(input)
        res, _ = km(nothing)
    else
        nsamples = size(km.data, 2)
        res_full, _ = km(permutedims(input))
        res = res_full[:, nsamples+1:end]
    end
    return permutedims(res)
end


text:
int main() {
  int x = 5;
  double y = 42.3;
  x = (int)y;
  printInt(x);
  y = (double)x;
  prD(y);
  x = (int)3423.56;
  printInt(x);
  y = (double)4;
  prD(y);
  printInt((int)6.7);
  prD((double)126234);
  prD(453.61233254);
  return 0;
}

void prD (double d) {
  printInt((int)d);
  int j = 0;
  d = d - floor(d);
  printString(".");
  while (j < 5) {
    d = d * 10.0;
    int x = (int)d;
    printInt((int)d);
    d = d - floor(d);
    j++;
  }
  return;
}

double floor(double d) {
  return (double)((int)d);
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
# functions that populate the initial conditions
# List of functions:

@doc """
### EulerEquationMod.ICTrigonometric

Sets all components of the solution to the free stream condition according
to the angle of attack and and Mach number.

Inputs:
mesh
sbp
eqn
opts

Inputs/Outputs:
u0: vector to populate with the solution

Aliasing restrictions: none.

"""->

function ICPolynomial(mesh::AbstractMesh{Tmsh},
                      sbp::AbstractOperator{Tsbp},
                      eqn::NSData{Tsol, Tsol, 3},
                      opts,
                      u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  sigma = 0.01
  params = eqn.params
  gamma = params.euler_params.gamma
	gamma_1 = params.euler_params.gamma_1
	aoa = params.euler_params.aoa
  beta = eqn.params.sideslip_angle
	rhoInf = 1.0
  uInf = eqn.params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = eqn.params.euler_params.Ma * sin(beta) * -1
  wInf = eqn.params.euler_params.Ma * cos(beta) * sin(aoa)
	TInf = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode
  sol = zeros(Tsol, 5)

  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)
      x = coords_j[1]
      y = coords_j[2]
      z = coords_j[3]

      rho = (x-x*x)*(y-y*y)* (z - z*z)
      u   = (x-x*x)*(y-y*y)* (z - z*z)
      v   = (x-x*x)*(y-y*y)* (z - z*z)
      w   = (x-x*x)*(y-y*y)* (z - z*z)
      T   = (x-x*x)*(y-y*y)* (z - z*z)
      rho = (sigma*rho + 1.0)*rhoInf
      u   = (sigma*u + 1.0)*uInf
      v   = (sigma*v + 1.0)*vInf
      w   = (sigma*w + 1.0)*wInf
      T   = (sigma*T + 1.0)*TInf
      p   = rho*T/gamma
      E   = T/(gamma*gamma_1) + 0.5*(u*u + v*v + w*w)

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = rho*v
      u0[dofnums_j[5]] = T/(gamma*gamma_1) + 0.5*(u*u + v*v + w*w)
      u0[dofnums_j[5]] *= rho
    end
  end
end

function ICPolynomial(mesh::AbstractMesh{Tmsh},
                      sbp::AbstractOperator{Tsbp},
                      eqn::NSData{Tsol, Tsol, 2},
                      opts,
                      u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions

  sigma   = 0.01
  Tdim    = 2
  params  = eqn.params
  gamma   = params.euler_params.gamma
  gamma_1 = gamma - 1.0
  aoa     = params.euler_params.aoa
  q       = zeros(Float64, Tdim+2)
  qRef    = zeros(Float64, Tdim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode

  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      x = coords_j[1]
      y = coords_j[2]

      # calcFreeStream(eqn.params, coords_j, sol)

      q[1] = (x-x*x)*(y-y*y)
      q[2] = (x-x*x)*(y-y*y)
      q[3] = (x-x*x)*(y-y*y)
      q[4] = (x-x*x)*(y-y*y)
      q[1] = (sigma*q[1] + 1.0)*qRef[1]
      q[2] = (sigma*q[2] + 1.0)*qRef[2]
      q[3] = (sigma*q[3] + 1.0)*qRef[3]
      q[4] = (sigma*q[4] + 1.0)*qRef[4]

      u0[dofnums_j[1]] = q[1]
      u0[dofnums_j[2]] = q[1]*q[2]
      u0[dofnums_j[3]] = q[1]*q[3]
      u0[dofnums_j[4]] = q[4]/(gamma*gamma_1) + 0.5*(q[2]*q[2] + q[3]*q[3])
      u0[dofnums_j[4]] *= q[1]
    end
  end

  return nothing

end

function ICChannel(mesh::AbstractMesh{Tmsh},
                   operator::AbstractOperator{Tsbp},
                   eqn::NSData{Tsol, Tsol, 3},
                   opts,
                   u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions
  Tdim = 3
  sigma = 0.01
  pi = 3.14159265358979323846264338
  params = eqn.params
  gamma = params.euler_params.gamma
	gamma_1 = params.euler_params.gamma_1
	aoa = params.euler_params.aoa
  beta = params.sideslip_angle
	rhoInf = 1.0
  uInf = params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = params.euler_params.Ma * sin(beta) * -1
  wInf = params.euler_params.Ma * cos(beta) * sin(aoa)
	TInf = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement

  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      x = coords_j[1]
      y = coords_j[2]
      z = coords_j[3]

      rho = rhoInf * (1 + sigma*x*y*z)
      ux = sin(pi*x) + 1
      uy = sin(pi*y) + 1
      uz = sin(pi*z) + 1
      u  = (1 + sigma*ux * uy * uz )* uInf
      vx = sin(pi*x) + 1
      vy = sin(pi*y) + 1
      vz = sin(pi*z) + 1
      v  = (1 + sigma*vx * vy * vz )* vInf
      wx = sin(pi*x) + 1
      wy = sin(pi*y) + 1
      wz = sin(pi*z) + 1
      w  = (1 + sigma*wx * wy * wz) * wInf
      T  = TInf

      if !params.isViscous
        u += 0.2 * uInf
        v += 0.2 * vInf
        w += 0.2 * wInf
      end

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = rho*w
      u0[dofnums_j[5]] = T/(gamma*gamma_1) + 0.5*(u*u + v*v + w*w)
      u0[dofnums_j[5]] *= rho
    end
  end

  return nothing

end  # end function

function ICChannel(mesh::AbstractMesh{Tmsh},
                   operator::AbstractOperator{Tsbp},
                   eqn::NSData{Tsol, Tsol, 2},
                   opts,
                   u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions
  dim = 2
  params = eqn.params
  sigma = 0.1
  pi = 3.14159265358979323846264338
  gamma = eqn.params.euler_params.gamma
  gamma_1 = gamma - 1.0
  aoa = eqn.params.euler_params.aoa
  qRef = zeros(Float64, dim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode

  for i=1 : numEl
    for j=1 : nnodes
      dofnums_j = sview(mesh.dofs, :, j, i)

      x = mesh.coords[1, j, i]
      y = mesh.coords[2, j, i]

      rho = qRef[1] * (sigma*exp(sin(0.5*pi*(x+y))) +  1.0)
      ux  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[2]
      uy  = exp(y) * sin(pi*y)
      u   = ux * uy
      vx  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[3]
      vy  = exp(y) * sin(pi*y)
      v   = vx * vy
      T   = (1 + sigma*exp(0.1*x+0.1*y)) * qRef[4]
      # T   = qRef[4]

      if !eqn.params.isViscous
        u += 0.2 * qRef[2]
      end

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = T/(gamma*gamma_1) + 0.5*(u*u + v*v)
      u0[dofnums_j[4]] *= rho
    end
  end

  return nothing

end  # end function

function ICDoubleSquare(mesh::AbstractMesh{Tmsh},
                        operator::AbstractOperator{Tsbp},
                        eqn::NSData{Tsol, Tsol, 2},
                        opts,
                        u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions
  sigma = 0.01
  pi = 3.14159265358979323846264338
  gamma = 1.4
  gamma_1 = gamma - 1.0
  aoa = eqn.params.euler_params.aoa
  rhoInf = 1.0
  uInf = eqn.params.euler_params.Ma*cos(aoa)
  vInf = eqn.params.euler_params.Ma*sin(aoa)
  TInf = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode
  sol = zeros(Tsol, 4)

  si = [0.5, 1.5]
  a = [-2.375, 16.875, -45.0, 55.0, -30.0, 6.0]
  # si = [0.5, 2.5]
  # a = [-0.220703125,1.46484375,-3.515625,3.59375,-1.40625,0.1875]
  # si = [0.5, 5.0]
  # a = [-0.01610526850581724,0.10161052685058161,-0.2235431590712794,0.19102779047909357,-0.04470863181425595,0.003251536859218615]
  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      x = coords_j[1]
      y = coords_j[2]

      calcFreeStream(eqn.params, coords_j, sol)
      E = sol[4]/sol[1]
      V2 = (sol[2]*sol[2] + sol[3]*sol[3]) / (sol[1]*sol[1])

      gx = 0.0
      gy = 0.0
      if x >= si[1] && x < si[2]
        gx = a[1] + a[2]*x + a[3]*x*x + a[4]*x^3 + a[5]*x^4 + a[6]*x^5
      elseif x >= si[2]
        gx = 1.0
      elseif x <= -si[1] && x > -si[2]
        gx = a[1] - a[2]*x + a[3]*x*x - a[4]*x^3 + a[5]*x^4 - a[6]*x^5
      elseif x <= -si[2]
        gx = 1.0
      end

      if y >= si[1] && y < si[2]
        gy = a[1] + a[2]*y + a[3]*y*y + a[4]*y^3 + a[5]*y^4 + a[6]*y^5
      elseif y >= si[2]
        gy = 1.0
      elseif y <= -si[1] && y > -si[2]
        gy = a[1] - a[2]*y + a[3]*y*y - a[4]*y^3 + a[5]*y^4 - a[6]*y^5
      elseif y <= -si[2]
        gy = 1.0
      end

      rho = rhoInf
      u   = uInf * (gx + gy - gx*gy)
      v   = vInf * (gx + gy - gx*gy)
      T   = TInf

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = T/(gamma*gamma_1) + 0.5*(u*u + v*v)
      u0[dofnums_j[4]] *= rho
    end
  end

  return nothing

end  # end function

function ICTrigonometric(mesh::AbstractMesh{Tmsh},
                         operator::AbstractOperator{Tsbp},
                         eqn::NSData{Tsol, Tsol, 2},
                         opts,
                         u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions
  sigma = 0.01
  gamma = 1.4
  gamma_1 = gamma - 1.0
  aoa = eqn.params.euler_params.aoa
  rhoInf = 1.0
  uInf = eqn.params.euler_params.Ma*cos(aoa)
  vInf = eqn.params.euler_params.Ma*sin(aoa)
  TInf = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode
  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      x = coords_j[1]
      y = coords_j[2]

      x2 = 2*x*pi
      y2 = 2*y*pi
      x3 = 3*x*pi
      y3 = 3*y*pi
      x4 = 4*x*pi
      y4 = 4*y*pi
      sx2 = sin(x2)
      sy2 = sin(y2)
      sx3 = sin(x3)
      sy3 = sin(y3)
      sx4 = sin(x4)
      sy4 = sin(y4)
      cx2 = cos(x2)
      cx3 = cos(x3)
      cx4 = cos(x4)
      cy2 = cos(y2)
      cy3 = cos(y3)
      cy4 = cos(y4)
      #
      # Exact solution in form of primitive variables
      #
      rho = sx2 * sy2
      u   = sx4 * sy4
      v   = sx3 * sy3
      T   = (1.0 - cx2) * (1.0 - cy2)
      rho = (sigma*rho + 1.0)*rhoInf
      u   = (sigma*u + 1.0)*uInf
      v   = (sigma*v + 1.0)*vInf
      T   = (sigma*T + 1.0)*TInf

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = T/(gamma*gamma_1) + 0.5*(u*u + v*v)
      u0[dofnums_j[4]] = u0[dofnums_j[4]] * rho
    end
  end

  return nothing

end  # end function

function ICTrigonometric(mesh::AbstractMesh{Tmsh},
                         operator::AbstractOperator{Tsbp},
                         eqn::NSData{Tsol, Tres, 3},
                         opts,
                         u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}
  # populate u0 with initial values
  # this is a template for all other initial conditions
  sigma = 0.0001
  pi = 3.14159265358979323846264338
  gamma = 1.4
  gamma_1 = gamma - 1.0
  aoa = eqn.params.euler_params.aoa
  beta = eqn.params.sideslip_angle
  rhoInf = 1.0
  uInf = eqn.params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = eqn.params.euler_params.Ma * sin(beta) * -1
  wInf = eqn.params.euler_params.Ma * cos(beta) * sin(aoa)
  TInf = 1.0

  numEl = mesh.numEl
  nnodes = mesh.numNodesPerElement
  dofpernode = mesh.numDofPerNode
  sol = zeros(Tsol, 5)
  for i=1:numEl
    for j=1:nnodes
      xyz = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      calcFreeStream(eqn.params, xyz, sol)

      xyz1 = 1 * pi * xyz
      xyz2 = 2 * pi * xyz
      xyz4 = 4 * pi * xyz
      sin_val_1 = sin(xyz1)
      cos_val_1 = cos(xyz1)
      sin_val_2 = sin(xyz2)
      cos_val_2 = cos(xyz2)
      sin_val_4 = sin(xyz4)
      cos_val_4 = cos(xyz4)
      #
      # Exact solution in form of primitive variables
      #
      rho = sin_val_2[1] * sin_val_2[2] * sin_val_2[3]
      u   = sin_val_4[1] * sin_val_4[2] * sin_val_4[3]
      v   = sin_val_2[1] * sin_val_2[2] * sin_val_2[3]
      w   = sin_val_1[1] * sin_val_1[2] * sin_val_1[3]
      T   = (1.0 - cos_val_4[1]) * (1.0 - cos_val_4[2]) * (1.0 - cos_val_4[3])

      rho = (sigma*rho + 1.0)*rhoInf
      u = (sigma*u + 1.0) * uInf
      v = (sigma*v + 1.0) * vInf
      w = (sigma*w + 1.0) * wInf
      T = (sigma*T + 1.0) * TInf
      vel2 = u*u + v*v + w*w

      u0[dofnums_j[1]] = rho
      u0[dofnums_j[2]] = rho*u
      u0[dofnums_j[3]] = rho*v
      u0[dofnums_j[4]] = rho*w
      u0[dofnums_j[5]] = rho*(T/(gamma*gamma_1) + 0.5*vel2)
    end
  end

  return nothing
end  # end function

global const ICDict = Dict{String, Function}(
  "ICTrigonometric" => ICTrigonometric,
  "ICPolynomial" => ICPolynomial,
  "ICChannel" => ICChannel,
  "ICDoubleSquare" => ICDoubleSquare,
)


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module Theory

export gen_Clnn_theory
export Clnn2CnlmNLM, sum_m_lmeqLM
export add_local_average_effect_nowindow
export calc_NobsA, calc_CobsA, calc_CNobsA

# deprecated:
export add_local_average_effect, calc_CNAnlmNLM

using ..Modes
using QuadGK

function gen_Clnn_theory(pk, cmodes)
    # Currently this only supports the simplest model. Be ready to support more.
    Clnn = fill(0.0, getlnnsize(cmodes))
    for i=1:length(Clnn)
        l, k1, k2 = getlkk(cmodes, i)
        if k1 == k2
            Clnn[i] = pk(k1)
        end
    end
    return Clnn
end

######################### Clnn <--> CnlmNLM #######################

function Clnn2CnlmNLM(Clnn, cmodes)
    amodes = cmodes.amodes
    nlmsize = getnlmsize(amodes)
    CnlmNLM = fill(0.0, nlmsize, nlmsize)
    for i′=1:nlmsize, i=1:nlmsize
        n, l, m = getnlm(amodes, i)
        n′, l′, m′ = getnlm(amodes, i′)
        if l′==l && m′==m
            if isvalidlnn(cmodes, l, n, n′)
                idx = getidx(cmodes, l, n, n′)
                CnlmNLM[i,i′] = Clnn[idx]
            end
        end
    end
    return CnlmNLM
end

function sum_m_lmeqLM(A_nlm_NLM, modes::AnlmModes)
    clnn = fill(NaN, modes.lmax+1, modes.nmax, modes.nmax)
    for n1=1:modes.nmax, n2=1:modes.nmax
        b1 = getnlmsize(modes, n1 - 1) + 1
        b2 = getnlmsize(modes, n2 - 1) + 1
        lmax = minimum(modes.lmax_n[[n1,n2]])
        for l=0:lmax
            L = l
            b = b1 + getlmsize(l-1)
            B = b2 + getlmsize(L-1)
            cl = real(A_nlm_NLM[b,B])
            for m=1:l
                M = m
                cl += 2 * real(A_nlm_NLM[b+m,B+M])
            end
            clnn[l+1,n1,n2] = cl / (2*l + 1)
        end
    end
    return clnn
end

sum_m_lmeqLM(anlm, cmodes::ClnnModes) = clnn2clΔnn̄(sum_m_lmeqLM(anlm, cmodes.amodes), cmodes)

function clnn2clΔnn̄(clnn, cmodes::ClnnModes)
    lnnsize = getlnnsize(cmodes)
    clΔnn̄ = fill(NaN, lnnsize)
    for i=1:lnnsize
        l, n1, n2 = getlnn(cmodes, i)
        clΔnn̄[i] = (clnn[l+1,n1,n2] + clnn[l+1,n2,n1]) / 2
    end
    return clΔnn̄
end

function clΔnn̄2clnn(clΔnn̄, cmodes::ClnnModes)
    lmax = cmodes.amodes.lmax
    nmax = cmodes.amodes.nmax
    lnnsize = getlnnsize(cmodes)
    clnn = fill(NaN, lmax + 1, nmax, nmax)
    for i=1:lnnsize
        l, n1, n2 = getlnn(cmodes, i)
        clnn[l+1,n1,n2] = clΔnn̄[i]
        clnn[l+1,n2,n1] = clΔnn̄[i]
    end
    return clnn
end

#clnn2clΔnn̄(clnn, bcmodes::Modes.ClnnBinnedModesOld) = clnn2clΔnn̄(clnn, bcmodes.fake_cmodes)
#clΔnn̄2clnn(clΔnn̄, bcmodes::Modes.ClnnBinnedModesOld) = clΔnn̄2clnn(clΔnn̄, bcmodes.fake_cmodes)

###################### Local average effect #########################################

function add_local_average_effect_nowindow(CNobs, cmodes, Veff)
    # currently assumes there is no window
    nmax = cmodes.amodes.nmax_l[1]
    dn00 = √(4*π) .* [quadgk(r->r^2 * cmodes.amodes.basisfunctions(n,0,r),
                             cmodes.amodes.rmin, cmodes.amodes.rmax)[1]
                      for n=1:nmax]
    C0nn = [CNobs[getidx(cmodes, 0, n, n)] for n=1:nmax]
    C011 = C0nn[1]
    A = sum(@. dn00^2 / Veff * C0nn / C011)
    @show A
    CNlocalobs = fill(NaN, size(CNobs))
    for i=1:length(CNlocalobs)
        l, n, n′ = getlnn(cmodes, i)
        Bnn′ = 0.0
        if l == 0
            Bnn′ = C0nn[n] + Clnn[n′] - A*C011 - 6*C0nn[n]*C0nn[n′]/Veff
        end
        Clnn′ = CNobs[getidx(cmodes, 0, n, n′)]
        Clnn′ = (1 + 3 * A * C011 / Veff) * Clnn′ - dn00[n] * dn00[n′] / Veff * Bnn′
        CNlocalobs[idx] = Clnn′
    end
    return CNlocalobs
end

function calc_dn00(cmodes)
    nmax = cmodes.amodes.nmax_l[1]
    dn00 = √(4*π) .* [quadgk(r->r^2 * cmodes.amodes.basisfunctions(n,0,r),
			     cmodes.amodes.rmin, cmodes.amodes.rmax)[1]
		      for n=1:nmax]
    return dn00
end

function calc_dn00obs(dn00, Wlnn, cmodes)
    nmax = length(dn00)
    dn00obs = fill(0.0, nmax)
    for n=1:nmax
	for N=1:nmax
	    if isvalidlnn(cmodes, 0, n, N)
		idx = getidx(cmodes, 0, n, N)
		dn00obs[n] += Wlnn[idx] * dn00[N]
	    end
	end
    end
    return dn00obs
end

function calc_DWlnn(cmix, cmodes, dn00)
    lnnsize = size(cmix,2)
    δdd = fill(0.0, lnnsize)
    for i=1:lnnsize
	l, n, n′ = getlnn(cmodes, i)
	if l == 0
	    δdd[i] = dn00[n] * dn00[n′]
	end
    end
    DWlnn = cmix * δdd
    return DWlnn
end

@doc raw"""
    calc_NobsA(Nobs_th, cmix, nbar, Veff, cmodes)

Calculate the observed shot noise including the local average effect for a
constant nbar.
"""
function calc_NobsA(Nobs_th, cmix, nbar, Veff, cmodes)
    dn00 = calc_dn00(cmodes)
    dn00obs = calc_dn00obs(dn00, nbar .* Nobs_th, cmodes)
    DWlnn = calc_DWlnn(cmix, cmodes, dn00 / √Veff)

    NobsA = Nobs_th + (-2 + dn00'dn00obs / Veff) * DWlnn / nbar
    return NobsA
end

function calc_Dlnnobs(cmix, dn00, cmodes)
    lnnsize = getlnnsize(cmodes)
    nmax = cmodes.amodes.nmax
    Dlnnobs = fill(0.0, lnnsize)
    for i=1:lnnsize
	l1, n1, n2 = getlnn(cmodes, i)
	for n1′=1:nmax, n2′=1:nmax
	    if isvalidlnn(cmodes, 0, n1′, n2′)
		i′ = getidx(cmodes, 0, n1′, n2′)
		Dlnnobs[i] += (2*l1 + 1) * cmix[i,i′] * dn00[n1′] * dn00[n2′]
	    end
	end
    end
    return Dlnnobs
end

function calc_W3l1n1n2(L, N, N′, dn00dn00V, cmodes, wk_cache)
    lnnsize = getlnnsize(cmodes)
    nmax = cmodes.amodes.nmax
    ell = [L, -1, 0]
    na = [N, 0, 0]
    nb = [N′, 0, 0]
    W3nn = fill(NaN, nmax, nmax)
    W3l1n1n2 = fill(NaN, lnnsize)
    for i1=1:lnnsize
	l1, n1, n2 = getlnn(cmodes, i1)
	ell[2] = l1
	na[2] = n1
	nb[2] = n2
	for n=1:nmax
	    na[3] = n
	    for n′=1:nmax
		nb[3] = n′
		W3nn[n,n′] = window_chain(ell, na, nb, wk_cache)
		if na[1] != nb[1]
		    na[1], nb[1] = nb[1], na[1]
		    W3nn[n,n′] += window_chain(ell, na, nb, wk_cache)
		else
		    W3nn[n,n′] *= 2
		end
	    end
	end
	W3l1n1n2[i1] = tr(dn00dn00V'W3nn)
    end
    return W3l1n1n2
end

@doc raw"""
    calc_CobsA(Clnn, Nobs_th, cmix, nbar, Veff, cmodes, wk_cache=nothing;
        method=:fast, Lmax=1)

Calculate the observed power spectrum including the local average effect for a
constant nbar.

The `method` is by default an approximate formula.
"""
function calc_CobsA(Clnn, Nobs_th, cmix, nbar, Veff, cmodes, wk_cache=nothing;
        method=:fast, Lmax=1)
    dn00 = calc_dn00(cmodes)
    dn00obs = calc_dn00obs(dn00, nbar .* Nobs_th, cmodes)
    DWlnn = calc_DWlnn(cmix, cmodes, dn00 / √Veff)
    D̃lnnobs = calc_Dlnnobs(cmix, dn00 / √Veff, cmodes)
    lnnsize = getlnnsize(cmodes)

    CW3 = fill(0.0, lnnsize)
    if method != :fast
	LNNidx = Int[]
	if method == :exact_sparse
	    # only calculate ℓ < Δℓ to save time.
	    Lmax = min(cmodes.amodes.lmax, Lmax)
	    for l=0:Lmax, n=1:2, n′=n:2
		if isvalidlnn(cmodes, l, n, n′)
		    push!(LNNidx, getidx(cmodes, l, n, n′))
		end
	    end
	else
	    LNNidx = collect(1:lnnsize)
	end
	dn00dn00V = dn00 * dn00' / Veff
	for LNN=LNNidx
	    L, N, N′ = getlnn(cmodes, LNN)
	    # expensive:
	    W3l1n1n2 = calc_W3l1n1n2(L, N, N′, dn00dn00V, cmodes, wk_cache)
	    CW3[LNN] = Clnn'W3l1n1n2 / (2*L+1)
	end
    else
	# approximation
	CW3 .= 2 * D̃lnnobs'Clnn * DWlnn
    end

    ClnnobsA = (cmix * Clnn
                + D̃lnnobs'Clnn * DWlnn  # only lnn=011 really matters
		- CW3)
    return ClnnobsA
end

@doc raw"""
    calc_CNobsA(Clnn, Nobs_th, cmix, nbar, Veff, cmodes, wk_cache=nothing; kwargs...)

Calculate the observed power spectrum including shot noise with the local
average effect for a constant nbar.

`kwargs` are passed to [calc_CobsA](@ref).
"""
function calc_CNobsA(Clnn, Nobs_th, cmix, nbar, Veff, cmodes, wk_cache=nothing; kwargs...)
    NobsA = calc_NobsA(Nobs_th, cmix, nbar, Veff, cmodes)
    CobsA = calc_CNobsA(Clnn, Nobs_th, cmix, nbar, Veff, cmodes, wk_cache=nothing; kwargs...)
    return CobsA .+ NobsA
end

#################### obsolete local average effect functions ###################

function add_local_average_effect(CNlnn, cmix, Wlnn, cmodes, Veff)
    base.depwarn("'add_local_average_effect()' is deprecated. Use 'calc_NobsA()' and 'calc_CobsA()' instead.", :add_local_average_effect)
    # Note 1: We only implement the δᴷₗ₀ terms, as the others are negligible
    # when Veff is large.
    #
    # Note 2: We assume there is no bandwidth binning.

    nmax = cmodes.amodes.nmax_l[1]

    # dn00
    dn00 = √(4*π) .* [quadgk(r->r^2 * cmodes.amodes.basisfunctions(n,0,r),
                             cmodes.amodes.rmin, cmodes.amodes.rmax)[1]
                      for n=1:nmax]

    # dn00obs
    dn00obs = fill(0.0, size(dn00))
    for n=1:nmax
        for N=1:nmax
            if isvalidlnn(cmodes, 0, n, N)
                idx = getidx(cmodes, 0, n, N)
                dn00obs[n] += Wlnn[idx] * dn00[N]
            end
        end
    end

    # sum_dn00obs_C0nn
    sum_dn00obs_C0nn = fill(0.0, nmax)
    for n=1:nmax
        for n′′=1:nmax
            if isvalidlnn(cmodes, 0, n, n′′)
                i′′ = getidx(cmodes, 0, n, n′′)
                sum_dn00obs_C0nn[n] += dn00obs[n′′] * CNlnn[i′′]
            end
        end
    end

    # Dlnnobs
    Dlnnobs = fill(0.0, size(CNlnn))
    for i=1:length(Dlnnobs)
        l1, n1, n2 = getlnn(cmodes, i)
        for n1′=1:nmax, n2′=1:nmax
            if isvalidlnn(cmodes, 0, n1′, n2′)
                i′ = getidx(cmodes, 0, n1′, n2′)
                Dlnnobs[i] += (2*l1 + 1) * cmix[i,i′] * dn00[n1′] * dn00[n2′]
            end
        end
    end

    # sum_Dobs_CN
    sum_Dobs_CN = Dlnnobs' * CNlnn

    # debug
    C011 = CNlnn[getidx(cmodes, 0, 1, 1)]
    A = sum_Dobs_CN / Veff / C011
    #@show dn00 dn00obs Dlnnobs sum_dn00obs_C0nn sum_Dobs_CN Veff C011 A
    @show "first" Dlnnobs[1:10] Dlnnobs[11:20] sum_Dobs_CN Veff C011 A
    @show sum_dn00obs_C0nn[1]

    # CNAlnn
    CNAlnn = deepcopy(CNlnn)
    for n=1:nmax, n′=1:nmax
        if isvalidlnn(cmodes, 0, n, n′)
            i = getidx(cmodes, 0, n, n′)
            CNAlnn[i] -= dn00[n′] / Veff * sum_dn00obs_C0nn[n]
            @show n′,dn00[n′],n,sum_dn00obs_C0nn[n]
            CNAlnn[i] -= dn00[n] / Veff * sum_dn00obs_C0nn[n′]
            CNAlnn[i] += dn00[n] * dn00[n′] / Veff^2 * sum_Dobs_CN
        end
    end
    return CNAlnn
end

function calc_CNAnlmNLM(CNlnn::AbstractVector, wmix, cmodes, Veff)
    base.depwarn("'calc_CNAnlmNLM()' is deprecated. Use 'calc_NobsA()' and 'calc_CobsA()' instead.", :calc_CNAnlmNLM)
    # Note 1: We only implement the δᴷₗ₀ terms, as the others are negligible
    # when Veff is large.
    #
    # Note 2: We assume there is no bandwidth binning.

    amodes = cmodes.amodes
    nmax = cmodes.amodes.nmax_l[1]

    # dn00
    dn00 = √(4*π) .* [quadgk(r->r^2 * cmodes.amodes.basisfunctions(n,0,r),
                             cmodes.amodes.rmin, cmodes.amodes.rmax)[1]
                      for n=1:nmax]

    # dnlm
    dnlm = fill(0.0, getnlmsize(amodes))
    for i=1:length(dnlm)
        n, l, m = getnlm(amodes, i)
        if l == 0 && m == 0
            dnlm[i] = dn00[n]
        end
    end

    dnlmobs = wmix * dnlm

    # sum_dnlmobs_Clnn
    sum_dnlmobs_Clnn = fill(0.0im, getnlmsize(amodes))
    for i=1:length(sum_dnlmobs_Clnn)
        n, l, m = getnlm(amodes, i)
        nmaxl = cmodes.amodes.nmax_l[l+1]
        for n′′=1:nmaxl
            if isvalidlnn(cmodes, l, n, n′′)
                i′′ = getidx(amodes, n′′, l, m)
                idx′′ = getidx(cmodes, l, n, n′′)
                sum_dnlmobs_Clnn[i] += dnlmobs[i′′] * CNlnn[idx′′]
            end
        end
    end

    # Dlnnobs
    Dlnnobs = fill(0.0, getlnnsize(cmodes))
    for i=1:length(Dlnnobs)
        l1, n1, n2 = getlnn(cmodes, i)
        i1base = getidx(amodes, n1, l1, 0)
        i2base = getidx(amodes, n2, l1, 0)
        D = real(conj(dnlmobs[i1base]) * dnlmobs[i2base])
        for m=1:l1
            D += 2 * real(conj(dnlmobs[i1base+m]) * dnlmobs[i2base+m])
        end
        Dlnnobs[i] = D
    end

    # sum_Dobs_CN
    sum_Dobs_CN = Dlnnobs' * CNlnn

    # debug
    C011 = CNlnn[getidx(cmodes, 0, 1, 1)]
    A = sum_Dobs_CN / Veff / C011
    #@show dnlm dnlmobs Dlnnobs sum_dnlmobs_Clnn sum_Dobs_CN Veff C011 A
    @show "second" Dlnnobs[1:10] Dlnnobs[11:20] sum_Dobs_CN Veff C011 A
    @show sum_dnlmobs_Clnn[getidx.(amodes, 1, 0, 0)]

    # CWnlmNLM
    CWnlmNLM = fill(im*NaN, getnlmsize(amodes), getnlmsize(amodes))
    for i′=1:size(CWnlmNLM,2), i=1:size(CWnlmNLM,1)
        n, l, m = getnlm(amodes, i)
        n′, l′, m′ = getnlm(amodes, i′)
        CW = 0.0im
        if l==l′ && m==m′
            if isvalidlnn(cmodes, l, n, n′)
                idx = getidx(cmodes, l, n, n′)
                CW += CNlnn[idx]
            end
        end
        if l′==0 && m′==0
            CW -= dn00[n′] / Veff * sum_dnlmobs_Clnn[i]
            if l==l′ && m==m′
                @show n′,l′,m′,dn00[n′],n,l,m,sum_dnlmobs_Clnn[i]
            end
        end
        if l==0 && m==0
            CW -= dn00[n] / Veff * conj(sum_dnlmobs_Clnn[i′])
        end
        if l′==0 && m′==0 && l==0 && m==0
            CW += dn00[n] * dn00[n′] / Veff^2 * sum_Dobs_CN
        end
        CWnlmNLM[i,i′] = CW
    end
    return CWnlmNLM
end

end

# vim: set sw=4 et sts=4 :


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 1cae6098-7665-11eb-171c-39b12a001162
begin
	using Empirikos
	using DataFrames
	using MosekTools
	using LaTeXStrings
	using Plots
	pgfplotsx()
end

# ╔═╡ 33ff43d4-7665-11eb-38bd-b74e7d842b2e
md"""# Preliminaries

Load packages
"""

# ╔═╡ c20e1e20-7665-11eb-31d0-a9bbe2952d95
md"Load the Moving to Opportunities neighborhoods dataset (Chetty and Hendren, 2018) and remove missing entries"

# ╔═╡ 61ace78c-7665-11eb-14c0-0964539796fa
nbhood_csv = Empirikos.Neighborhoods.load_table() |> DataFrame |> dropmissing;

# ╔═╡ ea0c5878-8188-11eb-205c-c95e406496f9
first(nbhood_csv,4)

# ╔═╡ d0a1ce7a-7665-11eb-12be-a3510bc83fc1
md"Wrap the data in a type that represents the Gaussian likelihood:"

# ╔═╡ 65596c02-7665-11eb-24a6-53194a8ef2ec
Zs = NormalSample.(nbhood_csv.p25_coef, nbhood_csv.p25_se)

# ╔═╡ 75a30c76-7665-11eb-1a42-5fb73c3fd435
md" Let us define the empirical Bayes targets we want to estimate. We will estimate $\mathbb E[\mu \mid Z=z, σ = σ]$ where we vary $\sigma \in \{0.5, 1, 2\}$ and $z \in \{-3\sigma, -2.9\sigma,\dotsc, 2.9\sigma, 3\sigma\}$."

# ╔═╡ 6a5dce64-7665-11eb-0640-074cb387b172
σs = [0.5 1.0 2.0]

# ╔═╡ 6d2de886-7665-11eb-38b2-918706170dab
targets = PosteriorMean.(NormalSample.(collect(-3:0.1:3) .* σs, σs))

# ╔═╡ 70a14292-7665-11eb-37bc-078e81fad4ce
target_names = [L"E[\mu \mid Z=z, \sigma=%$σ]" for σ in σs]

# ╔═╡ 86f6f028-7665-11eb-01df-b3f0fb425fc5
md"""# Smooth prior class

Throughout this Section we first assume that the true prior $G$ is smooth and lies in the class of mixtures of

$$\left\{\mathcal{N}(\mu, 0.25^2), \mu \in \{-3,-2.9,\dotsc,3\} \right\}.$$
"""

# ╔═╡ 4738f08e-7666-11eb-1f8d-494095706ad9
gcal_smooth = MixturePriorClass(Normal.(-3:0.01:3, 0.25))

# ╔═╡ 4b7a72ee-7666-11eb-16cd-d9dcdfce85c2
md"We start by fitting the NPMLE to that prior class (this is not needed for the intervals)."

# ╔═╡ 5c57e074-7666-11eb-31f7-df6bc0db85e9
npmle_smooth_fit = fit(NPMLE(gcal_smooth, Mosek.Optimizer), Zs);

# ╔═╡ 67e70304-7666-11eb-3fac-a5b23821d5c4
plot(
    -3:0.001:3,
    x -> pdf(npmle_smooth_fit.prior, x),
    xguide = L"\mu",
    yguide = L"g(\mu)",
    title = "NPMLE prior estimate",
    label = nothing,
	size = (400,200)
)

# ╔═╡ 133d68e0-7667-11eb-02e6-5fd8f00ae23a
md" Start with confidence interval construction (compound DKW-F-Localization)"

# ╔═╡ 7641e23c-7666-11eb-33e6-554c4e2d3bc6
floc_method_smooth = FLocalizationInterval(;
    flocalization = DvoretzkyKieferWolfowitz(0.05),
    convexclass = gcal_smooth,
    solver = Mosek.Optimizer,
);

# ╔═╡ ab5dd40a-7666-11eb-39da-81b36a92c705
postmean_cis_smooth = confint.(floc_method_smooth, targets,Zs);

# ╔═╡ 7ddc4a5e-7667-11eb-0990-135e43d168b3
md" Also evaluate the plugin estimates based on the NPMLE and finally plot everything (interval bands in blue, NPMLE estimate in dashed purple and z ↦ z dotted black)."

# ╔═╡ afff51f8-7666-11eb-326c-a1d8ec9d6bb0
targets_smooth_npmle = targets.(npmle_smooth_fit.prior);

# ╔═╡ b6500f2a-7666-11eb-018b-33a641ed3650
begin
plots_smooth = plot.(Ref(-3:0.1:3),eachcol(postmean_cis_smooth), xlabel = L"z / \sigma", label = "")
plot!.(
    plots_smooth,
    Ref(-3:0.1:3.0),
    eachcol(response.(location.(targets))),
    label = "",
    seriestype = :line,
    linestyle = :dot,
    seriescolor = :gray,
)
plot!.(
    plots_smooth,
    Ref(-3:0.1:3.0),
    eachcol(targets_smooth_npmle),
    label = "",
    seriestype = :line,
    linecolor = :purple,
    linestyle = :dash,
)
plot(plots_smooth..., ylabel = target_names, layout=(1,3), size=(600,200))
end

# ╔═╡ 8850d39c-7667-11eb-27b1-73efcdae0f3e
md"""# Discrete class
Now instead we make a more standard nonparametric EB assumption, i.e. that $G$ is an arbitrary discrete distribution supported on the grid $-3,-2.9,\dotsc,2.9,3$ and repeat the same steps as we did for the smooth class above.
"""

# ╔═╡ ab9429c6-7667-11eb-3c9b-e188eddd8702
discrete_class = Empirikos.DiscretePriorClass(-3:0.01:3)

# ╔═╡ bfd3b20c-7668-11eb-26df-15c438553ebe
npmle_discrete_fit = fit(NPMLE(discrete_class, Mosek.Optimizer), Zs);

# ╔═╡ f57e856e-7669-11eb-2b81-5906651d88dc
plot(
    support(npmle_discrete_fit.prior),
    probs(npmle_discrete_fit.prior),
    seriestype = :sticks,
    xguide = L"\mu",
    yguide = L"g(\mu)",
    title = "NPMLE prior estimate",
    label = nothing,
	size = (400,200)
)

# ╔═╡ 2eeedd3c-766a-11eb-18c2-eb97bade511e
floc_method_discrete = FLocalizationInterval(;
    flocalization = DvoretzkyKieferWolfowitz(0.05),
    convexclass = discrete_class,
    solver = Mosek.Optimizer,
);

# ╔═╡ bce4edc8-7667-11eb-0d6c-5d7b99fead6c

# ╔═╡ c2316fb8-7667-11eb-00a9-1bef1ce72652
postmean_cis_discrete = confint.(floc_method_discrete, targets, Zs);

# ╔═╡ c52a446a-7667-11eb-344b-4124556ca32a
targets_discrete_npmle = targets.(npmle_discrete_fit.prior);

# ╔═╡ ca99b43a-7667-11eb-04f9-ff24c8341d77
begin
plots_discrete = plot.(Ref(-3:0.1:3),
	eachcol(postmean_cis_discrete),
	xlabel = L"z/\sigma",
	label = ""
)
plot!.(
    plots_discrete,
    Ref(-3:0.1:3.0),
    eachcol(response.(location.(targets))),
    label = "",
    seriestype = :line,
    linestyle = :dot,
    seriescolor = :gray,
)
plot!.(
    plots_discrete,
    Ref(-3:0.1:3.0),
    eachcol(targets_discrete_npmle),
    label = "",
    seriestype = :line,
    linecolor = :purple,
    linestyle = :dash,
)
plot(plots_discrete..., ylabel = target_names, layout=(1,3), size=(600,200))
end

# ╔═╡ faf674d6-816d-11eb-2e3e-335dd6bd5008
md"""# AMARI intervals

Here we also show how the AMARI intervals may be computed. We assume that the prior lies in the discrete prior class and also compare AMARI intervals to the compound-F-localization intervals (for the posterior mean of the neighborhood causal effect of Yuma).
"""

# ╔═╡ 0ddf12d8-816e-11eb-2cea-35144d2145e2
yuma_df = filter(:Name_of_CZ => isequal("Yuma"), nbhood_csv)

# ╔═╡ b533cb1e-816e-11eb-09a7-f3b3f26111fe
yuma_posterior_mean = PosteriorMean(NormalSample(yuma_df.p25_coef[1],
		                                         yuma_df.p25_se[1]))

# ╔═╡ 37dad434-8176-11eb-2266-9b8f19780ca0
md"Let us first look at the plug-in estimate from the NPMLE."

# ╔═╡ 48932a94-8174-11eb-0fe4-c15fcd5c7468
yuma_posterior_mean(npmle_discrete_fit.prior)

# ╔═╡ 43a83f5e-8176-11eb-2078-d38c094a90f8
md"Next we compute the compound-F-localization interval:"

# ╔═╡ d2c67da2-816e-11eb-1f57-0109f10e166a
ci_floc = confint(floc_method_discrete, yuma_posterior_mean, Zs)

# ╔═╡ 50cc3230-8176-11eb-39e3-33ffbb2add47
md"And finally the AMARI interval."

# ╔═╡ ef375466-816e-11eb-31e4-b7028998408d
amari = AMARI(; convexclass = discrete_class,
				flocalization = DvoretzkyKieferWolfowitz(0.01),
		        plugin_G = npmle_discrete_fit,
	            solver=Mosek.Optimizer);

# ╔═╡ 64cc087a-816f-11eb-327a-91ab08149d0e
ci_amari = confint(amari, yuma_posterior_mean, Zs)

# ╔═╡ 259e5064-8173-11eb-2a63-cbdf55cbe313
md"The AMARI interval does not contain zero (and so we can reject the null that the posterior mean for Yuma is nonnegative), while the F-localization intervals do contain zero."

# ╔═╡ Cell order:
# ╟─33ff43d4-7665-11eb-38bd-b74e7d842b2e
# ╠═1cae6098-7665-11eb-171c-39b12a001162
# ╟─c20e1e20-7665-11eb-31d0-a9bbe2952d95
# ╠═61ace78c-7665-11eb-14c0-0964539796fa
# ╠═ea0c5878-8188-11eb-205c-c95e406496f9
# ╟─d0a1ce7a-7665-11eb-12be-a3510bc83fc1
# ╠═65596c02-7665-11eb-24a6-53194a8ef2ec
# ╟─75a30c76-7665-11eb-1a42-5fb73c3fd435
# ╠═6a5dce64-7665-11eb-0640-074cb387b172
# ╠═6d2de886-7665-11eb-38b2-918706170dab
# ╠═70a14292-7665-11eb-37bc-078e81fad4ce
# ╟─86f6f028-7665-11eb-01df-b3f0fb425fc5
# ╠═4738f08e-7666-11eb-1f8d-494095706ad9
# ╟─4b7a72ee-7666-11eb-16cd-d9dcdfce85c2
# ╠═5c57e074-7666-11eb-31f7-df6bc0db85e9
# ╠═67e70304-7666-11eb-3fac-a5b23821d5c4
# ╟─133d68e0-7667-11eb-02e6-5fd8f00ae23a
# ╠═7641e23c-7666-11eb-33e6-554c4e2d3bc6
# ╠═ab5dd40a-7666-11eb-39da-81b36a92c705
# ╟─7ddc4a5e-7667-11eb-0990-135e43d168b3
# ╠═afff51f8-7666-11eb-326c-a1d8ec9d6bb0
# ╠═b6500f2a-7666-11eb-018b-33a641ed3650
# ╟─8850d39c-7667-11eb-27b1-73efcdae0f3e
# ╠═ab9429c6-7667-11eb-3c9b-e188eddd8702
# ╠═bfd3b20c-7668-11eb-26df-15c438553ebe
# ╠═f57e856e-7669-11eb-2b81-5906651d88dc
# ╠═2eeedd3c-766a-11eb-18c2-eb97bade511e
# ╟─bce4edc8-7667-11eb-0d6c-5d7b99fead6c
# ╠═c2316fb8-7667-11eb-00a9-1bef1ce72652
# ╠═c52a446a-7667-11eb-344b-4124556ca32a
# ╠═ca99b43a-7667-11eb-04f9-ff24c8341d77
# ╟─faf674d6-816d-11eb-2e3e-335dd6bd5008
# ╠═0ddf12d8-816e-11eb-2cea-35144d2145e2
# ╠═b533cb1e-816e-11eb-09a7-f3b3f26111fe
# ╟─37dad434-8176-11eb-2266-9b8f19780ca0
# ╠═48932a94-8174-11eb-0fe4-c15fcd5c7468
# ╟─43a83f5e-8176-11eb-2078-d38c094a90f8
# ╠═d2c67da2-816e-11eb-1f57-0109f10e166a
# ╟─50cc3230-8176-11eb-39e3-33ffbb2add47
# ╠═ef375466-816e-11eb-31e4-b7028998408d
# ╠═64cc087a-816f-11eb-327a-91ab08149d0e
# ╟─259e5064-8173-11eb-2a63-cbdf55cbe313


text:
## Sample Space: Omega
## ==================

# Omega mathematically is defined as a unit n dimensional Euclidean hypercube.
# It is useful in different contexts to have types for
# 1. concrete elements ω in Omega - vectors or Real values
# 2. abstract values in Omega, which are subsets of Omega.
# 3. random values of Omega.
# Because we wan't to allow users to write moels without worrying about Omega,
# we use LazyBoxes to generate needed dimensions of Omega when required.

LazyOmega() = LazyBox(Float64)
LazyOmega{T<:Real}(T2::Type{T}) = LazyBox(T2)

"Create a Euclidean unit box with dimensions `dims`"
function unit_box{T<:Real}(::Type{LazyBox{T}}, dims::Set{Int})
  abox = LazyBox(T)
  for dim in dims
    abox[dim] = Interval{T}(zero(T), one(T))
  end
  abox
end

"All kinds of Omega"
typealias Omega{T} Union{Vector{T},HyperBox{T}, LazyBox{T}, LazyRandomVector{T}, Dict{Int,T}}

"Abstract representations of sample space - euclidean box"
typealias AbstractOmega{T} Union{HyperBox{T}, LazyBox{T}}

"A concrete (i.e. not abstract) element ω in sample space Ω"
typealias ConcreteOmega{T} Union{Vector{T}, LazyRandomVector{T}}

"Build an omega of `n` dimensions"
abstract_omega(n::Int) = HyperBox([Interval(0.0,1.0) for i = 1:n])


text:
const MatOrVec{T} = Union{Matrix{T}, Vector{T}} where T

Base.@irrational sqrt2invpi 0.79788456080286535587 sqrt(big(2)/big(π))

#from here https://stackoverflow.com/questions/46671965/printing-variable-subscripts-in-julia
subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))

#count the number of type vars. there's probably a better way.
function my_unwrap_unionall(@nospecialize a)
    count = 0
    while isa(a, UnionAll)
        a = a.body
        count += 1
    end
    return count
end

@inline function to_corr(Σ)
	D = sqrt(Diagonal(Σ))
    iD = inv(D)
    R = iD * Σ * iD
    R = (R + R') / 2
end

#=
    analytical_shrinkage(X::Matrix)
Analytical nonlinear shrinkage estimator of the covariance matrix. Based on the
Matlab code from [1]. Translated to Julia and used here under MIT license by
permission from the authors.

[1] Ledoit, O., and Wolf, M. (2018), "Analytical Nonlinear Shrinkage of
Large-Dimensional Covariance Matrices", University of Zurich Econ WP 264.
https://www.econ.uzh.ch/static/workingpapers_iframe.php?id=943
=#
function analytical_shrinkage(X)
n, p = size(X)
@assert n >= 12 # important: sample size n must be >= 12
sample = Symmetric(X'*X) / n
E = eigen(sample)
lambda = E.values
u = E.vectors

# compute analytical nonlinear shrinkage kernel formula
lambda = lambda[max(1, p-n+1):p]
L = repeat(lambda, 1, min(p, n))
h = n^(-1/3) # Equation (4.9)
H = h*L'
x = (L-L') ./ H
ftilde = (3/4/sqrt(5)) * mean(max.(1 .- x.^2 ./ 5, 0) ./ H, dims=2) # Equation (4.7)
Hftemp = (-3/10/pi) * x + (3/4/sqrt(5)/pi) * (1 .- x.^2 ./ 5) .* log.(abs.((sqrt(5).-x) ./ (sqrt(5).+x))) # Equation (4.8)
Hftemp[abs.(x) .== sqrt(5)] .= (-3/10/pi) .* x[abs.(x) .== sqrt(5)]
Hftilde = mean(Hftemp./H, dims=2)
if p<=n
    dtilde = lambda ./ ((pi * (p/n) *lambda .* ftilde).^2  + (1 .- (p/n) .- pi * (p/n) * lambda .* Hftilde).^2) # Equation (4.3)
else
    Hftilde0 = (1/pi) * (3/10/h^2 + 3/4/sqrt(5)/h*(1-1/5/h^2) * log((1+sqrt(5)*h)/(1-sqrt(5)*h)))*mean(1 ./ lambda) # Equation (C.8)
    dtilde0 = 1 / (pi * (p-n) / n * Hftilde0) # Equation (C.5)
    dtilde1 = lambda ./ (pi^2*lambda.^2 .* (ftilde.^2 + Hftilde.^2)) # Eq. (C.4)
    dtilde = [dtilde0*ones(p-n,1); dtilde1]
end
return u * Diagonal(dtilde[:]) * u'
end


text:
module TestEvaluation

import Wordle
import Test: @test, @test_throws

# From https://twitter.com/jmhl/status/1483965720845512705
test_cases = (
    ("abcde", "abcde", "🟩🟩🟩🟩🟩"),
    ("abcde", "fghij", "⬜⬜⬜⬜⬜"),
    ("abcde", "axxxx", "🟩⬜⬜⬜⬜"),
    ("abcde", "bxxxx", "⬜🟨⬜⬜⬜"),
    ("abcde", "bbxxx", "⬜🟩⬜⬜⬜"),
    ("abcde", "edcba", "🟨🟨🟩🟨🟨"),
    ("reels", "rebus", "🟩🟩⬜⬜🟩"),
    ("reels", "rebus", "🟩🟩⬜⬜🟩"),
    ("reeds", "redes", "🟩🟩🟨🟨🟩"),
    ("asset", "rebus", "⬜🟨⬜🟨⬜"),
    ("igloo", "boost", "⬜⬜⬜🟨🟨"),
    ("shoot", "boost", "🟨⬜🟩🟨🟩"),
    ("ovolo", "boost", "🟨⬜🟩⬜⬜"),
)

for test_case in test_cases
    trial = Wordle.try_candidate(
        Wordle.Candidate(test_case[1]),
        Wordle.Candidate(test_case[2]),
    )
    @test trial.result === Wordle.Result(test_case[3])
end

end # module


text:
function td_0_update_V!(V, ω, s, r, sn, γ)
    ΔV = ω * (r + γ * V[sn] - V[s])
    V[s] = V[s] + ΔV
    return ΔV
end

function td_sarsa_update_Q!(Q, ω, s, a, r, sn, an, γ)
    ΔQ = ω * (r + γ * Q[sn, an] - Q[s, a])
    Q[s, a] = Q[s, a] + ΔQ
    return ΔQ
end

function td_ql_update_Q!(Q, ω, s, a, r, sn, γ)
    Qmax = max(Q[sn, :]...)
    ΔQ = ω * (r + γ * Qmax - Q[s, a])
    Q[s, a] = Q[s, a] + ΔQ
    return ΔQ
end


text:
using NCDatasets
using Formatting

fn_o = ARGS[1]

zs = collect(Float64, 0:-50:-1000)
Nz = length(zs) - 1
N_zs = length(zs)

Dataset(fn_o, "c") do ds

    defDim(ds, "Nz", Nz)
    defDim(ds, "zs", N_zs)

    defVar(ds, "zs", zs, ("zs",))

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using BSplineKit
using BandedMatrices
using LinearAlgebra
using SparseArrays
using Test

using BSplineKit.Recombinations:
    num_constraints,
    num_recombined

function test_galerkin_analytical()
    # Compare Galerkin mass matrix against analytical integrals for k = 2
    # (easy!).
    # Test with non-uniform grid (Chebyshev points).
    gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]
    x = gauss_lobatto_points(42)
    B = BSplineBasis(2, copy(x))
    t = knots(B)
    G = galerkin_matrix(B)

    let i = 8
        @test G[i, i] ≈ (t[i + 2] - t[i]) / 3
        @test G[i - 1, i] ≈ (t[i + 1] - t[i]) / 6
    end

    let B3 = BSplineBasis(3, copy(x))
        # "all bases must share the same parent B-spline basis"
        @test_throws ArgumentError galerkin_matrix((B, B3))
    end

    nothing
end

function test_galerkin_projection(B)
    # Idea: project a spline with all coefficients set to c_i = 1.
    # The resulting projection should be equal to the sum of the rows (or
    # columns?) of the Galerkin matrix.
    @testset "Projection | $deriv" for deriv ∈ (Derivative(0), Derivative(1))
        M = @inferred galerkin_matrix(B, (Derivative(0), deriv))
        S = @inferred Spline(B, ones(length(B)))
        φ = @inferred galerkin_projection(S, B, deriv)
        @test φ ≈ vec(sum(M; dims = 1))
    end

    # Length of output must be equal to length of basis
    u = Vector{Float32}(undef, length(B) + 3)
    @test_throws DimensionMismatch galerkin_projection!(sin, u, B)

    nothing
end

function test_galerkin_recombined()
    gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]
    knots_base = gauss_lobatto_points(41)

    B = BSplineBasis(6, copy(knots_base))

    @testset "Basis recombination" begin
        R0 = RecombinedBSplineBasis(Derivative(0), B)
        R1 = RecombinedBSplineBasis(Derivative(1), B)
        R2 = RecombinedBSplineBasis(Derivative(2), B)

        @testset "3D tensor" begin
            test_galerkin_tensor(R0)
            test_galerkin_tensor(R1)

            # "the first two bases must have the same lengths"
            @test_throws DimensionMismatch galerkin_tensor(
                (B, R0, R0), Derivative.((0, 0, 0)),
            )
        end

        N = length(B)
        Ñ = length(R0)
        @test Ñ == N - 2

        @inferred Galerkin._quadrature_prod(Val(order(B)))
        G = @inferred galerkin_matrix(B)
        G0 = @inferred galerkin_matrix(R0)
        G1 = @inferred galerkin_matrix(R1)
        G2 = @inferred galerkin_matrix(R2)
        @test all(isposdef.((G, G0, G1, G2)))

        let M2 = copy(G2), M = copy(G)
            @test M isa Hermitian && M2 isa Hermitian

            # "matrix will not be symmetric with deriv = $deriv"
            @test_throws ArgumentError galerkin_matrix!(M, (B, B), Derivative.((0, 1)))
            # "matrix will not be symmetric if bases are different"
            B_alt = BSplineBasis(6, knots_base .+ 0.1)
            @test_throws ArgumentError galerkin_matrix!(M, (B, B_alt))

            # "wrong dimensions of Galerkin matrix"
            @test_throws DimensionMismatch galerkin_matrix!(M2, (B, B))
        end

        # Due to the partition of unity property, the sum of all elements
        # must be the size of the domain (= ∫ 1 dx).
        L = -(boundaries(B)...) * -1
        @test sum(G) ≈ L
        @test sum(G0) < L  # Dirichlet: ϕ_1 = b_2 => PoU is not satisfied
        @test sum(G1) ≈ L  # Neumann: ϕ_1 = b_1 + b_2 => PoU is satisfied
        @test !(sum(G2) ≈ L)  # here it's not true!

        # Similarly, the sum along each column (or row) is
        #       ∫ b_i dx = (t[i + k] - t[i]) / k
        t = knots(B)
        k = order(B)
        @test all(1:N) do i
            sum(view(G, :, i)) ≈ (t[i + k] - t[i]) / k
        end

        Sym{M} = Hermitian{T,A} where {T, A<:M}

        # Some symmetric matrices
        @testset "Symmetric: $(constraints(R))" for R in (B, R0, R1, R2)
            let M = @inferred galerkin_matrix(R, Derivative.((1, 1)))
                @test M isa Sym{BandedMatrix}
            end
            let M = @inferred galerkin_matrix(R, SparseMatrixCSC{Float32})
                @test M isa Sym{SparseMatrixCSC{Float32}}
            end
        end

        # Some non-symmetric matrices
        @inferred galerkin_matrix(R1, Derivative.((0, 1)))
        @test galerkin_matrix(R1, Derivative.((0, 1))) isa BandedMatrix
        @test galerkin_matrix(R1, Derivative.((0, 1)), Matrix{Float64}) isa Matrix{Float64}
        @inferred galerkin_matrix((B, R1))
        @test galerkin_matrix((B, R1)) isa BandedMatrix

        @test size(G) == (N, N)
        @test size(G0) == size(G1) == (Ñ, Ñ)

        let r = 2:(N - 1)
            # G0 is simply the interior elements of G.
            @test G0 == @view G[r, r]

            # First row
            @test G1[1, 1] ≈ G[1, 1] + 2G[1, 2] + G[2, 2]  # (b1b1 + 2b1b2 + b2b2)
            @test @views G1[1, 2:(Ñ - 1)] ≈ G[1, 3:Ñ] .+ G[2, 3:Ñ]

            # Last row
            @test G1[Ñ, Ñ] ≈ G[N - 1, N - 1] + 2G[N, N - 1] + G[N, N]
            @test @views G1[Ñ, 2:(Ñ - 1)] ≈ G[N - 1, 3:Ñ] + G[N, 3:Ñ]

            # Interior
            @test @views G1[2:(Ñ - 1), 2:(Ñ - 1)] ≈ G[3:Ñ, 3:Ñ]
        end

        @testset "Integration by parts" begin
            # We test the relation between the H = ⟨ ϕᵢ, ϕⱼ″ ⟩ and the
            # Q = ⟨ ϕᵢ′, ϕⱼ′ ⟩ matrices, associated to the Laplacian term
            # for instance in the heat equation. These are clearly related by
            # integration by parts.
            # If the (recombined) basis satisfies homogeneous Dirichlet or
            # Neumann boundary conditions, the term ϕᵢϕⱼ′ appearing in the IBP
            # should vanish at the boundaries, in which case the matrices are
            # simply related as H = -Q.
            # We test both boundary conditions.
            @testset "BC: $D" for D in Derivative.((0, 1))
                R = RecombinedBSplineBasis(D, B)
                H = galerkin_matrix(R, Derivative.((0, 2)))
                Q = galerkin_matrix(R, Derivative.((1, 1)))
                N = length(R)
                ε = N^2 * eps(eltype(Q))  # the noise seems to scale as N^2
                @test norm(H + Q, Inf) < norm(Q, Inf) * ε

                # Similar thing for the tensor F_{ijk} = ⟨ ϕᵢ, ϕⱼ ϕₖ″ ⟩.
                @testset "3D tensor" begin
                    T″ = galerkin_tensor(R, Derivative.((0, 0, 2)))
                    T1 = galerkin_tensor(R, Derivative.((1, 0, 1)))
                    T2 = galerkin_tensor(R, Derivative.((0, 1, 1)))
                    @test T2 ≈ permutedims(T1, (2, 1, 3))
                    @test norm(T″ + T1 + T2, Inf) < norm(T1, Inf) * ε
                end
            end
        end

        bcs = (Derivative.((0, )), Derivative.((0, 1)))
        @testset "Combining bases ($bc)" for bc in bcs
            R = RecombinedBSplineBasis(bc, B)
            M = galerkin_matrix((R, B))
            @testset "3D tensor" begin
                test_galerkin_tensor(R)
            end
            @test M ≈ galerkin_matrix((B, R))'
            let (δl, δr) = num_constraints(R)
                n, m = size(M)
                @test n + δl + δr == m
                rl, rr = num_recombined(R)
                @assert rl == rr  # only supported case for now
                r = rl
                # Verify that except from the borders, the matrix is the same as
                # the "original" matrix constructed from the original B-spline
                # basis.
                M_base = galerkin_matrix(B)
                I = (δl + r + 1):(m - (δr + r))
                @test view(M, (r + 1):(n - r), I) ≈ view(M_base, I, I)
            end
        end
    end

    nothing
end

function test_galerkin_tensor(R::RecombinedBSplineBasis,
                              derivs=Derivative.((0, 1, 0)))
    B = parent(R) :: BSplineBasis

    if derivs[3] === Derivative(0)
        # Because of partition of unity, ∑_k T_{ijk} = M_{ij}.
        let M = galerkin_matrix(B, (derivs[1], derivs[2]))
            T = galerkin_tensor(B, derivs)
            Tsum = sum(T, dims=3)
            @test dropdims(Tsum, dims=3) ≈ M
        end
    end

    # The first two bases must have the same lengths
    @test_throws DimensionMismatch galerkin_tensor((B, R, R), derivs)

    A_base = galerkin_tensor((B, B, B), derivs)
    A = galerkin_tensor((B, B, R), derivs)
    @inferred galerkin_tensor((B, B, R), derivs)

    N1, N2, N3 = size(A)
    δl, δr = num_constraints(R)
    @test size(A_base) == (N1, N2, N3 + δl + δr)

    # Verify that except from the borders, the tensor is the same as the
    # "original" tensor constructed from the original B-spline basis.
    rl, rr = num_recombined(R)
    @assert rl == rr  # only supported case for now
    r = rl
    m, n = N2, N3
    I = (δl + r + 1):(m - (δr + r))
    @test view(A, I, I, (r + 1):(n - r)) == view(A_base, I, I, I)

    # "BandedTensor3D must have bandshift = (0, 0, 0)"
    let A = BandedTensor3D{Float64}(undef, size(A_base), Val(order(B) - 1);
                                    bandshift = (0, 0, 2))
        @test_throws ArgumentError galerkin_tensor!(A, (B, B, B), derivs)
    end

    let A = BandedTensor3D{Float64}(undef, size(A_base), Val(order(B)),
                                    bandshift=bandshift(A_base))
        # "BandedTensor3D must have bandwidth = $(k - 1)"
        @test_throws ArgumentError galerkin_tensor!(A, (B, B, B), derivs)
    end

    let k = order(B)
        Bp = BSplineBasis(k - 1, knots(B); augment=Val(false))
        # "wrong dimensions of Galerkin tensor"
        @test_throws DimensionMismatch galerkin_tensor!(A_base, (Bp, Bp, Bp), derivs)
    end

    nothing
end

@testset "Galerkin" begin
    test_galerkin_analytical()

    gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]

    knots_base = gauss_lobatto_points(41)
    B = BSplineBasis(5, copy(knots_base))
    test_galerkin_projection(B)

    test_galerkin_recombined()
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
# Julia wrapper for header: /usr/include/libavutil/fifo.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_fifo_alloc,
    av_fifo_alloc_array,
    av_fifo_free,
    av_fifo_freep,
    av_fifo_reset,
    av_fifo_size,
    av_fifo_space,
    av_fifo_generic_peek_at,
    av_fifo_generic_peek,
    av_fifo_generic_read,
    av_fifo_generic_write,
    av_fifo_realloc2,
    av_fifo_grow,
    av_fifo_drain,
    av_fifo_peek2

function av_fifo_alloc(size::Integer)
    ccall((:av_fifo_alloc, libavutil), Ptr{AVFifoBuffer}, (UInt32,), size)
end

function av_fifo_alloc_array(nmemb::Csize_t, size::Csize_t)
    ccall((:av_fifo_alloc_array, libavutil), Ptr{AVFifoBuffer}, (Csize_t, Csize_t), nmemb, size)
end

function av_fifo_free(f)
    ccall((:av_fifo_free, libavutil), Cvoid, (Ptr{AVFifoBuffer},), f)
end

function av_fifo_freep(f)
    ccall((:av_fifo_freep, libavutil), Cvoid, (Ptr{Ptr{AVFifoBuffer}},), f)
end

function av_fifo_reset(f)
    ccall((:av_fifo_reset, libavutil), Cvoid, (Ptr{AVFifoBuffer},), f)
end

function av_fifo_size(f)
    ccall((:av_fifo_size, libavutil), Cint, (Ptr{AVFifoBuffer},), f)
end

function av_fifo_space(f)
    ccall((:av_fifo_space, libavutil), Cint, (Ptr{AVFifoBuffer},), f)
end

function av_fifo_generic_peek_at(f, dest, offset::Integer, buf_size::Integer, func)
    ccall((:av_fifo_generic_peek_at, libavutil), Cint, (Ptr{AVFifoBuffer}, Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}), f, dest, offset, buf_size, func)
end

function av_fifo_generic_peek(f, dest, buf_size::Integer, func)
    ccall((:av_fifo_generic_peek, libavutil), Cint, (Ptr{AVFifoBuffer}, Ptr{Cvoid}, Cint, Ptr{Cvoid}), f, dest, buf_size, func)
end

function av_fifo_generic_read(f, dest, buf_size::Integer, func)
    ccall((:av_fifo_generic_read, libavutil), Cint, (Ptr{AVFifoBuffer}, Ptr{Cvoid}, Cint, Ptr{Cvoid}), f, dest, buf_size, func)
end

function av_fifo_generic_write(f, src, size::Integer, func)
    ccall((:av_fifo_generic_write, libavutil), Cint, (Ptr{AVFifoBuffer}, Ptr{Cvoid}, Cint, Ptr{Cvoid}), f, src, size, func)
end

function av_fifo_realloc2(f, size::Integer)
    ccall((:av_fifo_realloc2, libavutil), Cint, (Ptr{AVFifoBuffer}, UInt32), f, size)
end

function av_fifo_grow(f, additional_space::Integer)
    ccall((:av_fifo_grow, libavutil), Cint, (Ptr{AVFifoBuffer}, UInt32), f, additional_space)
end

function av_fifo_drain(f, size::Integer)
    ccall((:av_fifo_drain, libavutil), Cvoid, (Ptr{AVFifoBuffer}, Cint), f, size)
end

function av_fifo_peek2(f, offs::Integer)
    ccall((:av_fifo_peek2, libavutil), Ptr{UInt8}, (Ptr{AVFifoBuffer}, Cint), f, offs)
end


text:
#=##############################################################################
# DESCRIPTION
    Simulation postprocessing.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Mar 2021
  * License   : MIT
=###############################################################################

"""
Generate a grid around the vehicle where velocity and vorticity fields are
probed and outputted as VTK files.
"""
function postprocessing_fluiddomain(read_path, save_path, nums;
                                    # PROCESSING OPTIONS
                                    fdx=1/80,                # Scaling of cell length
                                    static_particles=false,  # Add static particles
                                    fsgm=0.35,               # Scaling of smoothing radii
                                    maxvpmparticles=200000,
                                    # OUTPUT OPTIONS
                                    run_name="vahana2",
                                    prompt=true, paraview=false, debug=false,
                                    verbose=true, v_lvl=0)

    fsgm          = 0.35                             # Scaling of smoothing radii

    # Maneuver parameters
    ttot          = 30.0                             # Total time to perform maneuver
    nsteps        = 4*5400                           # Time steps of complete maneuver
    dt            = ttot/nsteps
    Vref          = 0.25 * 125*0.44704               # Reference velocity magnitude

    maneuver      = generate_maneuver_vahana()

    # Functions for calculating position and orientation of vehicle over time
    V(t)          = Vref*maneuver.Vvehicle(t/ttot)
    Xt(t)         = QuadGK.quadgk(V, 0, t, rtol=1e-3)[1]
    Xn(n)         = Xt(n*dt)
    angles(n)     = maneuver.anglevehicle(n/nsteps)

    # ----------------- SETUP FLUID DOMAIN GRID ----------------------------------------
    b             = 5.86                             # Some characteristic length
    L             = 1.4*b                            # Length of fluid domain
    dx            = fdx*b                            # Cell length

    Lx, Ly, Lz    = 1.0*L, L, 0.5*L                  # Length in each direction
    dx, dy, dz    = dx, 1*dx, 0.5*dx                 # Cell length in each direction
    nx, ny, nz    = Lx/dx, Ly/dy, Lz/dz              # Number of cells in each direction

    X0            = [0.4*L, 0, 0.125*L]              # Center of fluid domain

    P_min         = -[Lx, Ly, Lz]/2
    P_max         = [Lx, Ly, Lz]/2
    NDIVS         = Int.(ceil.([nx, ny, nz]))

    # Genereate fluid domain grid
    if verbose; println("\t"^v_lvl*"Defining fluid domain..."); end;
    grid_fdom     = gt.Grid(P_min, P_max, NDIVS)

    if verbose; println("\t"^(v_lvl+1)*"Number of nodes:\t$(grid_fdom.nnodes)"); end;

    # Visualize computational and fluid domain
    str           = save_path

    if debug
        if verbose; println("\t"^v_lvl*"Saving grid..."); end;
        str *= gt.save(grid_fdom, "fdom.$(num)"; path=save_path)
    end

    # ----------------- CALCULATE FLUID DOMAIN -----------------------------------------
    settings_fname = run_name*"_pfield_settings.jld"   # Partiel field settings file
    maxparticles   = grid_fdom.nnodes + maxvpmparticles# Maximum number of particle in field

    UJ             = vpm.UJ_fmm                      # Method for UJ evaluation

    overwrite_settings = (
                            (:maxparticles, maxparticles),
                            (:sgsmodel, vpm.sgs_none),
                            (:Uinf, (args...) -> zeros(3)),
                            (:sgsscaling, vpm.sgs_scaling_none)
                         )

    # Generate particle field
    pfield = vpm.generate_particlefield(settings_fname;
                                            path=read_path,
                                            overwrite_settings=overwrite_settings)

    prevT = zeros(3)
    prevM = Float64[i==j for i in 1:3, j in 1:3]

    if verbose; println("\t"^(v_lvl)*"Iterating over steps"); end;
    for (numi, num) in enumerate(nums)

        if verbose
            println("\t"^(v_lvl+1)*"Time step: $(num)\t($(numi) out of $(length(nums)))")
        end

        # # Bring grid back to origin
        # gt.lintransform!(grid_fdom, Float64[i==j for i in 1:3, j in 1:3], -prevT; reset_fields=true)
        #
        # # Reorient grid with global axes
        # gt.lintransform!(grid_fdom, collect(prevM'), zeros(3); reset_fields=true)
        #
        # # Bring grid back to origin and reorient with global axes at the same time
        # gt.lintransform!(grid_fdom, collect(prevM'), -prevM*prevT; reset_fields=true)
        #
        # # Translate and rotate grid
        # this_X = Xn(num)
        # T = X0 + this_X                           # Translation of fluid domain
        # M = gt.rotation_matrix(-angles(num)...)   # Orientation of fluid domain
        # gt.lintransform!(grid_fdom, M, T; reset_fields=true)

        # Bring grid back to origin and reorient with global axes
        # and translate and rotate to new position; all at the same time
        T = X0 + Xn(num)                          # Translation of fluid domain
        M = gt.rotation_matrix(-angles(num)...)   # Orientation of fluid domain
        gt.lintransform!(grid_fdom, collect(prevM')*M,
                                T - collect(M')*prevM*prevT; reset_fields=true)

        prevT, prevM = T, M

        # Read particle field
        vpm.read!(pfield, run_name*"_pfield.$(num).h5";
                                                path=read_path, overwrite=true)

        # Read embedded particles
        if static_particles
            vpm.read!(pfield, run_name*"_staticpfield.$(num).h5";
                                                path=read_path, overwrite=false)
        end

        np = vpm.get_np(pfield)

        # Rescale smoothing radii
        meansigma = 0
        meannormGamma = 0
        for P in vpm.iterate(pfield)
            P.sigma[1] *= fsgm
            meansigma += P.sigma[1]
            meannormGamma += sqrt(P.Gamma[1]^2 + P.Gamma[2]^2 + P.Gamma[3]^2)
        end
        meansigma /= np
        meannormGamma /= np

        # Add grid nodes to the particle field
        Gamma = 1e-4*meannormGamma*ones(3)
        sigma = 0.5*meansigma

        for ni in 1:grid_fdom.nnodes
            vpm.add_particle(pfield, view(grid_fdom.nodes, 1:3, ni), Gamma, sigma)
        end

        # Evaluate UJ
        vpm._reset_particles(pfield)
        t = @elapsed UJ(pfield)
        if verbose; println("\t"^(v_lvl+2)*"Evaluate UJ:\t$(round(t, digits=1)) s"); end;

        t = @elapsed begin
        Us = collect(vpm.get_U(pfield, i) for i in np+1:vpm.get_np(pfield))
        Ws = collect(vpm.get_W(pfield, i) for i in np+1:vpm.get_np(pfield))

        gt.add_field(grid_fdom, "U", "vector", Us, "node")
        gt.add_field(grid_fdom, "W", "vector", Ws, "node")

        # Save fluid domain
        gt.save(grid_fdom, "fdom.$(num)"; path=save_path, time=num)
        end

        if verbose; println("\t"^(v_lvl+2)*"Save VTK:\t$(round(t, digits=1)) s"); end;

        if debug
            str *= vpm.save(pfield, "debug_pfield"; path=save_path)
        end
    end

    str *= "fdom...vtk"

    if paraview || debug
        run(`paraview --data=$(str)`)
    end

    return str
end

"""
Generate inputs for PSU-WOPWOP and run it.
"""
function postprocessing_aeracoustics()

    # Path from where to read aerodynamic solution
    # read_path       = "temps/vahana2_sim16/"
    read_path       = "/home/flowlab/ealvarez/simulations/vahana2_sim16/"
    # Path where to save PSU-WOPWOP's outputs
    save_path       = "/home/flowlab/ealvarez/simulations/vahana2_sim16-psw00/"
    # Path to PSU-WOPWOP binary (not included in FLOWUnsteady)
    wopwopbin       = "/home/flowlab/ealvarez/code/wopwop3_serial"
    # Run name (prefix of rotor files to read)
    run_name        = "vahana2"

    # ------------ PARAMETERS --------------------------------------------------
    # NOTE: Make sure that these parameters match what was used in the
    #       aerodynamic solution.

    # rotorsystems[si][ri] is the number of blades of the ri-th rotor in the si-th system
    rotorsystems    = [[5, 5], [2, 2], [2, 2], [5, 5, 5, 5]]
    nsteps          = 6480-4               # when the aero solution ended
    tend            = nsteps/21600 * 30.0  # (s) when the aero solution ended

    # Simulation parameters
    RPM             = 60                   # RPM is just a reference value to go from nrevs to simulation time
    rho             = 1.225                # (kg/m^3) air density
    speedofsound    = 342.35               # (m/s) speed of sound

    # Solver parameters
    ww_nrevs        = tend*(RPM/60)        # Number of revolutions in PSU-WOPWOP
    ww_nsteps_per_rev = nsteps/ww_nrevs    # Number of steps per revolution in PSU-WOPWOP
    num_min         = 0                    # Start reading aero files from this number
    periodic        = false                # Periodic aerodynamic solution
    const_geometry  = false                # Whether to run PSW on constant geometry from num_min
    CW              = nothing              # Clock-wise rotation of constant geometry
    highpass        = 0*0.0001             # High pass filter. Set to >0 to get rid of 0th freq in OASPL

    # Observer definition: Circular array of microphones
    sph_R           = 1.905                # (m) radial distance from rotor hub
    sph_nR          = 0
    sph_nphi        = 0
    sph_ntht        = 72                   # Number of microphones
    sph_thtmin      = 0                    # (deg) first microphone's angle
    sph_thtmax      = 360                  # (deg) last microphone's angle
    sph_phimax      = 180
    sph_rotation    = [90, 0, 0]           # Rotation of grid of microphones
    obs_name = "circle_mic_array"          # Observer file name

    # Observer definition: Single microphone
    Rmic = 0*1.905                         # (m) radial distance from rotor hub
    anglemic = 90*pi/180                   # (rad) microphone angle from plane of rotation (- below, + above)
                                           # 0deg is at the plane of rotation, 90deg is upstream
    # microphoneX = nothing                  # Comment and uncomment this to switch from array to single microphone
    microphoneX = Rmic*[-sin(anglemic), cos(anglemic), 0]

    # ------------ RUN PSU-WOPWOP ----------------------------------------------
    uns.run_noise_wopwop(read_path, run_name, RPM, rho, speedofsound, rotorsystems,
                            ww_nrevs, ww_nsteps_per_rev, save_path, wopwopbin;
                            # ---------- OBSERVERS -------------------------
                            sph_R=sph_R,
                            sph_nR=sph_nR, sph_ntht=sph_ntht,
                            sph_nphi=sph_nphi, sph_phimax=sph_phimax,
                            sph_rotation=sph_rotation,
                            sph_thtmin=sph_thtmin, sph_thtmax=sph_thtmax,
                            microphoneX=microphoneX,
                            # ---------- SIMULATION OPTIONS ----------------
                            periodic=periodic,
                            # ---------- INPUT OPTIONS ---------------------
                            num_min=num_min,
                            const_geometry=const_geometry,
                            axisrot="automatic",
                            CW=CW,
                            highpass=highpass,
                            # ---------- OUTPUT OPTIONS --------------------
                            verbose=true, v_lvl=0,
                            prompt=false, debug_paraview=false,
                            debuglvl=0,                     # WW debug level
                            observerf_name="observergrid",  # .xyz file with observer grid
                            case_name="runcase",            # Name of case to create and run
                            );
end


text:
module WeberCedrus

using PyCall
export Cedrus

using Weber
import Weber: keycode, iskeydown, iskeyup, addtrial, addpractice, poll_events
import Base: hash, isless, ==, show

mutable struct Cedrus <: Weber.Extension
  devices::PyObject
  trial_start::Float64
end

function InitExtension()
  pyxid = pyimport_conda("pyxid","pyxid","haberdashPI")
  pyxid[:use_response_pad_timer] = true
  cedrus = Cedrus(pyxid[:get_xid_devices](),0.0)
  for dev in cedrus.devices
    if dev[:is_response_device]()
      dev[:reset_base_timer]()
    end
  end
  cedrus
end

@Weber.event struct CedrusDownEvent <: Weber.ExpEvent
  code::Int
  port::Int
  time::Float64
end

@Weber.event struct CedrusUpEvent <: Weber.ExpEvent
  code::Int
  port::Int
  time::Float64
end

struct CedrusKey <: Weber.Key
  code::Int
end

# make sure the cedrus keys have a well defined ordering
hash(x::CedrusKey,h::UInt) = hash(CedrusKey,hash(x.code,h))
==(x::CedrusKey,y::CedrusKey) = x.code == y.code
isless(x::CedrusKey,y::CedrusKey) = isless(x.code,y.code)

# make sure cedrus keys are displayed in a easily readable form
function show(io::IO,x::CedrusKey)
  if 0 <= x.code <= 19
    write(io,"key\":cedrus$(x.code):\"")
  else
    write(io,"Weber.CedrusKey($(x.code))")
  end
end

merge!(Weber.str_to_code,Dict(":cedrus$i:" => CedrusKey(i) for i in 0:19))

keycode(e::CedrusDownEvent) = CedrusKey(e.code)
keycode(e::CedrusUpEvent) = CedrusKey(e.code)

iskeydown(event::CedrusDownEvent) = true
iskeydown(key::CedrusKey) = e -> iskeydown(e,key::CedrusKey)
iskeydown(event::CedrusDownEvent,key::CedrusKey) = event.code == key.code

iskeyup(event::CedrusUpEvent) = true
iskeyup(key::CedrusKey) = e -> iskeydown(e,key)
iskeyup(event::CedrusUpEvent,key::CedrusKey) = event.code == key.code

time(e::CedrusUpEvent) = e.time
time(e::CedrusDownEvent) = e.time

function reset_response(cedrus::Cedrus)
  old_len = length(cedrus.devices)
  cedrus.devices = pyxid[:get_xid_devices]()
  if old_len != length(cedrus.devices)
    warn("The number of available Cedrus devices changed from $old_len to "*
         "$(length(cedrus.devices)) in the middle of an experiment!")
  end

  for dev in cedrus.devices
    if dev[:is_response_device]()
      dev[:reset_rt_timer]()
    end
  end
  cedrus.trial_start = Weber.tick()
end

function addtrial(e::ExtendedExperiment{Cedrus},moments...)
  addtrial(next(e),moment(reset_response,extension(e)),moments...)
end

function addpractice(e::ExtendedExperiment{Cedrus},moments...)
  addpractice(next(e),moment(reset_response,extension(e)),moments...)
end

function poll_events(callback::Function,exp::ExtendedExperiment{Cedrus},time::Float64)
  poll_events(callback,next(exp),time)
  for dev in extension(exp).devices
    if dev[:is_response_device]()
      dev[:poll_for_response]()
      while dev[:response_queue_size]() > 0
        resp = dev[:get_next_response]()
        if resp["pressed"]
          callback(exp,
                   CedrusDownEvent(resp["key"],resp["port"],
                                   resp["time"] + extension(exp).trial_start))
        else
          callback(exp,
                   CedrusUpEvent(resp["key"],resp["port"],
                                 resp["time"] + extension(exp).trial_start))
        end
      end
    end
  end
end

end


text:
using Test
using IntervalArithmetic

@testset "construction" begin
@testset "Interval() strict w.r.t. invalid inputs" begin
    @test_throws AssertionError Interval(3.0, 2.0)
    @test_throws AssertionError Interval(-1.0, -2.0)
    @test_throws AssertionError Interval(Inf, Inf)
    @test_throws AssertionError Interval(-Inf, -Inf)
    @test_throws AssertionError Interval(NaN, 1.0)
end

@testset "interval() flexible w.r.t. invalid inputs" begin
    @test interval(3, 2) |> isempty
    @test interval(NaN) |> isnai
    @test interval(Inf, Inf) |> isempty
    @test interval(-Inf, -Inf) |> isempty
    @test interval(1.0, 2.0) isa Interval
end
end # construction

@testset "isempty() distinguishes empty and non-empty" begin
    @test emptyinterval() |> isempty
    @test interval(1, 2) |> !isempty
end

@testset "isnai() distinguishes intervals and non-intervals" begin
    @test nai() |> isnai
    @test interval(1, 2) |> !isnai
    @test emptyinterval() |> !isnai
end

@testset "unit ball" begin
    a = unit_ball(Float64)
    @test diam(a) == 2
    @test radius(a) == 1
    @test mid(a) == 0
    @test zero_in(a)
end

@testset "nonnegatives" begin
    a = nonnegatives(Float64)
    @test diam(a) |> isinf
    @test inf(a) == 0
    @test sup(a) |> isinf
    @test isunbounded(a)
    @test zero_in(a)
    @test intersect(a, -a) == interval(0)
end

@testset "entire real line" begin
    a = entireinterval(Float64)
    @test diam(a) |> isinf
    @test inf(a) |> isinf
    @test sup(a) |> isinf
    @test isunbounded(a)
    @test zero_in(a)
    @test intersect(a, -a) == a
end


text:
# The following block allows reading Union types written prior to v0.2
const LEGACY_H5TYPE_UNION = VariableLengthDatatype(H5TYPE_DATATYPE)

function jlconvert(::ReadRepresentation{Union, Vlen{DataTypeODR()}}, f::JLDFile,
                   ptr::Ptr, header_offset::RelOffset)
    v = Union{jlconvert(ReadRepresentation{DataType,Vlen{DataTypeODR()}}(), f, ptr, NULL_REFERENCE)...}
    track_weakref!(f, header_offset, v)
    v
end

constructrr(::JLDFile, ::Type{T}, dt::VariableLengthDatatype, ::Vector{ReadAttribute}) where {T<:Union} =
    dt == LEGACY_H5TYPE_UNION ? (ReadRepresentation{Union,Vlen{DataTypeODR()}}(), true) :
                         throw(UnsupportedFeatureException())

# The following definition is needed to correctly load Strings written
# with JLD2 with versions v0.1.12 - v0.3.1
function read_array(f::JLDFile, dataspace::ReadDataspace,
                    rr::FixedLengthString{String}, data_length::Int,
                    filter_id::UInt16, header_offset::RelOffset,
                    attributes::Union{Vector{ReadAttribute},Nothing})
    rrv = ReadRepresentation{UInt8,odr(UInt8)}()
    v = read_array(f, dataspace, rrv, data_length, filter_id, NULL_REFERENCE, attributes)
    String(v)
end

@static if VERSION < v"1.7.0-A"
    # Location of `mutable` flag is moved from datatype to typename in julia v1.7
    # Switch to using accessor function added in v1.7

    # Borrowed from julia base
    function ismutabletype(@nospecialize(t::Type))
        t = Base.unwrap_unionall(t)
        # TODO: what to do for `Union`?
        return isa(t, DataType) && t.mutable
    end
end


text:
##############################################################################################
module TestTimeSeriesGeneration
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.GenerateTimeSeries
using JLD, Random

##############################################################################################
##############################################################################################
# Test generation and loading of the L96 model time series in default localtion

function testGenL96()
    try
        args = (0, 40, 0.05, 5000, 1500, 0.00, 8.0)
        L96_time_series(args)
        true
    catch
        false
    end
end

function testLoadL96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        load(path * "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0_tanl_0.05" *
             "_nanl_05000_spin_1500_h_0.050.jld2")
        true
    catch
        false
    end
end

##############################################################################################
# Test generation and loading of the IEEE39 bus model time series in the default location

function testGenIEEE39bus()
    try
        args = (0, 0.01, 5000, 1500, 0.0)
        IEEE39bus_time_series(args)
        true
    catch
        false
    end
end

function testLoadIEEE39bus()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        load(path * "IEEE39bus_time_series_seed_0000_diff_0.000_tanl_0.01" *
             "_nanl_05000_spin_1500_h_0.010.jld2")
        true
    catch
        false
    end
end

#######################################################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks