text:
using LoopVectorization
using TriangularSolve: ldiv!
using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, checknonsingular, BLAS, LinearAlgebra, Adjoint, Transpose
using StrideArraysCore
using Polyester: @batch

# 1.7 compat
normalize_pivot(t::Val{T}) where T = t
to_stdlib_pivot(t) = t
if VERSION >= v"1.7.0-DEV.1188"
    normalize_pivot(::LinearAlgebra.RowMaximum) = Val(true)
    normalize_pivot(::LinearAlgebra.NoPivot) = Val(false)
    to_stdlib_pivot(::Val{true}) = LinearAlgebra.RowMaximum()
    to_stdlib_pivot(::Val{false}) = LinearAlgebra.NoPivot()
end

function lu(A::AbstractMatrix, pivot = Val(true), thread = Val(true); kwargs...)
    return lu!(copy(A), normalize_pivot(pivot), thread; kwargs...)
end

function lu!(A, pivot = Val(true), thread = Val(true); check=true, kwargs...)
    m, n  = size(A)
    minmn = min(m, n)
    F = if minmn < 10 # avx introduces small performance degradation
        LinearAlgebra.generic_lufact!(A, to_stdlib_pivot(pivot); check=check)
    else
        lu!(A, Vector{BlasInt}(undef, minmn), normalize_pivot(pivot), thread; check=check, kwargs...)
    end
    return F
end

for (f, T) in [(:adjoint, :Adjoint), (:transpose, :Transpose)], lu in (:lu, :lu!)
  @eval $lu(A::$T, args...; kwargs...) = $f($lu(parent(A), args...; kwargs...))
end

const RECURSION_THRESHOLD = Ref(-1)

# AVX512 needs a smaller recursion limit
function pick_threshold()
    RECURSION_THRESHOLD[] >= 0 && return RECURSION_THRESHOLD[]
    LoopVectorization.register_size() == 64 ? 48 : 40
end

recurse(::StridedArray) = true
recurse(_) = false

function lu!(
    A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
    pivot = Val(true), thread = Val(true);
    check::Bool=true,
    # the performance is not sensitive wrt blocksize, and 8 is a good default
    blocksize::Integer=length(A) ≥ 40_000 ? 8 : 16,
    threshold::Integer=pick_threshold()
) where T
    pivot = normalize_pivot(pivot)
    info = zero(BlasInt)
    m, n = size(A)
    mnmin = min(m, n)
    if recurse(A) && mnmin > threshold
        if T <: Union{Float32,Float64}
            GC.@preserve ipiv A begin
                info = recurse!(PtrArray(A), pivot, m, n, mnmin, PtrArray(ipiv), info, blocksize, thread)
            end
        else
            info = recurse!(A, pivot, m, n, mnmin, ipiv, info, blocksize, thread)
        end
    else # generic fallback
        info = _generic_lufact!(A, pivot, ipiv, info)
    end
    check && checknonsingular(info)
    LU{T, typeof(A)}(A, ipiv, info)
end

@inline function recurse!(A, ::Val{Pivot}, m, n, mnmin, ipiv, info, blocksize, ::Val{true}) where {Pivot}
  if length(A) * _sizeof(eltype(A)) > 0.92 * LoopVectorization.VectorizationBase.cache_size(Val(1))
    _recurse!(A, Val{Pivot}(), m, n, mnmin, ipiv, info, blocksize, Val(true))
  else
    _recurse!(A, Val{Pivot}(), m, n, mnmin, ipiv, info, blocksize, Val(false))
  end
end
@inline function recurse!(A, ::Val{Pivot}, m, n, mnmin, ipiv, info, blocksize, ::Val{false}) where {Pivot}
  _recurse!(A, Val{Pivot}(), m, n, mnmin, ipiv, info, blocksize, Val(false))
end
@inline function _recurse!(A, ::Val{Pivot}, m, n, mnmin, ipiv, info, blocksize, ::Val{Thread}) where {Pivot,Thread}
  info = reckernel!(A, Val(Pivot), m, mnmin, ipiv, info, blocksize, Val(Thread))
  @inbounds if m < n # fat matrix
    # [AL AR]
    AL = @view A[:, 1:m]
    AR = @view A[:, m+1:n]
    apply_permutation!(ipiv, AR, Val(Thread))
    ldiv!(UnitLowerTriangular(AL), AR, Val(Thread))
  end
  info
end

@inline function nsplit(::Type{T}, n) where T
    k = 512 ÷ (isbitstype(T) ? sizeof(T) : 8)
    k_2 = k ÷ 2
    return n >= k ? ((n + k_2) ÷ k) * k_2 : n ÷ 2
end

function apply_permutation!(P, A, ::Val{true})
    batchsize = cld(2000, length(P))
    @batch minbatch=batchsize for j in axes(A, 2)
        @inbounds for i in axes(P, 1)
            i′ = P[i]
            tmp = A[i, j]
            A[i, j] = A[i′, j]
            A[i′, j] = tmp
        end
    end
    nothing
end
_sizeof(::Type{T}) where {T} = Base.isbitstype(T) ? sizeof(T) : sizeof(Int)
Base.@propagate_inbounds function apply_permutation!(P, A, ::Val{false})
    for i in axes(P, 1)
        i′ = P[i]
        i′ == i && continue
        @simd for j in axes(A, 2)
            tmp = A[i, j]
            A[i, j] = A[i′, j]
            A[i′, j] = tmp
        end
    end
    nothing
end

function reckernel!(A::AbstractMatrix{T}, pivot::Val{Pivot}, m, n, ipiv, info, blocksize, thread)::BlasInt where {T,Pivot}
    @inbounds begin
        if n <= max(blocksize, 1)
            info = _generic_lufact!(A, Val(Pivot), ipiv, info)
            return info
        end
        n1 = nsplit(T, n)
        n2 = n - n1
        m2 = m - n1

        # ======================================== #
        # Now, our LU process looks like this
        # [ P1 ] [ A11 A21 ]   [ L11 0 ] [ U11 U12  ]
        # [    ] [         ] = [       ] [          ]
        # [ P2 ] [ A21 A22 ]   [ L21 I ] [ 0   A′22 ]
        # ======================================== #

        # ======================================== #
        # Partition the matrix A
        # [AL AR]
        AL = @view A[:, 1:n1]
        AR = @view A[:, n1+1:n]
        #  AL  AR
        # [A11 A12]
        # [A21 A22]
        A11 = @view A[1:n1, 1:n1]
        A12 = @view A[1:n1, n1+1:n]
        A21 = @view A[n1+1:m, 1:n1]
        A22 = @view A[n1+1:m, n1+1:n]
        # [P1]
        # [P2]
        P1 = @view ipiv[1:n1]
        P2 = @view ipiv[n1+1:n]
        # ========================================

        #   [ A11 ]   [ L11 ]
        # P [     ] = [     ] U11
        #   [ A21 ]   [ L21 ]
        info = reckernel!(AL, Val(Pivot), m, n1, P1, info, blocksize, thread)
        # [ A12 ]    [ P1 ] [ A12 ]
        # [     ] <- [    ] [     ]
        # [ A22 ]    [ 0  ] [ A22 ]
        Pivot && apply_permutation!(P1, AR, thread)
        # A12 = L11 U12  =>  U12 = L11 \ A12
        ldiv!(UnitLowerTriangular(A11), A12, thread)
        # Schur complement:
        # We have A22 = L21 U12 + A′22, hence
        # A′22 = A22 - L21 U12
        #mul!(A22, A21, A12, -one(T), one(T))
        schur_complement!(A22, A21, A12)
        # record info
        previnfo = info
        # P2 A22 = L22 U22
        info = reckernel!(A22, Val(Pivot), m2, n2, P2, info, blocksize, thread)
        # A21 <- P2 A21
        Pivot && apply_permutation!(P2, A21, thread)

        info != previnfo && (info += n1)
        @turbo warn_check_args=false for i in 1:n2
            P2[i] += n1
        end
        return info
    end # inbounds
end

function schur_complement!(𝐂, 𝐀, 𝐁)
    @tturbo warn_check_args=false for m ∈ 1:size(𝐀,1), n ∈ 1:size(𝐁,2)
        𝐂ₘₙ = zero(eltype(𝐂))
        for k ∈ 1:size(𝐀,2)
            𝐂ₘₙ -= 𝐀[m,k] * 𝐁[k,n]
        end
        𝐂[m,n] = 𝐂ₘₙ + 𝐂[m,n]
    end
end

#=
    Modified from https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/LinearAlgebra/src/lu.jl
    License is MIT: https://julialang.org/license
=#
function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where Pivot
    m, n = size(A)
    minmn = length(ipiv)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
              amax = abs(zero(eltype(A)))
              for i = k:m
                  absi = abs(A[i,k])
                  if absi > amax
                      kp = i
                      amax = absi
                  end
              end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    @simd for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                @turbo check_empty=true warn_check_args=false for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            k == minmn && break
            # Update the rest
            @turbo warn_check_args=false for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    return info
end


text:
using EffectiveWaves, Test
# using OffsetArrays

@testset "Integrated reflection" begin

    # physical parameters
    θin = 0.13
    k=1.; ho = 2

    medium = Acoustic(2; ρ=1.0, c=1.0)
    ms = MultipleScattering

    specie = Specie(Particle(
        Acoustic(2; ρ=0.1, c=0.5), ms.Circle(0.6));
        volume_fraction=0.25
    )

    ω = real(k*medium.c)

    # From effective wave theory
    k_eff0 = wavenumber_low_volumefraction(ω, medium, [specie]; basis_order = ho)
    max_x = 15.0*k/imag(k_eff0)
    x = 0.0:0.001:max_x

    normal = [-1.0,0.0] # an outward normal to the surface
    material = Material(Halfspace(normal),[specie])
    source = PlaneSource(medium, [cos(θin),sin(θin)])

    wave0 = WaveMode(ω, k_eff0, source, material; basis_order = ho)

    # wave0 = EffectivePlaneWaveMode(ω, k_eff0, medium, [specie]; θin = θin, basis_order = ho, tol=1e-8)
    wave_avg0 = DiscretePlaneWaveMode(x, wave0, material.shape)

    R = reflection_coefficient(ω, wave_avg0, source, material)
    R_eff = reflection_coefficient(ω, wave0, source, material)

    @test abs(R-R_eff) < 1e-6

    # the matched wave also gives the same reflection coefficient
    m1 = Int(round(length(x)/50))
    m2 = m1 + 100
    wave_avg1 = deepcopy(wave_avg0)
    wave_avg1.amplitudes = wave_avg0.amplitudes[1:m2,:,:]
    wave_avg1.x = wave_avg0.x[1:m2]

    match_wave = MatchPlaneWaveMode{Float64,2}([wave0], wave_avg1, x[m1:m2])
    R_m = reflection_coefficient(ω, match_wave, source, material)

    @test abs(R_m-R_eff)  < 1e-6

    num_wavenumbers = 4
    k_effs = wavenumbers(ω, medium, [specie]; tol = 1e-8,
        basis_order = ho, num_wavenumbers = num_wavenumbers)
    k_effs = k_effs[1:num_wavenumbers]

    rel_errors = map(k_effs[1:end]) do k_eff
        # println("k_eff: ", k_eff)
        wave = WaveMode(ω, k_eff, source, material;basis_order = ho)
        wave_avg = DiscretePlaneWaveMode(x, wave, material.shape)
        R = reflection_coefficient(ω, wave_avg, source, material)
        R_eff = reflection_coefficient(ω, wave, source, material)
        @test abs(R-R_eff) < 5e-5

        wave_avg1 = deepcopy(wave_avg)
        wave_avg1.amplitudes = wave_avg.amplitudes[1:m2,:,:]
        wave_avg1.x = wave_avg.x[1:m2]

        match_wave = MatchPlaneWaveMode{Float64,2}([wave], wave_avg1, x[m1:m2])
        R_m = reflection_coefficient(ω, match_wave, source, material)
        @test abs(R_m-R_eff) < 5e-5
        abs(R_m-R_eff)
    end

end

ext: jl

max_stars_repo_name: arturgower/EffectiveWaves


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a GCN-based test compiler, and generate reflection methods for it

function gcn_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(typeof(func), Base.to_tuple_type(types), kernel)
    target = GCNCompilerTarget(dev_isa="gfx900")
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function gcn_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = gcn_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function gcn_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = gcn_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function gcn_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = gcn_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function gcn_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = gcn_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    gcn_method = Symbol("gcn_$(method)")
    @eval begin
        $gcn_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $gcn_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function gcn_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = gcn_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end


text:
import Random
using LinearAlgebra
import Printf

include("../src/matrix.jl")

Random.seed!(25)

##### block inverse.
D_x = 3
D_y = D_x

N = 4

y_set = collect( randn(D_y, D_x) for n = 1:N )

Y = makeblockdiagonalmatrix(y_set)

inv_y_set = similar(y_set)
invertblockmatrix!(inv_y_set, y_set)
inv_y =  makeblockdiagonalmatrix(inv_y_set)

Y_inv = inv(Y)

println("Inverse:")
println("discrepancy between inv_y and Y_inv is ", norm(Y_inv - inv_y))
println()

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
using IndexedTables
using StructArrays
using ProgressMeter
using Flux

"""
`μlawcmp(x::T; clip = 10f0, mu = 128f0) where T <: AbstractArray`

* __input:__ array of Float32 with a mean of 0
* __args:__
  * `clip` lets you remove outliers beyond a threshold
  * `mu` lets you specify the quantization after compression e.g. -128:127
* __output:__ array of Float32
"""
function μlawcmp(x::T; clip = 1f0, mu = 128f0) where T <: AbstractArray
    μlcmp.(clamp.(x, -clip, clip)./clip, mu)
end
μlcmp(x, μ = 128f0) = sign(x) * log(1 + μ * abs(x)) / log(1 + μ)

"""
`μlawexp(x::T; clip = 10f0, mu = 128f0) where T <: AbstractArray`

* __input:__ array of Float32 with a mean of 0
* __args:__
  * `clip` lets you remove outliers beyond a threshold
  * `mu` lets you specify the quantization after expansion e.g. -128:127
* __output:__ array of Float32
"""
function μlawexp(x::T; clip = 1f0, mu = 128f0) where T <: AbstractArray
    μlexp.(x, mu) .* clip
end
μlexp(y, μ = 128f0) = sign(y) * (1/μ) * ((1 + μ)^abs(y) - 1)

"""
`preptbl(tbl; dt, lc, gc, samples, batchsize, stride, lim)`

* __input:__
  * `tbl` a JuliaDB table
* __args:__
  * `dt` the symbol for the column containing data
  * `lc` a tuple of symbols for columns with local conditioning data
  * `gc` a tuple of symbols of global conditioning classes
  * `samples` set the size of the window that you want the model to train with
  * `batchsize` set the batch size (this is the 4th dimension)
  * `stride` set the data augmentation level
        (e.g. a stride:samples ratio of 64:256 will quadruple the training data)
  * `lim` limit the size of the resulting dataset (so it fits in memory)

  Example: Wave.preptbl(tbl, dt=:sin, lc=(:lc1, :lc2), gc=(:sin, :sin2, :sin4),
            samples=1024, batchsize=8, stride=128)

* __output:__ an object ready for consumption by Flux
The format of the data created by `Wave.preptbl` below is:
[                            <== an array:    dat
    [                        <== a batch:     dat[1]
        [                    <== the data:    dat[1][1]
            x1,              <== timeseries:  dat[1][1][1]
            x2               <== condvars:    dat[1][1][2]
        ]
        [
            y                <== labels:      dat[1][2]
        ]
    ]
    ...
]
"""
function preptbl(tbl; dt, lc, gc, samples, batchsize, stride, lim = 4194304)
    tbl_ = tbl
    out = []
    println("Data points in input: $(length(tbl_))")
    println("Augmenting data in strides of $stride")
    @showprogress for i in 0:samples÷stride-1
        tbl_ = transform(tbl_, Symbol("stride_$i") =>
               [(j+i*stride)÷(samples+1) for j in 1:length(tbl_)])
        push!(out,
            select(
                groupby((idx = x->rand(), cols = x->identity(x)),
                    tbl_, Symbol("stride_$i"); select=(dt, lc...)),
            (:idx, :cols)))
    end
    sleep(0.1); println("Merging, sorting, and filtering...")
    tbl_ = reduce(merge, out)
    tbl_ = sort(tbl_, :idx)
    tbl_ = rows(filter(x -> length(x.cols) == samples+1, tbl_))
    out = []
    println("Data points after augmentation: $(length(tbl_) * samples+1)")
    println("Grouping into batches of $batchsize")
    gcvec = zeros(Float32, samples, 1, length(gc), batchsize)
    gcvec[:, :, findfirst(dt .== gc), :] .= 1f0
    @showprogress for i in 1:batchsize:length(tbl_)-length(tbl_)%batchsize
        vec = cat([reduce((x, y) -> cat(x, y, dims=3), fieldarrays(row.cols))
                  for row in tbl_[i:i+batchsize-1]]..., dims=4)
        push!(out,
            (
                (
                μlawcmp(vec[1:end-1,:,1:1,:]),
                cat(vec[1:end-1,:,2:end,:], gcvec, dims=3)
                ),
            reduce((x, y) -> hcat(x, y), map(x -> Flux.onehot(x, -128:127), round.(Int8, μlawcmp(vec[end,1,1,:]) .* 127)))
            )
        )
        if sizeof(out) > lim
            println("-- Size limit reached after $i of ",
                    "$(length(tbl_)-length(tbl_)%batchsize) rows --"); break
        end
    end
    sleep(0.1)
    println("Size of dataset: $(round(Int, sizeof(out)/1024)) MB")
    println("Training samples: $(size(out, 1) * batchsize)")
    return out
end

"""
`testmodel(m, data; dt, lc, gc)`

* __input:__
  * `m` a Wave net model
  * `data` a JuliaDB table
* __args:__
  * `dt` the symbol for the column containing data
  * `lc` a tuple of symbols for columns with local conditioning data
  * `gc` a tuple of symbols of global conditioning classes
* __output:__ an object ready for consumption by Flux
"""
function testmodel(m, data; dt, lc, gc)
    field = 2^length(m.layers)
    tbl_ = select(data, (dt, lc...))
    gcvec = zeros(Float32, field, 1, length(gc), 1)
    gcvec[:, :, findfirst(dt .== gc), :] .= 1f0
    dat = []
    lbl = Array{Float32, 1}(undef, 0)
    for i in 1:field÷2
        vec = reshape(cat(columns(tbl_[i:i+field])..., dims=3), field+1, 1, :, 1)
        push!(dat, (
            μlawcmp(vec[1:end-1,:,1:1,:]),
            cat(vec[1:end-1,:,2:end,:], gcvec, dims=3)
        ))
        push!(lbl, μlawcmp(vec[end:end,1,1,1])[1])
    end

    pred = Array{Float32, 1}(undef, 0)
    @showprogress for sample in dat
        s = vcat(reshape(sample[1][1:end-length(pred)], :), pred)
        push!(pred, cpu(Flux.onecold(
                        m((gpu(reshape(s, :, 1, 1, 1)), gpu(sample[2]))),
                    -128:127))[1] / 128)
    end
    return pred, lbl
end


text:
function HypothesisTests.ChisqTest(contab::ConTab)
    HypothesisTests.ChisqTest(contab.tab)
end

function HypothesisTests.MultinomialLRTest(contab::ConTab)
    HypothesisTests.MultinomialLRTest(contab.tab)
end

function HypothesisTests.FisherExactTest(contab::ConTab)
    if !(size(contab, 1) == size(contab, 2) == 2)  error("FisherExactTest can be applied only for 2X2 tables") end
    HypothesisTests.FisherExactTest(contab.tab[1,1], contab.tab[1,2], contab.tab[2,1], contab.tab[2,2])
end


text:
module _CoSiEn
export CoSiEn
using Statistics: std, mean, median
using LinearAlgebra: Diagonal, UpperTriangular
    """
        CoSi, Bm = CoSiEn(Sig)

    Returns the cosine similarity entropy (`CoSi`) and the corresponding
    global probabilities estimated from the data sequence (`Sig`) using the
    default parameters:   embedding dimension = 2, time delay = 1,
    angular threshold = .1,  logarithm = base 2,

        CoSi, Bm = CoSiEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, r::Real=.1, Logx::Real=2, Norm::Int=0)

    Returns the cosine similarity entropy (`CoSi`) estimated from the data
    sequence (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 1   \n
    `tau`   - Time Delay, a positive integer    \n
    `r`     - Angular threshold, a value in range [0 < r < 1]   \n
    `Logx`  - Logarithm base, a positive scalar (enter 0 for natural log) \n
    `Norm`  - Normalisation of `Sig`, one of the following integers:    \n
            [0]  no normalisation - default
            [1]  normalises `Sig` by removing median(`Sig`)
            [2]  normalises `Sig` by removing mean(`Sig`)
            [3]  normalises `Sig` w.r.t. SD(`Sig`)
            [4]  normalises `Sig` values to range [-1 1]

    # See also `PhasEn`, `SlopEn`, `GridEn`, `MSEn`, `cMSEn`

    # References:
        [1] Theerasak Chanwimalueang and Danilo Mandic,
            "Cosine similarity entropy: Self-correlation-based complexity
            analysis of dynamical systems."
            Entropy
            19.12 (2017): 652.

    """
    function CoSiEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1,
        r::Real=.1, Logx::Real=2, Norm::Int=0)

    Logx == 0  ? Logx = exp(1) : nothing

    N = size(Sig,1)
    (N > 10) ? nothing : error("Sig:   must be a numeric vector")
    (m > 1) ? nothing :  error("m:     must be an integer > 1")
    (tau>0) ? nothing :  error("tau:   must be an integer > 0")
    (0<r<1) ? nothing :  error("r:     must be a scalar in range [0 1]")
    (Logx>0) ? nothing : error("Logx:  must be a positive number > 0")
    (Norm in collect(0:4)) ? nothing : error("Norm:   must be an integer in range [0 4]")

    if Norm == 1
        Xi = Sig .- median(Sig);
    elseif Norm == 2
        Xi = Sig .- mean(Sig);
    elseif Norm == 3
        Xi = (Sig .- mean(Sig))/std(Sig,corrected=false)
    elseif Norm == 4
        Xi = (2*(Sig .- minimum(Sig))/(maximum(Sig)-minimum(Sig))) .- 1;
    else
        Xi = Sig;
    end

    Nx = N-((m-1)*tau);
    Zm = zeros(Nx,m);
    for n = 1:m
        Zm[:,n] = Xi[(n-1)*tau+1:Nx+(n-1)*tau]
    end

    Num = Zm*transpose(Zm);
    Mag = sqrt.(sum(Diagonal(Num),dims=1))[:]
    Den = Mag*transpose(Mag)
    AngDis = round.(acos.(round.(Num./Den,digits=6))/pi,digits=6)
    if maximum(imag.(AngDis)) < (10^-5)
        Bm = (sum(UpperTriangular(AngDis .< r))-Nx)/(Nx*(Nx-1)/2)
    else
        Bm = (sum(UpperTriangular(real.(AngDis) .< r))-Nx)/(Nx*(Nx-1)/2)
        @warn("Complex values ignored.")
    end
    if Bm == 1 || Bm == 0
        CoSi = 0
    else
        CoSi = -(Bm*log(Logx, Bm)) - ((1-Bm)*log(Logx, 1-Bm))
    end

    return CoSi, Bm
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using Immerse, Images, Gtk, DataFrames, ZipFile
import ImageView

const testdir = splitdir(@__FILE__)[1]
const facesdir = joinpath(testdir, "orl_faces")
const orl_url = "http://www.cl.cam.ac.uk/Research/DTG/attarchive/pub/data/att_faces.zip"

include("faces_utilities.jl")

if !isdir(facesdir)
    fn = download(orl_url)
    mkpath(facesdir)
    unzip(fn, facesdir)
end

function load_faces(parentdir=facesdir)
    imgs = Any[]
    group = Int[]
    for i = 1:40
        childdir = joinpath(parentdir, string("s", i))
        for j = 1:10
            img = load(joinpath(childdir, string(j, ".pgm")))
            push!(imgs, Images.data(img))
            push!(group, i)
        end
    end
    imgs, group
end

function run_lda(imgs, group)
    M = zeros(Float32, length(imgs[1]), length(imgs))
    for i = 1:length(imgs)
        M[:,i] = vec(imgs[i])
    end
    evec, eval = lda(M, group)
    evec[:,1:2]'*M
end

function showimg(imgs, indx)
    img = imgs[indx]
    c = @GtkCanvas()
    f = @GtkAspectFrame(c, "", 0.5, 0.5, size(img,1)/size(img,2))
    win = @GtkWindow(f, string("Face ", indx))
    showall(win)
    c.draw = function(widget)
        copy!(widget, imgs[indx])
    end
end

function showimgs(imgs, indexes=1:length(imgs))
    imgm = grayim(cat(3, imgs[indexes]...))
    ImageView.view(imgm, pixelspacing=[1,1])
end

imgs, group = load_faces()
proj = run_lda(imgs, group)
df = DataFrame(Any[vec(proj[1,:]),vec(proj[2,:]),group], [:comp1,:comp2,:group])

include("faces_run.jl")


text:
@testset "MatrixAliasSampler" begin
  @testset "Rates similar to empirical results" begin
    rng = MersenneTwister(13)
    N = 5
    example_rates = example_probs = rand(rng, N, N) * 10
    matrix_sampler = MocosSim.MatrixAliasSampler(example_probs)

    num_samples = 10^7
    counts = zeros(Int, N, N)

    for src in 1:N
      for _ in 1:num_samples
        tgt = MocosSim.sample(rng, matrix_sampler, src)
        counts[tgt, src] += 1
      end
    end

    empirial_probs = counts ./ num_samples
    theoretical_probs = example_rates ./ sum(example_rates, dims=1)

    @test maximum(abs.(empirial_probs .- theoretical_probs)) < 3/sqrt(num_samples)
  end
end


text:
function localSearch3(estimator::SBMEstimator, dataset::Dataset)::Tuple{SBM, Matrix{Int}, OptResults}
    """ Local Search 3
    1. Initialize with random assignments;
    2. Solve analytically for the optimal probability matrix w;
    3. Solve quadratic program for the optimal assignments, given w;
    4. Repeat step 2, until no more improvement can be made.

    Parameters
    ----------
    estimator : SBMEstimator
        Optimization method specifications.
    dataset : Dataset
        Dataset representing an observed graph.

    Returns
    -------
    sbm : SBM
        Stochastic Block Model
    x : Array{Int64,2}
        Matrix of assignments of nodes to groups
    opt_results : OptResults
        Results of the optimization process
    """
    timeLimit = estimator.time_limit

    start = time(); # Start counting time

    # Initialize with random assignments
    x = randomAssignments(dataset, seed=estimator.seed)
    w = optimalProbMatrix(dataset, x)

    # Calculate objective value
    obj_value = calculateObjective(dataset, w, x)
    best_obj = obj_value

    status = nothing
    iterations = 0
    improved = true
    while true
        ## Base cases
        # Break if there is no improvement
        if ~improved
            status = :LocalOptimum_LS3
            break
        end
        improved = false

        # Break if time limit is exceeded
        availableTime = timeLimit - (time() - start)
        if availableTime <= 0
            status = :UserLimit
            break
        end

        # Optimal Assignments
        sbm = SBM(w, "poisson")
        assignment_results, x_opt = optimalAssignments(dataset, sbm, time_limit=availableTime,
                                                       verbose=estimator.verbose)
        if assignment_results.status == :UserLimit
            status = :UserLimit
            break
        elseif assignment_results.status == :Optimal
            x = x_opt
        end

        # Update w with optimal value
        w = optimalProbMatrix(dataset, x)

        # Evaluate move
        obj_value = calculateObjective(dataset, w, x)

        # Check if there was an improvement
        if obj_value < best_obj
            best_obj = obj_value
            improved = true
        end
        iterations += 1 # Update counter
    end
    solvetime = time() - start;

    # SBM
    sbm = SBM(w, "poisson") # assumes a poisson distribution

    # Build OptResults
    obj_lb = -Inf
    obj_ub = obj_value
    nodecount = nothing
    lazycount = nothing
    opt_results = OptResults(obj_lb, obj_ub, status, solvetime, iterations, nodecount, lazycount)
    if estimator.verbose
        displayResults(opt_results)
    end

    return sbm, x, opt_results
end


text:
"""
    solve_pattern_generation(d::StudentHousingData)

Solve generalized assignment problem with column generation.
If `suboptimal_rc` is true, we generate the first column we find with negative
reduced cost.
"""
function solve_pattern_generation(d::StudentHousingData, suboptimal_rc::Bool=false)

    nhouses = length(d.houses)
    npatterns = length(d.patterns_allow)

    p_init = collect(1:npatterns)
    # Create an initial set of columns by randomly picking the first entity
    # allowed for any house
    U = zeros(nhouses, npatterns)
    for p = 1:npatterns
        # If there is a pattern that doesn't match any houses, skip it
        if isempty(d.patterns_allow[p]) || (d.demands[p, 1, 1] == 0)
            continue
        end
        # A naive initial assignment
        # Find the first house allowed
        for i = 1:nhouses
            if house_fits_pattern(i, p, d)
                U[i, p] = 1.0
                break
            end
            # If we don't want to check all patterns, only add 1 col
            if suboptimal_rc
                break
            end
        end
    end

    # Our cost vector
    c = ones(nhouses)

    # =========================================================================
    # Master problem
    # =========================================================================
    m = Model(solver = GurobiSolver(OutputFlag=0))
    @variable(m, 0 <= λ[i in p_init] <= 1)
    @objective(m, Max, sum(λ[p] * dot(c, U[:, p]) for p in p_init))
    @constraints(m, begin
        beds[i = 1:nhouses], sum(λ[p] * U[i, p] * beds_needed(p, d) for p in p_init) <= beds_avail(i, d.houses)
        convexity[p in p_init], sum(λ[p]) <= 1
    end)

    # =========================================================================
    # The knapsack subproblem for each pattern
    # =========================================================================
    function get_kp(p::Int, π::Vector{Float64})
        sp = Model(solver = GurobiSolver(OutputFlag=0))
        @variable(sp, u[i = 1:nhouses] >= 0, Int)
        @constraints(sp, begin
            sum(u[i] for i = 1:nhouses) <= d.demands[p, 1, 1]
            [i = 1:nhouses], u[i] <= house_fits_pattern(i, p, d)
        end)
        @objective(sp, Max, dot((c - π * beds_needed(p, d)), u))
        sp
    end

    # Cache information about columns we generate
    U_generated = Vector{Float64}[]
    pattern_choice = Int[]
    λ_generated = Variable[]

    best_sp  = Model()
    best_p = 0
    best_rc = Inf
    iter = 1

    while best_rc > 0
        # Solve the restricted master problem
        @assert solve(m) == :Optimal
        π = getdual(m[:beds])
        w = getdual(m[:convexity])

        # Solve a subproblem for each pattern, get the column with the best reduced
        # cost over all patterns
        best_rc = -Inf
        for p = 1:npatterns
            # Get the subproblem
            sp = get_kp(p, π)
            @assert solve(sp) == :Optimal
            col_rc = getobjectivevalue(sp) - w[p]
            if col_rc > best_rc
                best_rc = col_rc
                best_sp = sp
                best_p = p
            end
            # If we don't care about the optimal column
            if suboptimal_rc
                # ... and we found one with positive reduced cost
                if col_rc > 1e-4
                    # just add the one we found
                    break
                end
            end
        end

        best_u = getvalue(best_sp[:u])
        push!(U_generated, best_u)
        push!(pattern_choice, best_p)

        # Add the best column
        push!(λ_generated, @variable(m,
                lowerbound = 0.0,
                upperbound = 1.0,
                basename="_λnew_$iter",
                objective = dot(c, best_u),
                inconstraints = [beds; convexity[best_p]],
                coefficients = [best_u * beds_needed(best_p, d); 1.0]
              )
        )
        iter += 1
    end
    # Solve one last time
    @assert solve(m) == :Optimal

    # Return the master problem and information about generated columns
    m, U_generated, λ_generated, pattern_choice
end


text:
using Polynomials
using Plots
include("funkcje.jl")

function makeSpline(coordinates)
    function CalculateCs()
        μ = [NaN64 for k in 1:n]
        α = vcat([0.], [3. * (A[k+1] - A[k])/h[k] - 3. * (A[k] - A[k-1])/h[k-1] for k in 2:n])

        C = [0. for k in 1:nPlus1]
        l = [0. for k in 1:nPlus1]
        z = [0. for k in 1:nPlus1]

        l[1] = 1.
        μ[1] = 0.
        z[1] = 0.

        for k in 2:n
            l[k] = 2 * (xs[k+1] - xs[k-1]) - h[k-1] * μ[k-1]
            μ[k] = h[k] / l[k]
            z[k] = (α[k] - h[k-1]*z[k-1]) / l[k]
        end

        l[nPlus1] = 1.
        z[nPlus1] = 0.
        C[nPlus1] = 0.

        for k in n:(-1):1
            C[k] = z[k] - μ[k]*C[k+1]
        end
        return C
    end
    nPlus1 = length(coordinates)
    n = nPlus1 - 1
    # x1 x2 ... xk+1
    xs = map(pair->pair[1], coordinates)
    # y1 y2 ... yk+1
    # ys = map(pair->pair[2], coordinates)
    A = map(pair->pair[2], coordinates)

    h = [xs[k+1] - xs[k] for k in 1:n]
    C = CalculateCs()

    B = [(A[k+1] - A[k]) / h[k] - h[k]*(C[k+1] + 2. * C[k]) / 3. for k in 1:n]
    D = [(C[k+1] - C[k]) / (3. * h[k]) for k in 1:n]

    outputSet = [Poly(NaN) for k in 1:n]
    intervals = myzip(xs,xs[2:end])

    for k in 1:n
        xk = xs[k]
        polynom = A[k] + B[k]*poly([xk]) + C[k]*(poly([xk])^2) + D[k]*(poly([xk]) ^ 3)
        outputSet[k] = polynom
    end
    return myzip(outputSet, intervals)
end

function findRoot(polyn, c)
    interval = polyn[2]
    xk = interval[1]
    xkPlus1 = interval[2]
    polynomial = polyn[1]

    rooty = roots(polynomial - c)
    filter!(z -> imag(z) == 0., rooty)
    rooty = map(z -> real(z), rooty)

    filter!(x -> xk <= x <= xkPlus1, rooty)
    return rooty
end

function findRoots(splines, c)
    Rootki = map(p -> findRoot(p,c), splines)
    filter!(r -> length(r) != 0, Rootki)
    return Rootki
end

function PlotMain(coor, splajn, pierwiastki, c, a, b, name, prec = 0.01)
    function PlotSpline()
        xs = []
        ys = []
        # ysPrim = []
        # ysPrimPrim = []
        for i in 1:length(splajn)
            xk = splajn[i][2][1]
            xkPlus1 = splajn[i][2][2]
            wielomian = splajn[i][1]
            xs = vcat(xs, [x for x in xk:prec:xkPlus1])
            ys = vcat(ys, map(wielomian, [x for x in xk:prec:xkPlus1]))
            # ysPrim = vcat(ysPrim, map(polyder(wielomian), [x for x in xk:prec:xkPlus1]))
            # ysPrimPrim = vcat(ysPrimPrim, map(polyder(polyder(wielomian)), [x for x in xk:prec:xkPlus1]))
        end
        plot!(plotka, xs, ys, lab = "Naturalna funkcja sklejana s(x)", linecolor=:blue)
        # plot!(plotka, xs, ysPrim, lab = "Naturalna funkcja sklejana s'(x)", linecolor=:green)
        # plot!(plotka, xs, ysPrimPrim, lab = "Naturalna funkcja sklejana s''(x)", linecolor=:black)
    end
    function PlotRoots()
        scatter!(plotka, pierwiastki, map(x->c, pierwiastki), lab = "Pierwiastki równania f(x) = c", color=:orange, markersize = 11)
    end
    function PlotC()
        plot!(plotka, [a,b], [c,c], lab = "y=c",  linecolor=:black, linewidth=5)
    end
    function PlotF()
        xs, ys = myunzip(coor)
        scatter!(plotka, xs, ys, lab = "Zadane punkty (x,f(x))", color=:green)
    end

    plotka = plot(legend = :topright)
    xlabel!(plotka, "Os x")
    ylabel!(plotka, "Os y")
    title!(string("Wykres s(x), y=c, pierwiastków s(x) = y", name))
    PlotSpline()
    PlotRoots()
    PlotC()
    PlotF()
    display(plotka)
    savefig(plotka, name)
end

function main(coor, c, name="")
    coor = sort(coor)
    splajn = makeSpline(coor)
    pierwiastki = VectorToList(findRoots(splajn, c))

    a = coor[1][1]
    b = coor[end][1]

    PlotMain(coor, splajn, pierwiastki, c, a, b, name)
    return pierwiastki
end

max_stars_repo_path: II rok/I semestr/ANM/Pracownia 2/prog/program.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
unshift!(LOAD_PATH, ".", "../src")

using Comm, Par
using Cells, Defs, Domains, Grids, Memoizeds, Norms

function main()
    println("Wave[FunHPC.jl]")

    println("Running on $(Comm.nprocs()) processes")

    filename = "wave.$(Comm.nprocs()).txt"
    file = open(filename, "w")

    println("Initialization")
    @time begin
        s = par(R=Domain) do; Domains.initial(Defs.tmin) end
        m = Memoized(0, s)
        s1 = rk2(m)             # to warm up
        info_output(STDOUT, m, is_first=true)
        file_output(file, m, is_first=true)
        wait(m)
        wait(s1[])
    end

    println("Evolution")
    # TODO: Don't check s.t at every iteration
    @time begin
        while ((Defs.nsteps<0 || m.n<Defs.nsteps) &&
               (Defs.tmax<0.0 || s[].t < Defs.tmax - 0.5*Defs.dt()))
            s = rk2(m)
            m = Memoized(m.n+1, s)
            info_output(STDOUT, m)
            file_output(file, m)
        end
        info_output(STDOUT, m, was_last=true)
        file_output(file, m, was_last=true)
    end

    close(file)

    println("Done.")
end

run_main(main)


text:
function step_size(x0, zero_tol = 1e-15)
    #=
    if abs(x0) <= zero_tol
        x0 = copysign(zero_tol, x0)
    end
    =#
    h = sqrt( eps(x0) )
    return h
end

function step_size_cbrt(x0, zero_tol = 1e-15)
    #=
    if abs(x0) <= zero_tol
        x0 = copysign(zero_tol, x0)
    end
    =#
    h = cbrt( eps(x0) )
    return h
end

function numdiff_fw(f, x0)
    # f: R -> R
    #=
    h = step_size(x0)
    x1 = x0 + h
    f0 = f(x0)
    f1 = f(x1)
    g = (f1 - f0) / h
    return g
    =#
    h = step_size(x0)
    numdiff_fw(f, x0, h)
end

function numdiff_fw(f, x0, h)
    # f: R -> R
    f0 = f(x0)
    numdiff_fw(f, x0, h, f0)
end

function numdiff_fw(f, x0, h, f0)
    # f: R -> R
    x1 = x0 + h
    f1 = f(x1)
    g = (f1 - f0) / h
end

function numdiff_bw(f, x0)
    # f: R -> R
     #=
    h = step_size(x0)
    x1 = x0 - h
    f0 = f(x0)
    f1 = f(x1)
    g = (f0 - f1) / h
    =#
    h = step_size(x0)
    numdiff_bw(f, x0, h)
end

function numdiff_bw(f, x0, h)
    f0 = f(x0)
    numdiff_bw(f, x0, h, f0)
end

function numdiff_bw(f, x0, h, f0)
    x1 = x0 - h
    f1 = f(x1)
    g = (f0 - f1) / h
end

function numdiff_central(f, x0)
    # f: R -> R
    h = step_size_cbrt(x0)
    xm1 = x0 - h
    xp1 = x0 + h
    fm1 = f(xm1)
    fp1 = f(xp1)
    g = (fp1 - fm1) / (2*h)
    return g
end

function gradient_fw(f, x0)
    # f: R^n -> R
    # g_i = df / dx_i
    h = step_size.(x0)
    return gradient_fw(f, x0, h)
end

function gradient_fw(f, x0, h)
    # f: R^n -> R
    # g_i = df / dx_i
    f0 = f(x0)
    return gradient_fw(f, x0, h, f0)
end

function gradient_fw(f, x0, h, f0)
    # f: R^n -> R
    # g_i = df / dx_i
    g = similar(x0)
    x1 = copy(x0)
    for kk = 1:length(x0)
        x1[kk] += h[kk]
        f1 = f(x1)
        g[kk] = (f1 - f0) / h[kk]
        x1[kk] -= h[kk]
    end
    return g
end

function directional_diff_fw(f, x0, d)
    h = step_size( norm(d) )
    return directional_diff_fw(f, x0, d, h)
end

function directional_diff_fw(f, x0, d, h)
    f0 = f(x0)
    return directional_diff_fw(f, x0, d, h, f0)
end

function directional_diff_fw(f, x0, d, h, f0)

    g = zero(eltype(x0))

    x1 = x0 + d * h
    f1 = f(x1)

    g = (f1 - f0) / h

    return g
end

function jacobian_fw(f, x0)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    h = step_size.(x0)  # n-by-1
    return jacobian_fw(f, x0, h)
end

function jacobian_fw(f, x0, h)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    f0 = f(x0)                  # m-by-1
    return jacobian_fw(f, x0, h, f0)   # m-by-n
end

function jacobian_fw(f, x0, h, f0)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    jac = zeros(eltype(f0), length(f0), length(x0)) # m-by-n
    x1 = copy(x0)
    for kk = 1:length(x0)
        x1[kk] += h[kk]
        f1 = f(x1)
        jac[:, kk] = (f1 .- f0) ./ h[kk]
        x1[kk] -= h[kk]
    end
    return jac
end

function jacobian_fw!(jac, f, x0)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    h = step_size.(x0)  # n-by-1
    return jacobian_fw!(jac, f, x0, h)
end

function jacobian_fw!(jac, f, x0, h)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    f0 = f(x0)                  # m-by-1
    return jacobian_fw!(jac, f, x0, h, f0)   # m-by-n
end

function jacobian_fw!(jac, f, x0, h, f0)
    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    f1 = similar(f0)
    x1 = copy(x0)
    for kk = 1:length(x0)
        x1[kk]     += h[kk]
        f1         .= f(x1)
        jac[:, kk] .= (f1 .- f0) ./ h[kk]
        x1[kk]     -= h[kk]
    end
    return jac
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
const SpatialDim = Union{XDim,YDim,ZDim}

"""
    Band <: Dimension

    Band(val=:)

Band [`Dimension`]($DDdimdocs) for multi-band rasters.

## Example:
```julia
banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)
```
"""
@dim Band


text:
module AMQPTestRPC

using AMQPClient, Test, Random

const JULIA_HOME = Sys.BINDIR

const QUEUE_RPC = "queue_rpc"

testlog(msg) = println(msg)
reply_queue_id = 0
server_id = 0

const NRPC_MSGS = 100
const NRPC_CLNTS = 4
const NRPC_SRVRS = 4

function test_rpc_client(;virtualhost="/", host="localhost", port=AMQPClient.AMQP_DEFAULT_PORT, auth_params=AMQPClient.DEFAULT_AUTH_PARAMS)
    # open a connection
    testlog("client opening connection...")
    conn = connection(;virtualhost=virtualhost, host=host, port=port, auth_params=auth_params)

    # open a channel
    testlog("client opening channel...")
    chan1 = channel(conn, AMQPClient.UNUSED_CHANNEL, true)

    # create a reply queue for a client
    global reply_queue_id
    reply_queue_id += 1
    queue_name = QUEUE_RPC * "_" * string(reply_queue_id) * "_" * string(getpid())
    testlog("client creating queue " * queue_name * "...")
    success, message_count, consumer_count = queue_declare(chan1, queue_name; exclusive=true)

    testlog("client testing rpc...")
    rpc_reply_count = 0
    rpc_fn = (rcvd_msg) -> begin
        rpc_reply_count += 1

        msg_str = String(rcvd_msg.data)
        println("client ", msg_str)

        basic_ack(chan1, rcvd_msg.delivery_tag)
    end

    # start a consumer task
    success, consumer_tag = basic_consume(chan1, queue_name, rpc_fn)
    @test success

    correlation_id = 0

    # publish NRPC_MSGS messages to the queue
    while correlation_id < NRPC_MSGS
        correlation_id += 1
        M = Message(Vector{UInt8}("hello from " * queue_name), content_type="text/plain", delivery_mode=PERSISTENT, reply_to=queue_name, correlation_id=string(correlation_id))
        basic_publish(chan1, M; exchange=default_exchange_name(), routing_key=QUEUE_RPC)
        # sleep a random time between 1 and 5 seconds between requests
        sleep(rand())
    end

    while (rpc_reply_count < NRPC_MSGS)
        sleep(1)
    end

    testlog("client closing down...")
    success, message_count = queue_purge(chan1, queue_name)
    @test success
    @test message_count == 0

    @test basic_cancel(chan1, consumer_tag)

    success, message_count = queue_delete(chan1, queue_name)
    @test success
    @test message_count == 0

    # close channels and connection
    close(chan1)
    AMQPClient.wait_for_state(chan1, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(chan1)

    close(conn)
    AMQPClient.wait_for_state(conn, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(conn)

    testlog("client done.")
end

function test_rpc_server(;virtualhost="/", host="localhost", port=AMQPClient.AMQP_DEFAULT_PORT, auth_params=AMQPClient.DEFAULT_AUTH_PARAMS)
    global server_id
    server_id += 1
    my_server_id = server_id

    # open a connection
    testlog("server $my_server_id opening connection...")
    conn = connection(;virtualhost=virtualhost, host=host, port=port, auth_params=auth_params)

    # open a channel
    testlog("server $my_server_id opening channel...")
    chan1 = channel(conn, AMQPClient.UNUSED_CHANNEL, true)

    # create queues (no need to bind if we are using the default exchange)
    testlog("server $my_server_id creating queues...")
    # this is the callback queue
    success, message_count, consumer_count = queue_declare(chan1, QUEUE_RPC)
    @test success

    # test RPC
    testlog("server $my_server_id testing rpc...")
    global rpc_count = 0
    rpc_fn = (rcvd_msg) -> begin
        global rpc_count
        rpc_count += 1

        @test :reply_to in keys(rcvd_msg.properties)
        reply_to = convert(String, rcvd_msg.properties[:reply_to])
        correlation_id = convert(String, rcvd_msg.properties[:correlation_id])

        resp_str = "$(my_server_id) received msg $(rpc_count) - $(reply_to): $(String(rcvd_msg.data))"
        println("server ", resp_str)

        M = Message(Vector{UInt8}(resp_str), content_type="text/plain", delivery_mode=PERSISTENT, correlation_id=correlation_id)
        basic_publish(chan1, M; exchange=default_exchange_name(), routing_key=reply_to)

        basic_ack(chan1, rcvd_msg.delivery_tag)
    end

    # start a consumer task
    success, consumer_tag = basic_consume(chan1, QUEUE_RPC, rpc_fn)
    @test success

    while (rpc_count < NRPC_MSGS*NRPC_CLNTS)
        sleep(1)
    end

    testlog("server $my_server_id closing down...")
    success, message_count = queue_purge(chan1, QUEUE_RPC)
    @test success
    @test message_count == 0

    @test basic_cancel(chan1, consumer_tag)

    success, message_count = queue_delete(chan1, QUEUE_RPC)
    @test success
    @test message_count == 0
    testlog("server $my_server_id deleted rpc queue")

    # close channels and connection
    close(chan1)
    AMQPClient.wait_for_state(chan1, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(chan1)

    close(conn)
    AMQPClient.wait_for_state(conn, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(conn)

    testlog("server $my_server_id done.")
    nothing
end

function runtests()
    testlog("testing multiple client server rpc")
    clients = Vector{Task}()
    servers = Vector{Task}()

    for idx in 1:NRPC_SRVRS
        push!(servers, @async test_rpc_server())
    end

    for idx in 1:NRPC_CLNTS
        push!(clients, @async test_rpc_client())
    end

    tasks_active = NRPC_CLNTS + NRPC_SRVRS
    while tasks_active > 0
        tasks_active = 0
        for idx in 1:NRPC_SRVRS
            istaskdone(servers[idx]) || (tasks_active += 1)
        end
        for idx in 1:NRPC_CLNTS
            istaskdone(clients[idx]) || (tasks_active += 1)
        end
        sleep(5)
    end
    testlog("done")
end
end # module AMQPTestRPC


text:
using Printf
using TimerOutputs
using Oceananigans
using Oceananigans.Advection
using Oceananigans.Utils

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Nt = 10  # Number of iterations to use for benchmarking time stepping.

   float_types = [Float64]       # Float types to benchmark.
         archs = [CPU()]         # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]  # Benchmark GPU on systems with CUDA-enabled GPUs.

Schemes = (CenteredSecondOrder, CenteredFourthOrder, UpwindBiasedThirdOrder, WENO5)

#####
##### Run benchmarks
#####

for arch in archs, FT in float_types, Scheme in Schemes
    N = arch isa CPU ? 64 : 256

    topo = (Periodic, Periodic, Periodic)
    grid = RegularCartesianGrid(FT, topology=topo, size=(N, N, N), halo=(3, 3, 3), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, float_type=FT, grid=grid, advection=Scheme())

    time_step!(model, 1)  # precompile

    bn =  benchmark_name((N, N, N), string(Scheme), arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(model, 1)
    end
end

#####
##### Print benchmark results
#####

println()
println(oceananigans_versioninfo())
println(versioninfo_with_gpu())
print_timer(timer, title="Advection scheme benchmarks", sortby=:name)
println()


text:
push!(LOAD_PATH, "..")

using
  Documenter,
  Literate,
  Plots,  # to not capture precompilation output
  GeophysicalFlows

# Gotta set this environment variable when using the GR run-time on Travis CI.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

examples = [
    "twodnavierstokes_decaying.jl",
    "twodnavierstokes_stochasticforcing.jl",
    "twodnavierstokes_stochasticforcing_budgets.jl",
    "barotropicqg_betadecay.jl",
    "barotropicqg_betaforced.jl",
    "barotropicqgql_betaforced.jl",
    "barotropicqg_decay_topography.jl",
    "multilayerqg_2layer.jl",
    "surfaceqg_decaying.jl",
]

for example in examples
  example_filepath = joinpath(EXAMPLES_DIR, example)
  withenv("GITHUB_REPOSITORY" => "FourierFlows/GeophysicalFlowsDocumentation") do
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true)
    Literate.notebook(example_filepath, OUTPUT_DIR, documenter=true)
    Literate.script(example_filepath, OUTPUT_DIR, documenter=true)
  end
end

#####
##### Build and deploy docs
#####

# Set up a timer to print a space ' ' every 240 seconds. This is to avoid Travis CI
# timing out when building demanding Literate.jl examples.
Timer(t -> println(" "), 0, interval=240)

format = Documenter.HTML(
  collapselevel = 2,
     prettyurls = get(ENV, "CI", nothing) == "true",
      canonical = "https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/"
)

makedocs(
 modules = [GeophysicalFlows],
 doctest = false,
   clean = true,
checkdocs = :all,
  format = format,
 authors = "Navid C. Constantinou and Gregory L. Wagner",
sitename = "GeophysicalFlows.jl",
   pages = Any[
            "Home"    => "index.md",
            "Examples" => [
              "TwoDNavierStokes" => Any[
                "generated/twodnavierstokes_decaying.md",
                "generated/twodnavierstokes_stochasticforcing.md",
                "generated/twodnavierstokes_stochasticforcing_budgets.md",
                ],
              "BarotropicQG" => Any[
                "generated/barotropicqg_betadecay.md",
                "generated/barotropicqg_betaforced.md",
                "generated/barotropicqg_decay_topography.md"
                ],
              "BarotropicQGQL" => Any[
                "generated/barotropicqgql_betaforced.md",
                ],
              "MultilayerQG" => Any[
                "generated/multilayerqg_2layer.md"
                ],
              "SurfaceQG" => Any[
                "generated/surfaceqg_decaying.md"
                ]
            ],
            "Modules" => Any[
              "modules/twodnavierstokes.md",
              "modules/barotropicqg.md",
              "modules/barotropicqgql.md",
              "modules/multilayerqg.md",
              "modules/surfaceqg.md"
            ],
            "Forcing" => "forcing.md",
            "DocStrings" => Any[
            "man/types.md",
            "man/functions.md"]
           ]
)

withenv("GITHUB_REPOSITORY" => "FourierFlows/GeophysicalFlowsDocumentation") do
  deploydocs(        repo = "github.com/FourierFlows/GeophysicalFlowsDocumentation.git",
                versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
               forcepush = true,
            push_preview = true
            )
end


text:
const SHIFT = 1
const REDUCEL = 2
const REDUCER = 3

# calculates action id for labeled reduce
# even numbers >= 2
reducel(label::Int) = label << 1
reducel() = REDUCEL

 # odd numbers >= 3
reducer(label::Int) = (label << 1) | 1
reducer() = REDUCER

# retrieve action type 1, 2 or 3
actiontype(action::Int) = action == 1 ? SHIFT : 2 + (action & 1)

# retrieve label id
tolabel(action::Int) = action >> 1

type State{T <: ParserType}
    step::Int
    score::Float64
    top::Int
    right::Int
    left
    lchild
    rchild
    lsibl
    rsibl
    tokens::Vector{Token}
    parser::DepParser
    prev::State{T}
    prevact::Int
    feat::Vector{Int}

    function State(step, score, top, right)
        new(step, score, top, right)
    end

    function State(step, score, top, right, left,
        lchild, rchild, lsibl, rsibl, tokens, parser)
        new(step, score, top, right, left,
            lchild, rchild, lsibl, rsibl, tokens, parser)
    end

    function State(step, score, top, right, left,
        lchild, rchild, lsibl, rsibl, tokens, parser, prev, prevact)
        new(step, score, top, right, left,
            lchild, rchild, lsibl, rsibl, tokens, parser, prev, prevact)
    end

end

function __nullstate()
    s = State(0, 0.0, 0, 1)
    s.left = s
    s.lchild, s.rchild = s, s
    s.lsibl, s.rsibl = s, s
    return s
end

const nullstate = __nullstate()

function State{T}(tokens::Vector{Token}, parser::DepParser{T})
    State{T}(1, 0.0, 0, 1, nothing, nothing,
        nothing, nothing, nothing, tokens, parser)
end

# to retrieve result
function heads(s::State)
    @assert isfinal(s)
    res = fill(-1, length(s.tokens))
    st = s
    while isdefined(st, :prev)
        st.lchild != nothing && ( res[st.lchild.top] = st.top )
        st.rchild != nothing && ( res[st.rchild.top] = st.top )
        st = st.prev
    end
    @assert all(h -> h >= 0, res)
    res
end

function labels(s::State)
    @assert isfinal(s)
    res = fill(-1, length(s.tokens))
    st = s
    while isdefined(st, :prev)
        acttype = actiontype(st.prevact)
        acttype == REDUCEL && ( res[st.lchild.top] = tolabel(st.prevact) )
        acttype == REDUCER && ( res[st.rchild.top] = tolabel(st.prevact) )
        st = st.prev
    end
    @assert all(h -> h >= 0, res)
    res
end

tokenat(s::State, i::Int) = get(s.tokens, i, rootword)
tokenat(s::State, t::State) = get(s.tokens, t.top, rootword)
tokenat(s::State, ::Void) = rootword

function tokenat(s::State, head::Symbol, tail::Symbol...)
    if isempty(tail) || s.(head) == nothing
        return tokenat(s, s.(head))
    else
        return tokenat(s.(head), tail...)
    end
end

function labelat(s::State, child::State)
    st = s
    while st.step != child.step + 1
        isdefined(st, :prev) || return -1
        st = st.prev
    end
    tolabel(st.prevact)
end

###################################################
#################### "expand"s ####################
###################################################

function transit{T}(s::State{T}, act::Int, top::Int,
    right::Int, left, lchild, rchild, lsibl, rsibl)
    score = s.score + s.parser.model(s, act)
    State{T}(s.step + 1, score, top, right, left, lchild,
        rchild, lsibl, rsibl, s.tokens, s.parser, s, act)
end

function expand(s::State, act::Int)
    acttype = actiontype(act)
    if acttype == SHIFT
        return transit(s, act, s.right, s.right+1, s, nothing, nothing, nothing, nothing)
    elseif acttype == REDUCEL
        return transit(s, act, s.top, s.right, s.left.left, s.left, s.rchild, s, s.rsibl)
    elseif acttype == REDUCER
        return transit(s, act, s.left.top, s.right, s.left.left, s.left.lchild, s, s.left.lsibl, s.left)
    else
        throw("Invalid action: $(act).")
    end
end

# check if buffer is empty
bufferisempty(s::State) = s.right > length(s.tokens)

# check if can perform Reduce, that is, length(stack) >= 2
reducible(s::State) = s.left != nothing

# check if the state is final state
isfinal(s::State) = bufferisempty(s) && !reducible(s)

function expand_reduce_state!(res::Vector{State{Labeled}}, s::State{Labeled}, act::Function)
    for label in 1:length(s.parser.labels)
        push!(res, expand(s, act(label)))
    end
end

function expand_reduce_state!(res::Vector{State{Unlabeled}}, s::State{Unlabeled}, act::Function)
    push!(res, expand(s, act()))
end

function expandpred{T}(s::State{T})
    isfinal(s) && return []
    res = State{T}[]
    if reducible(s)
        expand_reduce_state!(res, s, reducer)
        if s.left.top != 0
            expand_reduce_state!(res, s, reducel)
        end
    end
    bufferisempty(s) || push!(res, expand(s, SHIFT))
    res
end

function expand_reduce_state(s::State{Unlabeled}, act::Function)
    expand(s, act())
end

function expand_reduce_state(s::State{Labeled}, act::Function)
    child = act == reducel ? s.left.top :
            act == reducer ? s.top :
            throw("Action must be reduce")
    expand(s, act(tokenat(s, child).label))
end

function expandgold(s::State)
    isfinal(s) && return []
    if !reducible(s)
        return [expand(s, SHIFT)]
    else
        s0 = tokenat(s, s.top)
        s1 = tokenat(s, s.left)
        if s1.head == s.top
            return [expand_reduce_state(s, reducel)]
        elseif s0.head == s.left.top
            if all(i -> tokenat(s, i).head != s.top, s.right:length(s.tokens))
                return [expand_reduce_state(s, reducer)]
            end
        end
    end
    return [expand(s, SHIFT)]
end

###################################################
################## feature function ###############
###################################################

function featuregen(s::State)
    # n0  = tokenat(s, s.right)
    # n1  = tokenat(s, s.right+1)
    # s0  = tokenat(s, s.top)
    # s0l = tokenat(s, :lchild)
    # s0r = tokenat(s, :rchild)
    # s1 = tokenat(s, :left, :top)
    # s2 = tokenat(s, :left, :left)
    # s1l = tokenat(s, :left, :lchild)
    # s1r = tokenat(s, :left, :rchild)

    n0  = tokenat(s, s.right)
    n1  = tokenat(s, s.right+1)
    s0  = tokenat(s, s.top)
    s0l = tokenat(s, s.lchild)
    s0r = tokenat(s, s.rchild)
    if s.left == nothing
        s1, s2, s1l, s1r = rootword, rootword, rootword, rootword
    else
        s1 = tokenat(s, s.left.top)
        s2 = tokenat(s, s.left.left)
        s1l = tokenat(s, s.left.lchild)
        s1r = tokenat(s, s.left.rchild)
    end
    len = size(s.parser.model.weights, 1) # used in @template macro
    @template begin
        # template (1)
        (s0.word,)
        (s0.tag,)
        (s0.word, s0.tag)
        (s1.word,)
        (s1.tag,)
        (s1.word, s1.tag)
        (n0.word,)
        (n0.tag,)
        (n0.word, n0.tag)

        # additional for (1)
        (n1.word,)
        (n1.tag,)
        (n1.word, n1.tag)

        # template (2)
        (s0.word, s1.word)
        (s0.tag, s1.tag)
        (s0.tag, n0.tag)
        (s0.word, s0.tag, s1.tag)
        (s0.tag, s1.word, s1.tag)
        (s0.word, s1.word, s1.tag)
        (s0.word, s0.tag, s1.tag)
        (s0.word, s0.tag, s1.word, s1.tag)

        # additional for (2)
        (s0.tag, s1.word)
        (s0.word, s1.tag)
        (s0.word, n0.word)
        (s0.word, n0.tag)
        (s0.tag, n0.word)
        (s1.word, n0.word)
        (s1.tag, n0.word)
        (s1.word, n0.tag)
        (s1.tag, n0.tag)

        # template (3)
        (s0.tag, n0.tag, n1.tag)
        (s1.tag, s0.tag, n0.tag)
        (s0.word, n0.tag, n1.tag)
        (s1.tag, s0.word, n0.tag)

        # template (4)
        (s1.tag, s1l.tag, s0.tag)
        (s1.tag, s1r.tag, s0.tag)
        (s1.tag, s0.tag, s0r.tag)
        (s1.tag, s1l.tag, s0.tag)
        (s1.tag, s1r.tag, s0.word)
        (s1.tag, s0.word, s0l.tag)

        # template (5)
        (s2.tag, s1.tag, s0.tag)
    end
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
#####
# ListNode
#####

Base.@kwdef mutable struct ListNode{V}
    val::V = 0
    next::Union{Nothing,ListNode{V}} = nothing
end

ListNode(x) = ListNode(x, nothing)

function ListNode{V}(x::Vector{V}) where {V}
    head = cur = ListNode(x[1])
    for i in 2:length(x)
        cur.next = ListNode(x[i])
        cur = cur.next
    end
    return head
end

val(x::ListNode) = x.val
val!(x::ListNode, y) = x.val = y
next(x::ListNode) = x.next
next!(x::ListNode, y) = x.next = y

Base.:(==)(x::ListNode, y::ListNode) = x.val == y.val && x.next == y.next

#####
# TreeNode
#####

Base.@kwdef mutable struct TreeNode{V}
    val::V = 0
    left::Union{Nothing,TreeNode} = nothing
    right::Union{Nothing,TreeNode} = nothing
end

TreeNode(x) = TreeNode(; val=x)
TreeNode(x, left) = TreeNode(; val=x, left=left)

function Base.:(==)(t1::TreeNode, t2::TreeNode)
    return t1.val == t2.val && t1.left == t2.left && t1.right == t2.right
end

function Base.convert(::Type{TreeNode{V}}, xs::Vector) where {V}
    @assert length(xs) > 0
    root = TreeNode(xs[1])
    _build_tree!(root, xs, 1)
    return root
end

function _build_tree!(t::TreeNode, xs::Vector, i::Int=1)
    n = length(xs)
    i_left = i * 2
    if i_left <= n && !isnothing(xs[i_left])
        t.left = TreeNode(xs[i_left])
        _build_tree!(t.left, xs, i_left)
    else
        t.left = nothing
    end
    i_right = i * 2 + 1
    if i_right <= n && !isnothing(xs[i_right])
        t.right = TreeNode(xs[i_right])
        _build_tree!(t.right, xs, i_right)
    else
        t.right = nothing
    end
end


text:
using HorizonSideRobots

function cross_marking(r::Robot)
    for i = 0:3
        mark_n_back(r, HorizonSide(i))
    end
    putmarker!(r)
end

function mark_n_back(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
        putmarker!(r)
    end
    while ismarker(r) == true
        move!(r, HorizonSide(mod(Int(side) + 2, 4)))
    end
end


text:
using Test
@time @testset "Utilities" begin
    include("test_utilities.jl")
end
@time @testset "Shell resultants" begin
    include("test_resultants.jl")
end
@time @testset "Composite shell dynamics" begin
    include("test_composite_shell_dynamics.jl")
end
@time @testset "Composite shell statics" begin
    include("test_composite_shell_statics.jl")
end
@time @testset "Shell statics" begin
    include("test_shell_statics.jl")
end
@time @testset "Beam section" begin
    include("test_section.jl")
end
@time @testset "Beam mesh" begin
    include("test_mesh.jl")
end
@time @testset "Beam modal" begin
    include("test_modal.jl")
end
@time @testset "Beam buckling" begin
    include("test_buckling.jl")
end
@time @testset "Beam transient dynamics" begin
    include("test_dyn.jl")
end
@time @testset "Shell dynamics" begin
    include("test_shell_dynamics.jl")
end
@time @testset "Composite Layup" begin
    include("test_composite_layup.jl")
end


text:
uq1 = UQ1D(4, 8, -1.0, 1.0, "uniform", "collocation")
uq2 = UQ2D(4, 8, [-1.0, 1.0, -1.0, 1.0], ["uniform", "uniform"], "collocation")

pce1 = ran_chaos(rand(uq1.op.quad.Nquad), uq1)
pce2 = ran_chaos(rand(uq2.nq), uq2)

chaos_ran(pce1, uq1)
chaos_ran(pce2, uq2)

lambda_tchaos(pce1, 1.0, uq1)
lambda_tchaos(pce2, 1.0, uq2)


text:
module WNNMDenoise

using Base.Iterators
using LinearAlgebra
using OffsetArrays
using ImageCore
using ImageCore.MappedArrays
using ImageCore: NumberLike, GenericGrayImage, GenericImage
using ImageFiltering
using ImageDistances
using Statistics
using LowRankApprox

using ImageQualityIndexes # TODO: remove this
using ProgressMeter # TODO: remove this

include("utilities.jl")

export WNNM

"""

# References

[1] Gu, S., Zhang, L., Zuo, W., & Feng, X. (2014). Weighted nuclear norm minimization with application to image denoising. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (pp. 2862-2869).

"""
struct WNNM
    "Estimated gaussian noise level"
    noise_level::Float64
    "Number of WNNM iterations"
    K::Int
    "step value for each WNNM iteration"
    δ::Float64
    "number of patches in each WNNM iteration"
    num_patches::Vector{Int}
    "patch size in each WNNM iteration"
    patch_size::Vector{Int}
    "patch stride in each WNNM iteration"
    patch_stride::Vector{Int}
    "weight constant used to estimate the remained noise level of each patch"
    λ::Float64
    "weight constant in WNNM solver"
    C::Float64
    "non-local search size in block matching"
    window_size::Int
    "matrix rank that used to allow early-stop in svd approximation"
    svd_rank::Vector{Int}
end

function WNNM(noise_level;
              δ=0.1,
              C=2sqrt(2),
              window_size=60,
              patch_size=nothing,
              num_patches=nothing,
              K=nothing,
              λ=nothing,
              svd_rank=nothing,
              patch_stride=nothing)
    if noise_level <= 20
        isnothing(patch_size) && (patch_size = 6)
        isnothing(num_patches) && (num_patches = 70)
        isnothing(K) && (K = 8)
        isnothing(λ) && (λ = 0.56)
    elseif noise_level <= 40
        isnothing(patch_size) && (patch_size = 7)
        isnothing(num_patches) && (num_patches = 90)
        isnothing(K) && (K = 12)
        isnothing(λ) && (λ = 0.56)
    elseif noise_level <= 60
        isnothing(patch_size) && (patch_size = 8)
        isnothing(num_patches) && (num_patches = 120)
        isnothing(K) && (K = 14)
        isnothing(λ) && (λ = 0.58)
    else
        isnothing(patch_size) && (patch_size = 9)
        isnothing(num_patches) && (num_patches = 140)
        isnothing(K) && (K = 14)
        isnothing(λ) && (λ = 0.58)
    end

    patch_size isa Number && (patch_size = fill(patch_size, K))
    isnothing(patch_stride) && (patch_stride = @. max(1, patch_size ÷ 2 - 1))
    patch_stride isa Number && (patch_stride = fill(patch_stride, K))

    if isnothing(svd_rank)
        svd_rank = map(enumerate(patch_size)) do (i, sz)
            min(10 + 5i, 10 + sz^2 ÷ 2)
        end
    elseif svd_rank isa Number
        svd_rank = fill(svd_rank, K)
    end

    num_patches = fill(num_patches - 10, K)
    drop_freq = 2
    for k in 2:K
        # drop by 10 for every 2 iteration
        num_patches[k] = (k - 1) % drop_freq == 0 ? num_patches[k - 1] - 10 : num_patches[k - 1]
    end

    WNNM(noise_level, K, δ, num_patches, patch_size, patch_stride, λ, C, window_size, svd_rank)
end

## Implementation

function (f::WNNM)(imgₑₛₜ, imgₙ; clean_img=nothing)
    if imgₑₛₜ === imgₙ
        imgₑₛₜ = copy(imgₙ)
    else
        copyto!(imgₑₛₜ, imgₙ)
    end

    T = eltype(eltype(imgₑₛₜ))
    for iter in 1:f.K
        @. imgₑₛₜ = imgₑₛₜ + f.δ * (imgₙ - imgₑₛₜ) # This iteration can be done more sophisticatedly

        # The noise level for the first iteration is known (whether it is estimated outside or a
        # white noise). The noise is removed in each iteration, so we have to estimate a noise level
        # at a patch level; the denoising performance on each patch can be different, which means a
        # global noise level can be misleading.
        σₚ = iter == 1 ? f.noise_level : nothing

        # calculating svd using blas threads is not an optimal parallel strategy
        imgₑₛₜ .= with_blas_threads(1) do
            _estimate_img(imgₑₛₜ, imgₙ;
                noise_level=f.noise_level,
                patch_size=f.patch_size[iter],
                patch_stride=f.patch_stride[iter],
                num_patches=f.num_patches[iter],
                window_size=f.window_size,
                λ=T(f.λ),
                C=T(f.C),
                σₚ=σₚ,
                svd_rank=f.svd_rank[iter],
            )
        end

        # TODO: remove this logging part when it is ready
        if !isnothing(clean_img)
            @info "Result" iter psnr = assess_psnr(clean_img, imgₑₛₜ, 255) num_patches = f.num_patches[iter] svd_rank=f.svd_rank[iter]
            display(Gray.(imgₑₛₜ ./ 255))
            sleep(0.1)
        end
    end
    return imgₑₛₜ
end

function _estimate_img(imgₑₛₜ, imgₙ; patch_size, patch_stride, num_patches, kwargs...)
    patch_size = ntuple(_ -> patch_size, ndims(imgₑₛₜ))

    # only set key patch stride in one dimension, this gives better denoising performance (Why?)
    # This also follow the original implementation
    Δ = CartesianIndex(patch_stride, ntuple(_->1, ndims(imgₑₛₜ)-1)...)
    # Δ = CartesianIndex(ntuple(_->patch_stride, ndims(imgₑₛₜ)))

    r = CartesianIndex(patch_size .÷ 2)
    R = CartesianIndices(imgₑₛₜ)
    R = first(R) +r:Δ:last(R) -r

    imgₑₛₜ⁺ = zeros(eltype(imgₑₛₜ), axes(imgₑₛₜ))
    W = zeros(Int, axes(imgₑₛₜ))

    progress = Progress(length(R[1:patch_stride:end]))
    out_patches = [Matrix{eltype(imgₑₛₜ)}(undef, prod(patch_size), num_patches) for i in 1:Threads.nthreads()]
    Threads.@threads for p in R
        out = out_patches[Threads.threadid()]
        fill!(out, zero(eltype(out)))
        patch_q_indices = _estimate_patch!(out, imgₑₛₜ, imgₙ, p; patch_size=patch_size, num_patches=num_patches, kwargs...)

        view(W, patch_q_indices) .+= 1
        view(imgₑₛₜ⁺, patch_q_indices) .+= out
        next!(progress)
    end

    return imgₑₛₜ⁺ ./ max.(W, 1)
end

function _estimate_patch!(out, imgₑₛₜ, imgₙ, p;
                         noise_level,
                         patch_size::Tuple,
                         num_patches::Integer,
                         window_size,
                         λ,
                         C,
                         svd_rank,
                         σₚ=nothing)
    rₚ = CartesianIndex(patch_size .÷ 2)
    p_indices = p - rₚ:p + rₚ

    patch_q_indices = block_matching(imgₑₛₜ, p;
        num_patches=num_patches,
        patch_size=patch_size,
        search_window_size=window_size,
        patch_stride=1
    )
    patch_q_indices = hcat([indices[:] for indices in patch_q_indices]...)
    m = mean(@view(imgₑₛₜ[patch_q_indices]); dims=2)

    if isnothing(σₚ)
        # Try: use the mean estimated σₚ of each patch
        σₚ = _estimate_noise_level(view(imgₑₛₜ, p_indices), view(imgₙ, p_indices), noise_level; λ=λ)
    end
    out .= @view(imgₑₛₜ[patch_q_indices]) .- m
    WNNM_optimizer!(out, out, eltype(out)(σₚ); C=C, svd_rank=svd_rank)
    out .+= m

    return patch_q_indices
end

@doc raw"""
    WNNM_optimizer(Y, σₚ; C, rank, fixed_point_num_iters=3)

Optimizes the weighted nuclear norm minimization problem with a fixed point estimation

```math
    \min_X \lVert Y - X \rVert^2_{F} + \lVert X \rVert_{w, *}
```

The weight `w` is specially chosen so that it satisfies the condition of Corollary 1 in [1].

# References

[1] Gu, S., Zhang, L., Zuo, W., & Feng, X. (2014). Weighted nuclear norm minimization with application to image denoising. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (pp. 2862-2869).

"""
function WNNM_optimizer!(out, Y, σₚ; C, svd_rank, fixed_point_num_iters=3)
    # Apply Corollary 1 in [1] for image denoise purpose
    # Note: this solver is reserved to the denoising method and is not supposed to be used in other
    #       applications; it simply isn't designed so.

    # This is different from the original implementation. Here we use an approximate version of svd;
    # it gives better performance in both speed and denoising result.
    n = size(Y, 2)
    if svd_rank >= ceil(Int, 0.5minimum(size(Y)))
        U, ΣY, V = svd!(Y)
    else
        U, ΣY, V = LowRankApprox.psvd(Y; rank=svd_rank)
    end

    # For image denoising problems, it is natural to shrink large singular value less, i.e., to set
    # smaller weight to large singular value. For this reason, it uses `w = (C * sqrt(n))/(ΣX + eps())`
    # as the weights; inversely propotional to ΣX. With singular values ΣX sorted ascendingly, the
    # condition for Corollary 1 holds, and thus we could directly get the desired solution in a single
    # step.

    # Here we iterate more than once because we don't know what ΣX is; we have to iterate it a while
    # from ΣY to get an relatively good estimation of it.
    # TODO: could we set default σₚ as 0?
    # TODO: is this the best initialization we can get?
    ΣX = @. sqrt(max(ΣY^2 - n * σₚ^2, 0))
    for _ in 1:fixed_point_num_iters
        # the iterative algorithm proposed in section 2.2.2 in [1]
        # Step 1 in the iterative algorithm becomes trivial and a no-op

        # Step 2 degenerates to a soft thresholding; both P and Q are identity matrix.
        # all in one line to avoid unnecessary allocation for temporarily variable w
        @. ΣX = soft_threshold(ΣY, (C * sqrt(n) * σₚ^2) / (ΣX + eps()))
    end

    mul!(out, rmul!(U, Diagonal(ΣX)), V')
end

function _estimate_noise_level(patchₑₛₜ, patchₙ, σₙ; λ=0.56)
    # Estimate the noise level of given patch during the WNNM iteration
    # we still need to know the input noisy level σₙ to give an estimation
    λ * sqrt(abs(σₙ^2 - mse(patchₑₛₜ, patchₙ)))
end

end


text:
function starterdeckreader(dualcolors)::Vector{Int}
    monos = ["W","U","B","R","G"] .* ".txt"
    duals = dualcolors .* ".txt"
    files = vcat(monos, duals)
    dir = joinpath(@__DIR__, "..", "data", "starterdecks")
    decklists = deckreader_mtga_format.( joinpath.(dir, files) )
    starterdecks = deckinfo.(decklists)
    col = array_rep_all(card_db)
    for deck in starterdecks
        for elem in deck
            ix = first(elem.prints).index
            col[ix] += elem.amount
            if col[ix] > 4
                col[ix] = 4
            end
        end
    end
    col
end

const startersequences = Vector{Int}[
   starterdeckreader( ["UG", "WR", "BR", "GW", "UB"]),
   starterdeckreader( ["BR", "GW", "BW", "UG", "UR"]),
   starterdeckreader( ["BW", "UR", "RG", "UW", "GB"] ),
   starterdeckreader( ["RG", "GB", "UB", "UW", "WR"] )
]

const tenstarters = starterdeckreader( ["RG", "GB", "UB", "UW", "WR", "BW", "UR", "BR", "GW", "UG"])

const basiccards = let
    dir = joinpath(@__DIR__, "..", "data", "starterdecks")
    files = readdir(dir)
    filter!(x -> in(x, ["W.txt", "U.txt", "B.txt", "R.txt", "G.txt"]), files)
    decklists = deckreader_mtga_format.( joinpath.(dir, files) )
    starterdecks = deckinfo.(decklists)
    col = array_rep_all(card_db)
    for deck in starterdecks
        for elem in deck
            ix = first(elem.prints).index
            col[ix] += elem.amount
            if col[ix] > 4
                col[ix] = 4
            end
        end
    end
    col
end

const kldgrant = let
    file = joinpath(@__DIR__, "..", "data", "KLD_starter", "kld_starter.txt" )
    deck = deckreader_mtga_format( file )
    deck = deckinfo(deck)

    col = array_rep_all(card_db)
    for elem in deck
        ix = first(elem.prints).index
        col[ix] += elem.amount
    end

    for r in 1:4
        for s in [:KLD, :AER]
            for ix in indices_set_rarity(card_db, s, r)
                col[ix] += r > 3 ? 1 : 2
            end
        end
    end

    for i in eachindex(col)
        if col[i] > 4
            col[i] = 4
        end
    end
    col

end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
# TODO eventually wrap this in some UI to ask the questions
function create_registry(name::String; uuid = Base.Random.uuid4(), depot = DEFAULT_DEPOT[],
                         description::String = "", repo::String = "")

    path = joinpath(depot, "registries", name)
    if !isdir(path)
        mkdir(path)
    end

    if isfile(joinpath(path, "registry.toml"))
        error("Registry already exists. Aborting.")
    end

    # TODO learn the TOML tools
    open(joinpath(path, "registry.toml"), "w") do io
        write(io, "name = \"$name\"\n")
        write(io, "uuid = \"$uuid\"\n")
        write(io, "repo = \"$repo\"\n")
        write(io, "\n")
        write(io, "description = \"\"\"\n")
        write(io, description)
        write(io, "\"\"\"")
        write(io, "\n")
        write(io, "[packages]\n")
    end

    # TODO make this directory into a git repo (and try push upstream?)
end

"""
    registry_depot(registry::String)

Returns the path of the depot that registry named `registry` is in.
"""
function registry_depot(registry::String)
    for depot in Pkg3.depots()
        registry_path = joinpath(depot, "registries", registry)
        if isfile(joinpath(registry_path, "registry.toml"))
            return depot
        end
    end
    error("Unknown registry \"$registry\"")
end

"""
registry_depot(registry::String)

Returns the path containing the registry.toml file that the registry named `registry` is in.
"""
registry_path(registry::String) = joinpath(registry_depo(registry), "registries", registry)

# TODO clone a registry
function clone_registry(repo::String; depot = DEFAULT_DEPOT[])
    # Some LibGit2 stuff
end

# TODO push a registry to its upstream repo / perhaps make a PR for Uncurated on github


text:
module RadonAlgorithm

Base.Experimental.@optlevel 3

export Radon, RadonSquare
export AbstractFBP
export FBP, FBPFFTSquare, FBPAFFT, FBPAFFTSquare
export RadonInfo, FBPInfo

using ..Applicative
using ..Monads
using ..CTImages
using ..Geometry
using ..FanBeam: fan2para, para2fan
import ..Utils: ORI, linspace, _atype, half, width
import ..AbstractAlgorithms:
    radon,
    iradon,
    project_image,
    reconstruct_image,
    _alg_progress
using ..AbstractAlgorithms, ..Coordinates
using ..Interpolation: AbstractInterp2DOrNone, interpolate
using FFTW, ProgressMeter, IntervalSets

function _radon end
function _iradon end
function _alg_method end

@inline function _radon(
    a::AbstractProjectionAlgorithm,
    img::AbstractMatrix,
    args...;
    kwargs...
)
    _radon(_alg_method(a), a, img, args...; kwargs...)
end

@inline function _iradon(
    a::AbstractIRadonAlgorithm,
    sinog::AbstractMatrix,
    args...;
    kwargs...
)
    _iradon(_alg_method(a), a, sinog, args...; kwargs...)
end

macro _defradonfn(f::Symbol, body)
    esc(quote
        @inline function $f(
            image::AbstractMatrix{T},
            ts::AbstractVector{X},
            ϕs::AbstractVector{Y};
            ν::Real = 1,
            scale::Optional{Real} = nothing,
            τ::Optional{Real} = nothing,
            ratio::Optional{Real} = nothing,
            background::Optional{Z} = nothing,
            rescaled::Bool = false,
            interpolation::Optional{Interp} = nothing,
            progress::Bool = false,
        ) where {
            T <: Real,
            X <: Real,
            Y <: Real,
            Z <: Real,
            Interp <: AbstractInterp2DOrNone,
        }
            ν = maybe(ν, scale)
            @assert ν > 0
            rows, cols = size(image)
            τ = maybe(maybe(rows / cols, τ), ratio)
            @assert τ > 0
            nd = length(ts)
            nϕ = length(ϕs)
            scϕs = sincos.(ϕs)
            z::T = maybe(zero(T), background)
            sinog = similar(_atype(image), nd, nϕ)
            fill!(sinog, z)
            rimage = rescaled ? rescale(image) : image
            interp = isnothing(interpolation) ?
                interpolate(rimage) : interpolation(rimage)
            CTSinogram($body)
        end
    end)
end

macro _defiradonfn(f::Symbol, body)
    esc(quote
        @inline function $f(
            sinog::AbstractMatrix{T},
            xs::AbstractVector{U1},
            ys::AbstractVector{U2},
            ::Cartesian;
            ν::Real = 1,
            scale::Optional{Real} = nothing,
            ϕs::Optional{I} = nothing,
            angles::Optional{J} = nothing,
            background::Optional{U3} = nothing,
            filter::Optional{F} = nothing,
            interpolation::Optional{Interp} = nothing,
            progress::Bool = false,
        ) where {
            T <: Real,
            U1 <: Real,
            U2 <: Real,
            I <: Interval{:closed},
            J <: Interval{:closed},
            U3 <: Real,
            F <: AbstractCTFilter,
            Interp <: AbstractInterp2DOrNone,
        }
            ν = maybe(ν, scale)
            @assert ν > 0
            ϕs = maybe(ϕs, angles)
            nd, nϕ = size(sinog)
            cols = length(xs)
            rows = length(ys)
            filtered = apply(maybe(RamLak(), filter)) do f
                filter_freq = fft(sinog, 1) .* f(T, nd, nϕ)
                bfft(filter_freq, 1) |> real
            end
            interpolation = maybe(interpolate, interpolation)
            interp = interpolation(filtered)
            t₀::T = (nd + 1) / 2
            ϕs = maybe(ORI(0..2π), ϕs)
            scϕs = sincos.(linspace(T, ϕs, nϕ))
            z::T = maybe(zero(T), background)
            tomog = similar(_atype(sinog), rows, cols)
            fill!(tomog, z)
            CTTomogram($body)
        end
        @inline function $f(
            sinog::AbstractMatrix{T},
            xs::AbstractVector{X},
            ys::AbstractVector{Y};
            kwargs...
        ) where {T<:Real,X<:Real,Y<:Real}
            $f(sinog, xs, ys, Cartesian(); kwargs...)
        end
    end)
end

macro _defradonalgfn(A::Symbol, f::Symbol)
    esc(quote
        @inline function (a::$A)(
            image::AbstractMatrix,
            ts::AbstractVector,
            ϕs::AbstractVector;
            kwargs...
        )
            $f(image, ts, ϕs; kwargs...)
        end
        @inline function (a::$A)(
            image::AbstractMatrix,
            ts::Interval,
            ϕs::Interval;
            kwargs...
        )
            _radon(a, image, ts, ϕs; kwargs...)
        end
        @inline function (a::$A)(image::AbstractMatrix; kwargs...)
            _radon(a, image; kwargs...)
        end
    end)
end

macro _defiradonalgfn(A::Symbol, f::Symbol)
    esc(quote
        @inline function (a::$A)(
            sinog::AbstractMatrix,
            xs::AbstractVector,
            ys::AbstractVector,
            coo::AbstractCoordinates = Cartesian();
            kwargs...
        )
            $f(sinog, xs, ys, coo; kwargs...)
        end
        @inline function (a::$A)(
            sinog::AbstractMatrix,
            xs::Interval,
            ys::Interval,
            coo::AbstractCoordinates = Cartesian();
            kwargs...
        )
            _iradon(a, sinog, xs, ys, coo; kwargs...)
        end
        @inline function (a::$A)(
            sinog::AbstractMatrix,
            xs::AbstractVector,
            coo::AbstractCoordinates = Cartesian();
            kwargs...
        )
            $f(sinog, xs, xs, coo; kwargs...)
        end
        @inline function (a::$A)(sinog::AbstractMatrix; kwargs...)
            _iradon(a, sinog; kwargs...)
        end
    end)
end

function _radon_progress(n::Integer, p::Bool, dt::Real=0.2)
    _alg_progress(Progress, "Computing Radon transform...", n, p, dt)
end

function _iradon_progress(n::Integer, p::Bool, dt::Real=0.2)
    _alg_progress(Progress, "Computing inverse Radon transform...", n, p, dt)
end

abstract type AbstractFBP <: AbstractIRadonAlgorithm end

include("Filters.jl")
include("radon/radon.jl")
include("iradon/iradon.jl")

struct RadonInfo{
    G <: AbstractGeometry,
    A <: AbstractProjectionAlgorithm
} <: AlgorithmInfo{A}
    geometry::G
    algorithm::A
end

RadonInfo(g::AbstractGeometry) = RadonInfo(g, Radon())

struct FBPInfo{G<:AbstractGeometry,A<:AbstractFBP} <: AlgorithmInfo{A}
    geometry::G
    algorithm::A
end

FBPInfo(g::AbstractGeometry) = FBPInfo(g, FBP())

end # module


text:
"""
    NormInfinityBridge{T}

The `NormInfinityCone` is representable with LP constraints, since
``t \\ge \\max_i \\lvert x_i \\rvert`` if and only if
``t \\ge x_i`` and ``t \\ge -x_i`` for all ``i``.
"""
struct NormInfinityBridge{T,F,G} <:
       SetMapBridge{T,MOI.Nonnegatives,MOI.NormInfinityCone,F,G}
    constraint::CI{F,MOI.Nonnegatives}
end
function concrete_bridge_type(
    ::Type{<:NormInfinityBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormInfinityCone},
) where {T}
    F = MOIU.promote_operation(+, T, G, G)
    return NormInfinityBridge{T,F,G}
end

function map_set(::Type{<:NormInfinityBridge}, set::MOI.NormInfinityCone)
    return MOI.Nonnegatives(2 * (MOI.dimension(set) - 1))
end
function inverse_map_set(::Type{<:NormInfinityBridge}, set::MOI.Nonnegatives)
    return MOI.NormInfinityCone(div(MOI.dimension(set), 2) + 1)
end

function map_function(::Type{<:NormInfinityBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    t = scalars[1]
    lb = scalars[2:end]
    ub = MOIU.operate(-, T, lb)
    f_new = MOIU.operate(vcat, T, ub, lb)
    for i in 1:(2*(length(scalars)-1))
        MOIU.operate_output_index!(+, T, i, f_new, t)
    end
    return f_new
end
function inverse_map_function(::Type{<:NormInfinityBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    t = MOIU.operate!(/, T, sum(scalars), T(length(scalars)))
    d = div(length(scalars), 2)
    x = MOIU.operate!(
        /,
        T,
        MOIU.operate!(-, T, scalars[(d+1):end], scalars[1:d]),
        T(2),
    )
    return MOIU.operate(vcat, T, t, x)
end
# Given a_i is dual on t - x_i >= 0 and b_i is dual on t + x_i >= 0,
# the dual on (t, x) in NormInfinityCone is (u, v) in NormOneCone, where
# v_i = -a_i + b_i and u = sum(a) + sum(b).
function adjoint_map_function(::Type{<:NormInfinityBridge}, func)
    scalars = MOIU.eachscalar(func)
    t = sum(scalars)
    d = div(length(scalars), 2)
    x = (scalars[(d+1):end] - scalars[1:d])
    return vcat(t, x)
end
function inverse_adjoint_map_function(
    ::Type{<:NormInfinityBridge{T}},
    func::AbstractVector{T},
) where {T}
    # This is used by `MOI.ConstraintDualStart`.
    # The result should belong to `MOI.Nonnegatives` and the sum of the elements should
    # be `t`.
    # Only one of dual of the upper and lower bound should be active so one of the two
    # duals is zero. We know which one since they should be nonnegative.
    # Then if `t = sum abs(x_i)`, we will indeed have only one of them being zero.
    t = func[1]
    y = func[2:end]
    lb = [x > 0 ? x : zero(x) for x in y]
    ub = [x < 0 ? -x : zero(x) for x in y]
    x = [ub; lb]
    return x .+ (t - sum(x)) / length(x)
end

"""
    NormOneBridge{T}

The `NormOneCone` is representable with LP constraints, since
``t \\ge \\sum_i \\lvert x_i \\rvert`` if and only if there exists a vector y such that
``t \\ge \\sum_i y_i`` and ``y_i \\ge x_i``, ``y_i \\ge -x_i`` for all ``i``.
"""
struct NormOneBridge{T,F,G} <: AbstractBridge
    y::Vector{MOI.VariableIndex}
    nn_index::CI{F,MOI.Nonnegatives}
end
function bridge_constraint(
    ::Type{NormOneBridge{T,F,G}},
    model::MOI.ModelLike,
    f::G,
    s::MOI.NormOneCone,
) where {T,F,G}
    f_scalars = MOIU.eachscalar(f)
    d = MOI.dimension(s)
    y = MOI.add_variables(model, d - 1)
    ge = MOIU.operate(-, T, f_scalars[1], MOIU.operate(sum, T, y))
    lb = f_scalars[2:d]
    ub = MOIU.operate(-, T, lb)
    lb = MOIU.operate!(+, T, lb, MOI.VectorOfVariables(y))
    ub = MOIU.operate!(+, T, ub, MOI.VectorOfVariables(y))
    f_new = MOIU.operate(vcat, T, ge, ub, lb)
    nn_index = MOI.add_constraint(model, f_new, MOI.Nonnegatives(2d - 1))
    return NormOneBridge{T,F,G}(y, nn_index)
end

function MOI.supports_constraint(
    ::Type{NormOneBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormOneCone},
) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:NormOneBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:NormOneBridge{T,F}}) where {T,F}
    return [(F, MOI.Nonnegatives)]
end
function concrete_bridge_type(
    ::Type{<:NormOneBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.NormOneCone},
) where {T}
    S = MOIU.scalar_type(G)
    F = MOIU.promote_operation(
        vcat,
        T,
        MOIU.promote_operation(+, T, S, S),
        MOIU.promote_operation(-, T, S, S),
    )
    return NormOneBridge{T,F,G}
end

# Attributes, Bridge acting as a model
MOI.get(b::NormOneBridge, ::MOI.NumberOfVariables) = length(b.y)
MOI.get(b::NormOneBridge, ::MOI.ListOfVariableIndices) = b.y
function MOI.get(
    b::NormOneBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.Nonnegatives},
) where {T,F}
    return 1
end
function MOI.get(
    b::NormOneBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.Nonnegatives},
) where {T,F}
    return [b.nn_index]
end

# References
function MOI.delete(model::MOI.ModelLike, c::NormOneBridge)
    MOI.delete(model, c.nn_index)
    return MOI.delete(model, c.y)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    c::NormOneBridge{T,F,G},
) where {T,F,G}
    nn_func =
        MOIU.eachscalar(MOI.get(model, MOI.ConstraintFunction(), c.nn_index))
    t = MOIU.operate!(/, T, nn_func[1] + sum(nn_func), T(2))
    d = div(length(nn_func) - 1, 2)
    x = MOIU.operate!(
        /,
        T,
        MOIU.operate!(-, T, nn_func[(d+2):end], nn_func[2:(d+1)]),
        T(2),
    )
    return MOIU.convert_approx(
        G,
        MOIU.remove_variable(MOIU.operate(vcat, T, t, x), c.y),
    )
end
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet, c::NormOneBridge)
    dim = div(
        MOI.dimension(MOI.get(model, MOI.ConstraintSet(), c.nn_index)) + 1,
        2,
    )
    return MOI.NormOneCone(dim)
end
function MOI.supports(
    ::MOI.ModelLike,
    ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:NormOneBridge},
)
    return true
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::NormOneBridge{T},
    value,
) where {T}
    x_value = value[1 .+ (1:length(bridge.y))]
    y_value = abs.(x_value)
    for i in eachindex(bridge.y)
        MOI.set(model, MOI.VariablePrimalStart(), bridge.y[i], y_value[i])
    end
    nn_value = vcat(
        value[1] - reduce(+, y_value, init = zero(T)),
        y_value - x_value,
        y_value + x_value,
    )
    MOI.set(model, attr, bridge.nn_index, nn_value)
    return
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::NormOneBridge,
)
    nn_primal = MOI.get(model, attr, bridge.nn_index)
    t = (nn_primal[1] + sum(nn_primal)) / 2
    d = length(bridge.y)
    x = (nn_primal[(d+2):end] - nn_primal[2:(d+1)]) / 2
    return vcat(t, x)
end
# Given a_i is dual on y_i - x_i >= 0 and b_i is dual on y_i + x_i >= 0 and c is dual on t - sum(y) >= 0,
# the dual on (t, x) in NormOneCone is (u, v) in NormInfinityCone, where
# v_i = -a_i + b_i and u = c.
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::NormOneBridge,
)
    nn_dual = MOI.get(model, attr, bridge.nn_index)
    d = length(bridge.y)
    x = nn_dual[(d+2):end] - nn_dual[2:(d+1)]
    return vcat(nn_dual[1], x)
end
# value[1 + i] = nn_dual[1 + d + i] - nn_dual[1 + i]
# and `nn_dual` is nonnegative. By complementarity slackness, only one of each
# `nn_dual` can be nonzero (except if `x = 0`) so we can set
# depending on the sense of `value[1 + i]`.
function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintDualStart,
    bridge::NormOneBridge,
    value,
)
    d = length(bridge.y)
    nn_dual = zeros(eltype(value), 2d + 1)
    nn_dual[1] = value[1]
    for i in eachindex(bridge.y)
        v = value[1+i]
        if v < 0
            nn_dual[1+i] = -v
        else
            nn_dual[1+d+i] = v
        end
    end
    MOI.set(model, MOI.ConstraintDualStart(), bridge.nn_index, nn_dual)
    return
end


text:
###################################################
## rect.jl
## Rectangular virtual city
###################################################

"""
  `rectNetwork` : n*m points rectangular network
 - all road same type (5)
 - `distX` and `distY` are the width and height of a block (in meters)
"""
function rectNetwork(width::Int = 5, height::Int = 10; distX::Float64=200., distY::Float64 = 200.)
        #nodes are indexed as follow :
        #  12
        #  34 (for width=2,2)

        function coordToLoc(i,j) # matrix indexation, i= row j=column
            return j + (i-1)*width
        end
        function locToCoord(n)
            return (div((n-1),width) + 1, ((n-1) % width) + 1)
        end

        #Create nodes
        nodes = Array{Node}(width*height)
        for i in 1:height, j in 1:width
            nodes[coordToLoc(i,j)] = Node((j-(width+1)/2)*distX, ((height+1)/2 - i)*distY)
        end
        roads = Dict{Tuple{Int,Int},Road}()

        #Vertical roads
        for i in 1:(height-1), j in 1:width
            n1 = coordToLoc(i,j); n2 = coordToLoc(i+1,j)
            dist = distanceCoord(nodes[n1],nodes[n2])
            roads[(n1,n2)] = Road(n1,n2,dist, 5)
            roads[(n2,n1)] = Road(n2,n1,dist, 5)
        end

        #Horizontal roads
        for i in 1:height, j in 1:width - 1
            n1 = coordToLoc(i,j); n2 = coordToLoc(i,j+1)
            dist = distanceCoord(nodes[n1],nodes[n2])
            roads[(n1,n2)] = Road(n1,n2,dist, 5)
            roads[(n2,n1)] = Road(n2,n1,dist, 5)
        end
        return Network(nodes, roads)
end

squareNetwork(width::Int = 5; distance::Float64=200.) = rectNetwork(width, width, distX=distance, distY=distance)

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module CoinbasePro

using HTTP, JSON, DataFrames, Parameters, Dates, Base64, SHA, Printf, Decimals

const HOST = "https://api.pro.coinbase.com"
_dropmissing(x) = x
_dropmissing(d::Dict) = Dict(k => v for (k, v) in d if !ismissing(v))

struct AuthenticatedClient
    api_key::String
    secret_key::String
    passphrase::String
end

function auth_headers(c::AuthenticatedClient, method::String, endpoint::String, body::String)
    # CB-ACCESS-KEY The api key as a string.
    # CB-ACCESS-SIGN The base64-encoded signature (see Signing a Message).
    # CB-ACCESS-TIMESTAMP A timestamp for your request.
    # CB-ACCESS-PASSPHRASE The passphrase you specified when creating the API key.
    timestamp = @sprintf "%f" datetime2unix(now(UTC))
    message = string(timestamp, method, endpoint, body)
    hmac_key = base64decode(c.secret_key)
    signature = hmac_sha256(hmac_key, message)
    return Dict(
        "CB-ACCESS-KEY" => c.api_key,
        "CB-ACCESS-SIGN" => base64encode(signature),
        "CB-ACCESS-TIMESTAMP" => timestamp,
        "CB-ACCESS-PASSPHRASE" => c.passphrase,
        "Content-Type" => "application/json",
    )
end

function _handle_response(res::HTTP.Response)
    if res.status >= 300
        @error "Bad request $r"
    end
    return JSON.parse(String(res.body))
end

function _request(c::AuthenticatedClient, method::String, endpoint::String, params::Dict=Dict(), data::Dict=Dict())
    url = string(HOST, endpoint)
    body = length(data) > 0 ? JSON.json(data) : ""
    headers = auth_headers(c, method, endpoint, body)
    _handle_response(
        HTTP.request(method, url, headers, body, query=params)
    )
end

function _request(method::String, endpoint::String, params::Dict=Dict(), data::Dict=Dict())
    url = string(HOST, endpoint)
    _handle_response(
        HTTP.request(method, url, ["Content-Type" => "application/json"], JSON.json(data), query=params)
    )
end

get_products() = _request("GET", "/products")
get_product_order_book(product_id::String, level::Int=1) =
    _request("GET", string("/products/", product_id, "/book"), Dict(:level => level))
get_product_ticker(product_id::String) = _request("GET", "/products/$product_id/ticker")

# function get_product_trades(product_id::String, before::String="",
#     after::String="",
#     limit::Union{missing,Int}=missing,
#     result::Union{Missing,Int}=missing,)
#     _request(c, "GET", "/products/$product_id/ticker")
# end

function get_product_historic_rates(
        product_id::String;
        start::Union{Missing,String}=missing,
        stop::Union{Missing,String}=missing,
        granularity::Union{Missing,Int}=missing,
    )
    params = Dict()
    !ismissing(start) && setindex!(params, start, :start)
    !ismissing(stop) && setindex!(params, stop, :end)
    if !ismissing(granularity)
        accepted_grans = [60, 300, 900, 3600, 21600, 86400]
        if !(granularity in accepted_grans)
            error("Specified granularity is $granularity, must be in approved values: $accepted_grans")
        end
        params[:granularity] = granularity
    end
    @show params
    res = _request("GET", "/products/$product_id/candles", params)
    rows = []
    for row in res
        push!(rows, (dt = Dates.unix2datetime(row[1]), low = row[2], high = row[3], open = row[4], close = row[5], volume = row[6]))
    end
    DataFrame(rows)
end

get_product_24hr_stats(product_id::String) = _request("GET", "/products/$product_id/stats")
get_currencies() = _request("GET", "/currencies")
get_time() = _request("GET", "/time")

## auth methods
get_orders(c::AuthenticatedClient; product_id::Union{Missing,String}=missing, order_id::Union{Missing,String}=missing) =
    _request(c, "GET", "/orders", _dropmissing(Dict("product_id" => product_id, "order_id" => order_id)))

get_accounts(c::AuthenticatedClient) = _request(c, "GET", "/accounts/")

convert_stablecoin(c::AuthenticatedClient, from::String, to::String, amount::Union{Integer,Decimals.Decimal}) =
    _request(c, "POST", "/conversions", Dict(), Dict("from" => from, "to" => to, "amount" => amount))

end # module