text:
module TestUpdatableQRFactorizations
using LinearAlgebra
using LinearAlgebra: Givens
using Test
using UpdatableQRFactorizations

@testset "UpdatableQRFactorizations" begin

    n, m = 32, 8 # NOTE: a test of UpdatableQR assumes n = 4m
    A = randn(n, m)
    F = GivensQR(A)
    @test F isa GivensQR

    # testing size
    @test size(F) == (n, m)
    for i in 1:3
        @test size(F, 1) == size(A, 1)
    end

    # testing copy
    CF = copy(F)
    @test CF.Q == F.Q
    @test CF.R == F.R
    @test size(CF) == size(F)

    @testset "GivensQ" begin
        Q = Matrix(F.Q)

        @test Q'Q ≈ I(n)
        @test Q*Q' ≈ I(n)

        @test F.Q'*F.Q == I(n)
        @test F.Q'*F.Q == I(n)

        x = randn(n)
        Qx = Q*x
        @test Qx ≈ F.Q*x
        @test F.Q \ Qx ≈ x
        Qtx = Q'*x
        @test Qtx ≈ F.Q'*x
        @test F.Q' \ Qtx ≈ x

        r = 4
        X = randn(r, n)
        @test Q*X' ≈ F.Q*X'
        @test X*Q ≈ X*F.Q
        @test Q'*X' ≈ F.Q'*X'
        @test X*Q' ≈ X*F.Q'
        @test Q \ (Q*X') ≈ X'
        @test Q' \ (Q'*X') ≈ X'

        x = randn(n)
        y = similar(x)
        X = randn(n, r)
        Y = similar(X)
        @test mul!(y, F.Q, x) ≈ Q*x
        @test mul!(Y, F.Q, X) ≈ Q*X
        @test mul!(y, F.Q', x) ≈ Q'*x
        @test mul!(Y, F.Q', X) ≈ Q'*X

        # GivensQ with empty rotations vector behaves like identity
        Qid = GivensQ(Vector{Givens{Float64}}(undef, 0), n, m)
        @. y = x
        @test lmul!(Qid, x) ≈ y
        @test lmul!(Qid', x) ≈ y
        @test Matrix(Qid) ≈ I(n)
    end

    @testset "GivensQR" begin
        QA = F.Q'*A
        @test UpperTriangular(QA[1:m, 1:m]) ≈ F.R
        @test Matrix(F) ≈ A

        x = randn(m)
        Ax = A*x
        @test F \ (A*x) ≈ x
        r = 4
        X = randn(m, r)
        @test F \ (A*X) ≈ X
    end

    @testset "UpdatableQR" begin
        F = UpdatableQR(A)
        @test Matrix(F) ≈ A

        x = randn(m)
        Ax = A*x
        @test F \ (A*x) ≈ x
        r = 4
        X = randn(m, r)
        @test F \ (A*X) ≈ X

        # testing column addition
        F = UpdatableQR(A)
        b = randn(n)
        add_column!(F, b)
        @test Matrix(F) ≈ [A b]

        # testing column addition at particular index
        x = randn(m+1)
        X = randn(m+1, r)
        for i in 1:m
            F = UpdatableQR(A)
            add_column!(F, b, i)
            B = [A[:, 1:i-1] b A[:, i:end]]
            @test B ≈ Matrix(F)
            Q = Matrix(F.Q)
            @test Q'Q ≈ I(n) # maintains orthogonality

            # testing solves
            @test F \ (B*x) ≈ x
            r = 4
            @test F \ (B*X) ≈ X
        end

        # testing column addition
        F = UpdatableQR(A)
        remove_column!(F)
        @test Matrix(F) ≈ A[:, 1:m-1]

        # testing column removal
        for i in 1:m
            F = UpdatableQR(A)
            remove_column!(F, i)
            B = @view A[:, 1:m .!= i]
            @test B ≈ Matrix(F)
            Q = Matrix(F.Q)
            @test Q'Q ≈ I(n) # maintains orthogonality
        end

        # testing empty initialization
        F = UpdatableQR(n, m)
        @test size(F) == (n, 0)
        @test Matrix(F) ≈ zeros(n, 0)

        order = reverse(1:m) # adding columns of A in reverse order (but to front of factorization)
        for i in order
            add_column!(F, A[:, i], 1)
        end
        @test Matrix(F) ≈ A

        F = UpdatableQR(n, m)
        add_column!(F, A) # adding all columns in A simultaneously
        @test Matrix(F) ≈ A

        # try out full rank QR
        F = UpdatableQR(n, n) # NOTE: assumes n = 4m
        add_column!(F, A)
        add_column!(F, A)
        add_column!(F, A)
        add_column!(F, A)
        M = Matrix(F)
        @test size(M) == (n, n)
        for i in 1:4
            Mi = @view M[:, m*(i-1)+1:m*(i-1)+m]
            @test Mi ≈ A
        end
    end

    # multiply is very efficient!
    # @time Q*x
    # @time Q*x
    # @time F.Q*x
    # @time F.Q*x

    # IDEA: testing overdetermined system, irrelevant for now
    # n, m = 128, 64
    # A = randn(m, n)
    # @time F = GivensQR(A)

    # A = randn(n, m)
    # rotations = allocate_rotations(A)
    # F = GivensQR()
    # for i in 1:m
    #     Ai = @view A[:, 1:i]
    #     append_column!()
    # end
end

end # module


text:
module TestCategoricalAlgebra

using Test

@testset "Permutations" begin
  include("Permutations.jl")
end

end

ext: jl


text:
include("neural_network.jl")
include("tabular.jl")
include("linear.jl")


text:
function floor(x::DoubleFloat{T}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    # !isinteger(LO(x)), LO(x) is mixed or fractional
    if isinteger(HI(x))
        if isfractional(LO(x))
            if signbit(LO(x))
                DoubleFloat(HI(x)-one(T), zero(T))
            else
                DoubleFloat(HI(x), zero(T))
            end
        else # LO(x) is mixed: +/- int.frac
            DoubleFloat(HI(x), floor(LO(x)))
        end
    else # HI(x) is mixed or fractional
        DoubleFloat(floor(HI(x)), zero(T))
    end
end

floor(::Type{Int128}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int128(floor(x))
floor(::Type{Int64}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int64(floor(x))
floor(::Type{Int32}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int32(floor(x))
floor(::Type{Int16}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int16(floor(x))
floor(::Type{Integer}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int(floor(x))

function ceil(x::DoubleFloat{T}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    # !isinteger(LO(x)), LO(x) is mixed or fractional
    if isinteger(HI(x))
        if isfractional(LO(x))
            if signbit(LO(x))
                DoubleFloat(HI(x), zero(T))
            else
                DoubleFloat(HI(x)+one(T), zero(T))
            end
        else # LO(x) is mixed: +/- int.frac
            DoubleFloat(HI(x), ceil(LO(x)))
        end
    else # HI(x) is mixed or fractional
        DoubleFloat(ceil(HI(x)), zero(T))
    end
end

ceil(::Type{Int128}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int128(ceil(x))
ceil(::Type{Int64}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int64(ceil(x))
ceil(::Type{Int32}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int32(ceil(x))
ceil(::Type{Int16}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int16(ceil(x))
ceil(::Type{Integer}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int(ceil(x))

function trunc(x::DoubleFloat{T}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    if isneg(x)
        ceil(x)
    else
        floor(x)
    end
end

trunc(::Type{Int128}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int128(trunc(x))
trunc(::Type{Int64}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int64(trunc(x))
trunc(::Type{Int32}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int32(trunc(x))
trunc(::Type{Int16}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int16(trunc(x))
trunc(::Type{Integer}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int(trunc(x))

function round(x::DoubleFloat{T}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    if isnonneg(x)
        trunc(x + 0.5)
    else
        -trunc(-x + 0.5)
    end
end

round(::Type{Int128}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int128(round(x))
round(::Type{Int64}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int64(round(x))
round(::Type{Int32}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int32(round(x))
round(::Type{Int16}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int16(round(x))
round(::Type{Integer}, x::DoubleFloat{T}) where {T<:AbstractFloat} = Int(round(x))

function round(x::DoubleFloat{T}, ::RoundingMode{:Up}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    return ceil(x)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:Down}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    return floor(x)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:RoundToZero}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    return isneg(x) ? ceil(x) : floor(x)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:RoundFromZero}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    return isneg(x) ? floor(x) : ceil(x)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:Nearest}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    isneg(x) && return -round(-x, RoundNearest)
    a = trunc(x + 0.5)
    return iseven(a) ? a : trunc(x - 0.5)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:NearestTiesAway}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    isneg(x) && return -round(-x, RoundNearestTiesAway)
    !isinteger(x - 0.5) && return round(x, RoundNearest)
    return round(x + 0.5, RoundNearest)
end

function round(x::DoubleFloat{T}, ::RoundingMode{:NearestTiesUp}) where {T<:AbstractFloat}
    (isinteger(x) || !isfinite(x)) && return x
    isneg(x) && return -round(-x, RoundNearestTiesUp)
    !isinteger(x - 0.5) && return round(x, RoundUp)
    return round(x + 0.5, RoundNearest)
end

"""
     spread(x)
the nearest integer to x, away from zero
spread complements trunc()
"""
function spread(x::DoubleFloat{T}) where {T<:AbstractFloat}
    (!isfinite(x) || isinteger(x)) && return x
    if isinteger(HI(x))
        if signbit(LO(x))
              signbit(HI(x)) ? DoubleFloat(HI(x)-one(T), zero(T)) : DoubleFloat(HI(x), zero(T))
        else
              signbit(HI(x)) ? DoubleFloat(HI(x), zero(T)) : DoubleFloat(HI(x)+one(T), zero(T))
        end
    else
        signbit(HI(x)) ? DoubleFloat(floor(HI(x)), zero(T)) :  DoubleFloat(ceil(HI(x)), zero(T))
    end
end

"""
    tld(x, y)
Truncate the result of x/y.
""" tld

"""
    sld(x, y)
Spread the result of x/y.
""" sld

for (F,G) in ((:fld, :floor), (:cld, :ceil), (:tld, :trunc), (:sld, :spread))
    @eval begin
        function $F(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
            z = x / y
            return $G(z)
        end
    end
end


text:
let ntraces = 10, max_N = 20, dat = fill(0.0, max_N, ntraces)
    for (idx, seed) in enumerate((abs(rand(Int)) for _ in 1:ntraces))
        cols = gen_bands(Cumulene(θ = π/2 + π/100), xsymbol = :N, bounds=(1, max_N), nsamples = 20, seed=seed)
        gen_near_bands!(cols, ysymbol = :E, bounds=(-10, 0), nsamples = 200, f = TheoryOfCISS.calc_data3, threshold = 1.0)
        xs, ys, zs = TheoryOfCISS.datamap(x -> x[:p], TheoryOfCISS.smoothen(cols, σ=10, w=TheoryOfCISS.box, p = x -> fullpolarization_reduced(x, 3)))
        dat[:, idx] = zs[100, :]
    end
    dat
end

#plot(100 .* zs[4, :], xlabel="N", ylabel="polarization (%)")

let ntraces = 10, Δ = π/8, max_N = 5, dat = fill(0.0, max_N, ntraces)
    for (idx, θ) in enumerate(range(π/2 - Δ, π/2 + Δ, length=ntraces))
        cols = gen_bands(Cumulene(θ = θ), xsymbol = :N, bounds=(1, max_N), nsamples = max_N, rng=1.0I)
        gen_near_bands!(cols, ysymbol = :E, bounds=(-10, 0), nsamples = 200, f = TheoryOfCISS.calc_data3, threshold = 1.0)
        xs, ys, zs = TheoryOfCISS.datamap(x -> x[:p], TheoryOfCISS.smoothen(cols, σ=10, w=TheoryOfCISS.box, p = x -> fullpolarization_reduced(x, 3)))
        dat[:, idx] = zs[100, :]
    end
    dat
end


text:
function ProjectVec(Fun,CG,Param)
OrdPoly=CG.OrdPoly;
OrdPolyZ=CG.OrdPolyZ;
nz=Param.Grid.nz;
uS=zeros(CG.NumG,nz,1);
vS=zeros(CG.NumG,nz,1);
uLoc=zeros(OrdPoly+1,OrdPoly+1);
vLoc=zeros(OrdPoly+1,OrdPoly+1);
JC = Param.cache.JC
X = Param.cache.X
for iz=1:nz
  for iF=1:Param.Grid.NumFaces
    for j=1:OrdPoly+1
      for i=1:OrdPoly+1
        x=0.5*(X[i,j,1,:,iF,iz]+X[i,j,2,:,iF,iz]);
        det=JC[i,j,iF,iz];
        (uLoc[i,j],vLoc[i,j])=Fun(x,Param);
        uLoc[i,j]=uLoc[i,j]*det;
        vLoc[i,j]=vLoc[i,j]*det;
      end
    end
    uS[CG.Faces[iF].Glob,iz]=uS[CG.Faces[iF].Glob,iz]+reshape(uLoc,(OrdPoly+1)*(OrdPoly+1),1);
    vS[CG.Faces[iF].Glob,iz]=vS[CG.Faces[iF].Glob,iz]+reshape(vLoc,(OrdPoly+1)*(OrdPoly+1),1);
  end
end
uS=uS./CG.M;
vS=vS./CG.M;
return (uS,vS)
end


text:
module ThreeOmegaMethod

using QuadGK

export ThreeOmega, LayerTOM, HeaterGeometry, Source

export three_omega

abstract type LayerTOMInformation end
struct LayerTOM{T1} <: LayerTOMInformation where {T1<:Float64}
    ky::T1
    kxy::T1
    d::T1
    ρC::T1
end

abstract type HeaterInformation end
struct HeaterGeometry{T1} <: HeaterInformation where {T1<:Float64}
    b::T1
    l::T1
    ρh::Array{T1}
end

abstract type HeaterSource end
struct Source{T1} <: HeaterSource where {T1<:Float64}
    p::T1
    f::Array{T1}
end

abstract type OutputT end
struct ThreeOmega{T1} <: OutputT where {T1<:Float64}
    f::Array{T1}
    ΔThreal::Array{T1}
    ΔThimag::Array{T1}
    integerror::Array{T1}
end

function three_omega(
    Layers::Array{T1,N1},
    heater::T2,
    source::T3,
    thresistances::Array{T4,N4};
    int_limit::T5=1.0e6,
) where {T1<:LayerTOM, N1, T2<:HeaterGeometry, T3<:Source, T4<:Number, N4, T5<:Number}
    numLayers::Int64 = size(Layers, 1)
    # Check input
    numLayers > 1 || throw(DimensionMismatch("There should be at least one layer below the heater."))
    length(thresistances) == numLayers-1 || throw(DimensionMismatch("The number of interface-resistances should be lower than the number of layers: length(thresistances) == size(layers,2)-1."))
    # Power per unit length
    p_l::Float64 = source.p / heater.l
    # Reverse order of input for bottom-up calculation
    Layers = reverse(Layers, dims=1)
    thresistances = reverse(vec(thresistances), dims=1)
    # Lambda parameters
    λ = [lambda(Layers[i].ky, Layers[i].kxy, Layers[i].ky / Layers[i].ρC) for i = 1 : numLayers]
    # Construct B parameter for the system of LayerTOMs
    B = _concatenate_B(λ, Layers, thresistances, numLayers)
    # Construct integration term. The half term comes from the heater at the top of the stack
    F::Function = (η,ω) -> (1.0./λ[end](η,ω).*(B(η,ω)[1] + B(η,ω)[2])./(0.5*B(η,ω)[2] - 0.5*B(η,ω)[1]).*(sin.(heater.b*η)./(heater.b*η)).^2)
    # return integration
    sol = _integrate_temperature(
            source.f,
            int_limit,
            F,
            heater.ρh,
            p_l/π/2,
            p_l/2/heater.b,
            im*4*π.*source.f,
    )
    return sol
end # ThreeOmegaMethod

"""Lambda parameter."""
lambda(ky::T1, kxy::T1, α::T1) where {T1<:Float64} = (η,ω) -> (ky*sqrt.(kxy*η^2 + im*2*ω/α))

"""Interface matrix."""
function Λ(λ1::T0, λ2::T0, r::T1) where {T0<:Function, T1<:Float64}
    x11 = (η,ω) -> λ2(η,ω) + λ1(η,ω) + λ2(η,ω).*λ1(η,ω).*r
    x12 = (η,ω) -> λ2(η,ω) - λ1(η,ω) - λ2(η,ω).*λ1(η,ω).*r
    x21 = (η,ω) -> λ2(η,ω) - λ1(η,ω) + λ2(η,ω).*λ1(η,ω).*r
    x22 = (η,ω) -> λ2(η,ω) + λ1(η,ω) - λ2(η,ω).*λ1(η,ω).*r
    X = (η,ω) -> (0.5./λ2(η,ω).*[x11(η,ω) x12(η,ω);
                                 x21(η,ω) x22(η,ω)])
    return X
end

"""Transfer matrix."""
function U(λ::T0, ky::T1, d::T1) where {T0<:Function, T1<:Float64}
    X = (η,ω) -> ([exp.(-λ(η,ω)/ky*d) 0.;
                   0. exp.(λ(η,ω)/ky*d)])
    return X
end

"""Integration of the temperature term."""
function _integrate_temperature(
    f::Array{T0,N0},
    int_limit::T0,
    F::T1,
    ρh::Array{T0,N2},
    plπ::T0,
    plb::T0,
    h_param::Array{T2,N0},
) where {T0<:Float64, N0, T1<:Function, T2<:ComplexF64, N2}
    ΔTh = Array{ComplexF64,1}(undef, length(f))
    int_error = Array{Float64,1}(undef, length(f))
    @inbounds for i in eachindex(f)
        temp0::Tuple{ComplexF64,Float64} = quadgk((η) -> F(η, 2*π*f[i]), 0, int_limit, rtol=sqrt.(eps()))
        temp1::ComplexF64 = temp0[1]
        int_error[i] = temp0[2]
        ΔTh[i] = ( plπ*temp1 + ρh[1]*plb ) ./ ( 1 + ρh[2]*h_param[i] * (ρh[1] + temp1/plb)  )
    end
    sol = ThreeOmega(f, real(ΔTh), imag(ΔTh), int_error)
    return sol
end # function _integrate_temperature()

"""

    Concatenate B terms depending on the number of LayerTOMs. The boundary condition for
    the substrate is semi-infinite.

"""
function _concatenate_B(
    λ::Array{T0,N0},
    Layers::Array{T1,N1},
    thresistances::Array{T2,N2},
    numLayers::T3,
) where {T0<:Function, N0, T1<:LayerTOM, N1, T2<:Number, N2, T3<:Int64}
    B = Λ(λ[end-1], λ[end], thresistances[end])
    @inbounds for i = (numLayers-1):-1:2
        tempU::Function = U(λ[i], Layers[i].ky, Layers[i].d)
        tempΛ::Function = Λ(λ[i-1], λ[i], thresistances[i])
        B = let B = B
            (η,ω) -> B(η,ω) * tempU(η,ω) * tempΛ(η,ω)
        end # B = let B = B
    end # i = numLayers-1 : -1 : 2
    # Add semi-infinite boundary condition
    B = let B = B
        (η,ω) -> B(η,ω) * [0.; 1.]
    end # B = let B = B
    return B
end # function _concatenate_B()

end # module


text:
module ProfileView

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Profile
using FlameGraphs
using FlameGraphs.IndirectArrays
using Base.StackTraces: StackFrame
using MethodAnalysis
using InteractiveUtils
using Gtk.ShortNames, GtkObservables, Colors, FileIO, IntervalSets
import Cairo
using Graphics
using Preferences

using FlameGraphs: Node, NodeData
using Gtk.GConstants.GdkModifierType: SHIFT, CONTROL, MOD1

export @profview, warntype_last

const clicked = Ref{Any}(nothing)   # for getting access to the clicked bar

const _graphtype = Ref{Symbol}(Symbol(@load_preference("graphtype", "flame")))

function set_graphtype!(graphtype::String)
    if !(graphtype in ("flame", "icicle"))
        throw(ArgumentError("Invalid graphtype: $graphtype. Valid options are :flame or :icicle"))
    end
    @set_preferences! "graphtype" => graphtype
    _graphtype[] = Symbol(graphtype)
    @info "Default graphtype set to $(repr(graphtype))"
    nothing
end
set_graphtype!(graphtype::Symbol) = set_graphtype!(string(graphtype))

"""
    warntype_last()
    warntype_last(io::IO; kwargs...)

Show `code_warntype` for the most recently-clicked bar.

Optionally direct output to stream `io`. Keyword arguments are passed to `code_warntype`.
"""
function warntype_last(io::IO=stdout; kwargs...)
    st = clicked[]
    if st === nothing || st.linfo === nothing
        @warn "click on a non-inlined bar to see `code_warntype` info"
        return nothing
    end
    return code_warntype(io, call_type(st.linfo.specTypes)...; kwargs...)
end

mutable struct ZoomCanvas
    bb::BoundingBox  # in user-coordinates
    c::Canvas
end

"""
    @profview f(args...)

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.
"""
macro profview(ex)
    return quote
        # pause the eventloop while profiling
        Gtk.pause_eventloop() do
            Profile.clear()
            @profile $(esc(ex))
        end
        view()
    end
end

"""
    ProfileView.closeall()

Closes all windows opened by ProfileView.
"""
function closeall()
    for (w, _) in window_wrefs
        destroy(w)
    end
    empty!(window_wrefs)   # see precompile.jl's usage of closeall
    return nothing
end

const window_wrefs = WeakKeyDict{Gtk.GtkWindowLeaf,Nothing}()
const tabname_allthreads = Symbol("All Threads")
const tabname_alltasks = Symbol("All Tasks")

NestedGraphDict = Dict{Symbol,Dict{Symbol,Node{NodeData}}}
"""
    ProfileView.view([fcolor], data=Profile.fetch(); lidict=nothing, C=false, recur=:off, fontsize=14, windowname="Profile", kwargs...)

View profiling results. `data` and `lidict` must be a matched pair from `Profile.retrieve()`.
You have several options to control the output, of which the major ones are:

- `fcolor`: an optional coloration function. The main options are `FlameGraphs.FlameColors`
  and `FlameGraphs.StackFrameCategory`.
- `C::Bool = false`: if true, the graph will include stackframes from C code called by Julia.
- `recur`: on Julia 1.4+, collapse recursive calls (see `Profile.print` for more detail)
- `expand_threads::Bool = true`: Break down profiling by thread (true by default)
- `expand_tasks::Bool = false`: Break down profiling of each thread by task (false by default)
- `graphtype::Symbol = :default`: Control how the graph is shown. `:flame` displays from the bottom up, `:icicle` from
  from the top down. The default type can be changed via e.g. `ProfileView.set_graphtype!(:icicle)`, which
  is stored as a preference for the active environment via `Preferences.jl`.

See [FlameGraphs](https://github.com/timholy/FlameGraphs.jl) for more information.
"""
function view(fcolor, data::Vector{UInt64}; lidict=nothing, C=false, combine=true, recur=:off, pruned=FlameGraphs.defaultpruned,
                expand_threads::Bool=true, expand_tasks::Bool=false, kwargs...)
    g = flamegraph(data; lidict=lidict, C=C, combine=combine, recur=recur, pruned=pruned)
    g === nothing && return nothing
    # Dict of dicts. Outer is threads, inner is tasks
    # Don't report the tasks at the "all threads" level because their id is thread-specific, so it's not useful
    # to track them across thread TODO: Perhaps fix that in base, so tasks keep the same id across threads?
    gdict = NestedGraphDict(tabname_allthreads => Dict{Symbol,Node{NodeData}}(tabname_alltasks => g))
    if expand_threads && isdefined(Profile, :has_meta) && Profile.has_meta(data)
        for threadid in Profile.get_thread_ids(data)
            g = flamegraph(data; lidict=lidict, C=C, combine=combine, recur=recur, pruned=pruned, threads = threadid)
            gdict_inner = Dict{Symbol,Node{NodeData}}(tabname_alltasks => g)
            if expand_tasks
                taskids = Profile.get_task_ids(data, threadid)
                if length(taskids) > 1
                    # skip when there's only one task as it will be the same as "all tasks"
                    for taskid in taskids
                        g = flamegraph(data; lidict=lidict, C=C, combine=combine, recur=recur, pruned=pruned, threads = threadid, tasks = taskid)
                        gdict_inner[Symbol(taskid)] = g
                    end
                end
            end
            gdict[Symbol(threadid)] = gdict_inner
        end
    end
    return view(fcolor, gdict; data=data, lidict=lidict, kwargs...)
end
function view(fcolor; kwargs...)
    data, lidict = Profile.retrieve()
    view(fcolor, data; lidict=lidict, kwargs...)
end
function view(data::Vector{UInt64}; lidict=nothing, kwargs...)
    view(FlameGraphs.default_colors, data; lidict=lidict, kwargs...)
end
function view(; kwargs...)
    # pausing the event loop here to facilitate a fast retrieve
    data, lidict = Gtk.pause_eventloop() do
        Profile.retrieve()
    end
    view(FlameGraphs.default_colors, data; lidict=lidict, kwargs...)
end

# This method allows user to open a *.jlprof file
viewblank() = (FlameGraphs.default_colors, Node(NodeData(StackTraces.UNKNOWN, 0, 1:0)))
view(::Nothing; kwargs...) = view(viewblank()...; kwargs...)

function view(g::Node{NodeData}; kwargs...)
    view(FlameGraphs.default_colors, g; kwargs...)
end
function view(fcolor, g::Node{NodeData}; data=nothing, lidict=nothing, kwargs...)
    win, _ = viewgui(fcolor, g; data=data, lidict=lidict, kwargs...)
    Gtk.showall(win)
end
function view(g_or_gdict::Union{Node{NodeData},NestedGraphDict}; kwargs...)
    view(FlameGraphs.default_colors, g_or_gdict; kwargs...)
end
function view(fcolor, g_or_gdict::Union{Node{NodeData},NestedGraphDict}; data=nothing, lidict=nothing, kwargs...)
    win, _ = viewgui(fcolor, g_or_gdict; data=data, lidict=lidict, kwargs...)
    Gtk.showall(win)
end

function viewgui(fcolor, g::Node{NodeData}; kwargs...)
    gdict = NestedGraphDict(tabname_allthreads => Dict{Symbol,Node{NodeData}}(tabname_alltasks => g))
    viewgui(fcolor, gdict; kwargs...)
end
function viewgui(fcolor, gdict::NestedGraphDict; data=nothing, lidict=nothing, windowname="Profile", graphtype = :default, kwargs...)
    _c, _fdraw, _tb_open, _tb_save_as = nothing, nothing, nothing, nothing # needed to be returned for precompile helper
    if graphtype == :default
        graphtype = _graphtype[]
    end
    thread_tabs = collect(keys(gdict))
    nb_threads = Notebook() # for holding the per-thread pages
    Gtk.GAccessor.scrollable(nb_threads, true)
    Gtk.GAccessor.show_tabs(nb_threads, length(thread_tabs) > 1)
    sort!(thread_tabs, by = s -> something(tryparse(Int, string(s)), 0)) # sorts thread_tabs as [all threads, 1, 2, 3 ....]

    for thread_tab in thread_tabs
        gdict_thread = gdict[thread_tab]
        task_tabs = collect(keys(gdict_thread))
        sort!(task_tabs, by = s -> s == tabname_alltasks ? "" : string(s)) # sorts thread_tabs as [all threads, 0xds ....]

        nb_tasks = Notebook() # for holding the per-task pages
        Gtk.GAccessor.scrollable(nb_tasks, true)
        Gtk.GAccessor.show_tabs(nb_tasks, length(task_tabs) > 1)
        task_tab_num = 1
        for task_tab in task_tabs
            g = gdict_thread[task_tab]
            gsig = Observable(g)  # allow substitution by the open dialog
            c = canvas(UserUnit)
            set_gtk_property!(widget(c), :expand, true)

            f = Frame(c)
            tb = Toolbar()
            tb_open = ToolButton("gtk-open")
            Gtk.GAccessor.tooltip_text(tb_open, "open")
            tb_save_as = ToolButton("gtk-save-as")
            Gtk.GAccessor.tooltip_text(tb_save_as, "save")
            tb_zoom_fit = ToolButton("gtk-zoom-fit")
            Gtk.GAccessor.tooltip_text(tb_zoom_fit, "zoom to fit")
            tb_zoom_in = ToolButton("gtk-zoom-in")
            Gtk.GAccessor.tooltip_text(tb_zoom_in, "zoom in")
            tb_zoom_out = ToolButton("gtk-zoom-out")
            Gtk.GAccessor.tooltip_text(tb_zoom_out, "zoom out")
            tb_info = ToolButton("gtk-info")
            Gtk.GAccessor.tooltip_text(tb_info, "ProfileView tips")

            push!(tb, tb_open)
            push!(tb, tb_save_as)
            push!(tb, SeparatorToolItem())
            push!(tb, tb_zoom_fit)
            push!(tb, tb_zoom_out)
            push!(tb, tb_zoom_in)
            push!(tb, SeparatorToolItem())
            push!(tb, tb_info)
            # FIXME: likely have to do `allkwargs` in the open/save below (add in C, combine, recur)
            signal_connect(open_cb, tb_open, "clicked", Nothing, (), false, (widget(c),gsig,kwargs))
            signal_connect(save_as_cb, tb_save_as, "clicked", Nothing, (), false, (widget(c),data,lidict,g))
            signal_connect(info_cb, tb_info, "clicked", Nothing, (), false, ())

            bx = Box(:v)
            push!(bx, tb)
            push!(bx, f)
            # don't use the actual taskid as the tab as it's very long
            push!(nb_tasks, bx, task_tab_num == 1 ? task_tab : Symbol(task_tab_num - 1))
            fdraw = viewprof(fcolor, c, gsig, (tb_zoom_fit, tb_zoom_out, tb_zoom_in), graphtype; kwargs...)
            GtkObservables.gc_preserve(nb_threads, c)
            GtkObservables.gc_preserve(nb_threads, fdraw)
            _c, _fdraw, _tb_open, _tb_save_as = c, fdraw, tb_open, tb_save_as
            task_tab_num += 1
        end
        push!(nb_threads, nb_tasks, string(thread_tab))
    end

    bx = Box(:v)
    push!(bx, nb_threads)

    # Defer creating the window until here because Window includes a `show` that will unpause the Gtk eventloop
    win = Window(windowname, 800, 600)
    push!(win, bx)

    # Register the window with closeall
    window_wrefs[win] = nothing
    signal_connect(win, :destroy) do w
        delete!(window_wrefs, win)
    end

    # Ctrl-w and Ctrl-q destroy the window
    signal_connect(win, "key-press-event") do w, evt
        if evt.state == CONTROL && (evt.keyval == UInt('q') || evt.keyval == UInt('w'))
            @async destroy(w)
            nothing
        end
    end

    return win, _c, _fdraw, (_tb_open, _tb_save_as)
end

function viewprof(fcolor, c, gsig, zoom_buttons, graphtype; fontsize=14)
    obs = on(gsig) do g
        viewprof_func(fcolor, c, g, fontsize, zoom_buttons, graphtype)
    end
    gsig[] = gsig[]
    return obs
end

function viewprof_func(fcolor, c, g, fontsize, zoom_buttons, graphtype)
    if !in(graphtype, (:flame, :icicle))
        throw(ArgumentError("Invalid option for `graphtype`: `$(repr(graphtype))`. Valid options are `:flame` and `:icicle`"))
    end
    tb_zoom_fit, tb_zoom_out, tb_zoom_in = zoom_buttons
    # From a given position, find the underlying tag
    function gettag(tagimg, xu, yu)
        x = ceil(Int, Float64(xu))
        y = ceil(Int, Float64(yu))
        Y = size(tagimg, 2)
        x = max(1, min(x, size(tagimg, 1)))
        y = max(1, min(y, Y))
        tagimg[x,Y-y+1]
    end
    function device_bb(c)
        if graphtype == :icicle
            BoundingBox(0, Graphics.width(c), Graphics.height(c), 0)
        elseif graphtype == :flame
            BoundingBox(0, Graphics.width(c), 0, Graphics.height(c))
        end
    end

    isempty(g.data.span) && return nothing
    img = flamepixels(fcolor, g)
    tagimg = flametags(g, img)
    # The first column corresponds to the bottom row, which is our fake root node. Get rid of it.
    img, tagimg = img[:,2:end], discardfirstcol(tagimg)
    img24 = RGB24.(img)
    img24 = img24[:,end:-1:1]
    fv = XY(0.0..size(img24,1), 0.0..size(img24,2))
    zr = Observable(ZoomRegion(fv, fv))
    signal_connect(zoom_fit_cb, tb_zoom_fit, "clicked", Nothing, (), false, (zr))
    signal_connect(zoom_out_cb, tb_zoom_out, "clicked", Nothing, (), false, (zr))
    signal_connect(zoom_in_cb, tb_zoom_in, "clicked", Nothing, (), false, (zr))
    sigrb = init_zoom_rubberband(c, zr)
    sigpd = init_pan_drag(c, zr)
    sigzs = init_zoom_scroll(c, zr)
    sigps = init_pan_scroll(c, zr)
    surf = Cairo.CairoImageSurface(img24)
    append!(c.preserved, Any[sigrb, sigpd, sigzs, sigps])
    let tagimg=tagimg    # julia#15276
        sigredraw = draw(c, zr) do widget, r
            ctx = getgc(widget)
            set_coordinates(ctx, device_bb(ctx), BoundingBox(r.currentview))
            rectangle(ctx, BoundingBox(r.currentview))
            set_source(ctx, surf)
            p = Cairo.get_source(ctx)
            Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
            fill(ctx)
        end
        lasttextbb = Ref(BoundingBox(1,0,1,0))
        sigmotion = on(c.mouse.motion) do btn
            # Repair image from ovewritten text
            if c.widget.is_realized && c.widget.is_sized
                ctx = getgc(c)
                if Graphics.width(lasttextbb[]) > 0
                    r = zr[]
                    set_coordinates(ctx, device_bb(ctx), BoundingBox(r.currentview))
                    rectangle(ctx, lasttextbb[])
                    set_source(ctx, surf)
                    p = Cairo.get_source(ctx)
                    Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
                    fill(ctx)
                end
                # Write the info
                xu, yu = btn.position.x, btn.position.y
                sf = gettag(tagimg, xu, yu)
                if sf != StackTraces.UNKNOWN
                    str = string(basename(string(sf.file)), ", ", sf.func, ": line ", sf.line)
                    set_source(ctx, fcolor(:font))
                    Cairo.set_font_face(ctx, "sans-serif $(fontsize)px")
                    xi = zr[].currentview.x
                    xmin, xmax = minimum(xi), maximum(xi)
                    lasttextbb[] = deform(Cairo.text(ctx, xu, yu, str, halign = xu < (2xmin+xmax)/3 ? "left" : xu < (xmin+2xmax)/3 ? "center" : "right"), -2, 2, -2, 2)
                end
                reveal(c)
            end
        end
        # Left-click prints the full path, function, and line to the console
        # Right-click calls the edit() function
        sigshow = on(c.mouse.buttonpress) do btn
            if btn.button == 1 || btn.button == 3
                ctx = getgc(c)
                xu, yu = btn.position.x, btn.position.y
                sf = gettag(tagimg, xu, yu)
                clicked[] = sf
                if sf != StackTraces.UNKNOWN
                    if btn.button == 1
                        if sf.linfo isa Core.MethodInstance
                            println(sf.file, ':', sf.line, ", ", sf.linfo)
                        else
                            println(sf.file, ':', sf.line, ", ", sf.func, " [inlined]")
                        end
                    elseif btn.button == 3
                        edit(string(sf.file), sf.line)
                    end
                end
            end
        end
        append!(c.preserved, Any[sigredraw, sigmotion, sigshow])
    end
    return nothing
end

@guarded function open_cb(::Ptr, settings::Tuple)
    c, gsig, kwargs = settings
    selection = open_dialog("Load profile data", toplevel(c), ("*.jlprof","*"))
    isempty(selection) && return nothing
    return _open(gsig, selection; kwargs...)
end

function _open(gsig, selection; kwargs...)
    ret = load(selection)
    if isa(ret, Node{NodeData})
        gsig[] = ret
    else
        data, lidict = ret::Tuple{Vector{UInt64},Profile.LineInfoDict}
        gsig[] = flamegraph(data; lidict=lidict, kwargs...)
    end
    return nothing
end

@guarded function save_as_cb(::Ptr, profdata::Tuple)
    c, data, lidict, g = profdata
    selection = save_dialog("Save profile data as *.jlprof file", toplevel(c), ("*.jlprof",))
    isempty(selection) && return nothing
    if data === nothing && lidict === nothing
        return _save(selection, g)
    end
    return _save(selection, data, lidict)
end

function _save(selection, args...)
    FileIO.save(File{format"JLPROF"}(selection), args...)
    return nothing
end

@guarded function zoom_fit_cb(::Ptr, zr::Observable{ZoomRegion{T}}) where {T}
    zr[] = GtkObservables.reset(zr[])
    return nothing
end

@guarded function zoom_in_cb(::Ptr, zr::Observable{ZoomRegion{T}}) where {T}
    setindex!(zr, zoom(zr[], 1/2))
    return nothing
end

@guarded function zoom_out_cb(::Ptr, zr::Observable{ZoomRegion{T}}) where {T}
    setindex!(zr, zoom(zr[], 2))
    return nothing
end

@guarded function info_cb(::Ptr, ::Tuple)
    # Note: Keep this updated with the readme
    info = """
    ProfileView.jl Interface Tips
    ----------------------------------------------

    `Ctrl-q` and `Ctrl-w` close the window. You can also use `ProfileView.closeall()` to close all windows opened by ProfileView.

    Left-clicking on a bar will cause information about this line to be printed in the REPL. This can be a convenient way to "mark" lines for later investigation.

    Right-clicking on a bar calls the `edit()` function to open the line in an editor. (On a trackpad, use a 2-fingered tap.)

    CTRL-clicking and dragging will zoom in on a specific region of the image. You can also control the zoom level with CTRL-scroll (or CTRL-swipe up/down).

    CTRL-double-click to restore the full view.

    You can pan the view by clicking and dragging, or by scrolling your mouse/trackpad (scroll=vertical, SHIFT-scroll=horizontal).

    The toolbar at the top includes icons to load and save profile data. Clicking the save icon will prompt you for a filename; you should use extension *.jlprof for any file you save. Launching `ProfileView.view(nothing)` opens a blank window, which you can populate with saved data by clicking on the "open" icon.

    After clicking on a bar, you can type `warntype_last()` and see the result of `code_warntype` for the call represented by that bar.

    `ProfileView.view(windowname="method1")` allows you to name your window, which can help avoid confusion when opening several ProfileView windows simultaneously.

    On Julia 1.8 `ProfileView.view(expand_tasks=true)` creates one tab per task. Expanding by thread is on by default and can be disabled with `expand_threads=false`.

    Using the `graphtype` kwarg for `ProfileView.view` controls how the graph is shown. `:flame` displays from the bottom up, `:icicle` from the top down. The default type can be changed via e.g. `ProfileView.set_graphtype!(:icicle)`, which is stored as a preference for the active environment via `Preferences.jl`.
    """
    info_dialog(info)
    return nothing
end

discardfirstcol(A) = A[:,2:end]
discardfirstcol(A::IndirectArray) = IndirectArray(A.index[:,2:end], A.values)

if ccall(:jl_generating_output, Cint, ()) == 1
    include("precompile.jl")
    _precompile_()
end

end


text:
# ---
# title: 1620. Coordinate With Maximum Network Quality
# id: problem1620
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/coordinate-with-maximum-network-quality/description/>
# hidden: true
# ---
#
# You are given an array of network towers `towers` and an integer `radius`,
# where `towers[i] = [xi, yi, qi]` denotes the `ith` network tower with location
# `(xi, yi)` and quality factor `qi`. All the coordinates are **integral
# coordinates** on the X-Y plane, and the distance between two coordinates is
# the **Euclidean distance**.
#
# The integer `radius` denotes the **maximum distance** in which the tower is
# **reachable**. The tower is **reachable** if the distance is less than or
# equal to `radius`. Outside that distance, the signal becomes garbled, and the
# tower is **not reachable**.
#
# The signal quality of the `ith` tower at a coordinate `(x, y)` is calculated
# with the formula `⌊qi / (1 + d)⌋`, where `d` is the distance between the tower
# and the coordinate. The **network quality** at a coordinate is the sum of the
# signal qualities from all the **reachable** towers.
#
# Return _the integral coordinate where the **network quality** is maximum_. If
# there are multiple coordinates with the same **network quality** , return _the
# lexicographically minimum coordinate_.
#
# **Note:**
#
#   * A coordinate `(x1, y1)` is lexicographically smaller than `(x2, y2)` if either `x1 < x2` or `x1 == x2` and `y1 < y2`.
#   * `⌊val⌋` is the greatest integer less than or equal to `val` (the floor function).
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/22/untitled-diagram.png)
#
#
#
#     Input: towers = [[1,2,5],[2,1,7],[3,1,9]], radius = 2
#     Output: [2,1]
#     Explanation:
#     At coordinate (2, 1) the total quality is 13
#     - Quality of 7 from (2, 1) results in ⌊7 / (1 + sqrt(0)⌋ = ⌊7⌋ = 7
#     - Quality of 5 from (1, 2) results in ⌊5 / (1 + sqrt(2)⌋ = ⌊2.07⌋ = 2
#     - Quality of 9 from (3, 1) results in ⌊9 / (1 + sqrt(1)⌋ = ⌊4.5⌋ = 4
#     No other coordinate has higher quality.
#
# **Example 2:**
#
#
#
#     Input: towers = [[23,11,21]], radius = 9
#     Output: [23,11]
#
#
# **Example 3:**
#
#
#
#     Input: towers = [[1,2,13],[2,1,7],[0,1,9]], radius = 2
#     Output: [1,2]
#
#
# **Example 4:**
#
#
#
#     Input: towers = [[2,1,9],[0,1,9]], radius = 2
#     Output: [0,1]
#     Explanation: Both (0, 1) and (2, 1) are optimal in terms of quality but (0, 1) is lexicograpically minimal.
#
#
#
#
# **Constraints:**
#
#   * `1 <= towers.length <= 50`
#   * `towers[i].length == 3`
#   * `0 <= xi, yi, qi <= 50`
#   * `1 <= radius <= 50`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export projGN, projGNCG

"""
	function projGrad

	Projects gradient, i.e.,

					 | gc[i],  				xl[i] < x[i] < xh[i]
	gc[i]  = | max(gc[i],0)	 	xc[i] == xh[i]
					 | min(gc[i],0)   xc[i] == xl[i]

	Input:

		gc 				  - gradient vector
		mc 				  - model
		boundsLow   - lower bounds
		boundsHigh  - upper bounds
"""
function projGrad(gc,mc,boundsLow,boundsHigh)
	pgc = gc[:]
	pgc[mc.==boundsLow] = min.(gc[mc.==boundsLow],0)
	pgc[mc.==boundsHigh] = max.(gc[mc.==boundsHigh],0)
	return pgc
end

function dummy(mc,Dc,iter,pInv,pMis)
# this function does nothing and is used as a default for dumpResults().
end;

"""
	mc,Dc,outerFlag = projGN(mc,pInv::InverseParam,pMis, indFor = [], dumpResults::Function = dummy)

	(Projected) Gauss-Newton method.

	Input:

		mc::Vector          - intial guess for model
		pInv::InverseParam  - parameter for inversion
		pMis                - misfit terms
		indCredit           - indices of forward problems to work on
		dumpResults			- A function pointer for saving the results throughout the iterations.
							- We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
							- where mc is the recovered model, Dc is the predicted data.
							- If dumpResults is not given, nothing is done (dummy() is called).
		out::Int            - flag for output (-1: no output, 1: final status, 2: residual norm at each iteration)
		solveGN::Function   - solver for GN system, default: projPCG, other options: projGNexplicit
		                      The interface for this method is:
							            dm,his = solveGN(gc,pMis,pInv,sig,dsig,d2F,d2R,Active)

	Output:
		mc                  - final model
		Dc                  - data
		outerFlag           - flag for convergence
		His                 - iteration history

"""
function  projGN(mc,pInv::InverseParam,pMis;indCredit=[],
	dumpResults::Function = dummy,out::Int=2,solveGN::Function=projPCG)

	maxIter     = pInv.maxIter      #  Max. no. iterations.
	pcgMaxIter  = pInv.pcgMaxIter   #  Max cg iters.
	pcgTol      = pInv.pcgTol       #  CG stopping tolerance.
	stepTol     = pInv.minUpdate    #  Step norm stopping tol.
	maxStep     = pInv.maxStep
	low         = pInv.boundsLow
	high        = pInv.boundsHigh
	alpha       = pInv.alpha

	His = getGNhis(maxIter,pcgMaxIter)
	#---------------------------------------------------------------------------
	#  Initialization.
	#---------------------------------------------------------------------------

	Active = (mc .<=low) .| (mc.>=high)  # Compute active set

	## evaluate function and derivatives
	sig,dsig = pInv.modelfun(mc)
	if isempty(indCredit)
		Dc,F,dF,d2F,pMis,tMis = computeMisfit(sig,pMis,true)
	else
		Dc,F,dF,d2F,pMis,tMis,indDebit = computeMisfit(sig,pMis,true,indCredit=indCredit)
	end
	dF = dsig'*dF
	dumpResults(mc,Dc,0,pInv,pMis) # dump initial model and dpred0

	# compute regularizer

	tReg = @elapsed R,dR,d2R = computeRegularizer(pInv.regularizer,mc,pInv.mref,pInv.MInv,alpha)

	# objective function
	Jc  = F  + R
	gc  = dF + dR

	F0 = F; J0 = Jc
	############################################################################
	##  Outer iteration.                                                        #
	############################################################################
	iter = 0
	outerFlag = -1; stepNorm=0.0

	outStr = @sprintf("\n %4s\t%08s\t%08s\t%08s\t%08s\t%08s\n",
					  	"i.LS", "F", "R","alpha[1]","Jc/J0","#Active")
	updateHis!(0,His,Jc,norm(projGrad(gc,mc,low,high)),F,Dc,R,alpha[1],count(!iszero, Active),0.0,-1,tMis,tReg)

	if out>=2; print(outStr); end
	f = open("jInv.out", "w")
	write(f, outStr)
	close(f)

	while outerFlag == -1

		iter += 1
		outStr = @sprintf("%3d.0\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3d\n",
		         iter, F, R,alpha[1],Jc/J0,count(!iszero, Active))
		if out>=2; print(outStr); end
		f = open("jInv.out", "a")
		write(f, outStr)
		close(f)

		# solve linear system to find search direction
		His.timeLinSol[iter+1] += @elapsed  delm,hisLinSol = solveGN(gc,pMis,pInv,sig,dsig,d2F,d2R,Active)
		push!(His.hisLinSol,hisLinSol)

		# scale step
		if maximum(abs.(delm)) > maxStep; delm = delm./maximum(abs.(delm))*maxStep; end

		# take gradient direction in the active cells
		ga = -gc[mc .== low]
		if !isempty(ga)
			if maximum(abs.(ga)) > maximum(abs.(delm)); ga = ga./maximum(abs.(ga))*maximum(abs.(delm)); end
			delm[mc .== low] = ga
		end
		ga = -gc[mc .== high]
		if !isempty(ga)
			if maximum(abs.(ga)) > maximum(abs.(delm)); ga = ga./maximum(abs.(ga))*maximum(abs.(delm)); end
			delm[mc .== high] = ga
		end

		## Begin projected Armijo line search
		muLS = 1; lsIter = 1; mt = zeros(size(mc)); Jt = Jc
		while true
			mt = mc + muLS*delm
			mt[mt.<low]  = low[mt.<low]
			mt[mt.>high] = high[mt.>high]

			His.timeReg[iter+1] += @elapsed R,dR,d2R = computeRegularizer(pInv.regularizer,mt,pInv.mref,pInv.MInv,alpha)

			if R!=Inf # to support barrier functions.
				## evaluate function

				sigt, = pInv.modelfun(mt)
				if isempty(indCredit)
					Dc,F,dF,d2F,pMis,tMis = computeMisfit(sigt,pMis,false)
				else
					Dc,F,dF,d2F,pMis,tMis,indDebit = computeMisfit(sigt,false,indCredit=indCredit)
				end
				His.timeMisfit[iter+1,:]+=tMis

				# objective function
				Jt  = F  + R
				if out>=2;
					println(@sprintf( "   .%d\t%3.2e\t%3.2e\t\t\t%3.2e",
						lsIter, F,       R,       Jt/J0))
				end

				if Jt < Jc
					break
				end
			else
				Jt  = R;
				F = 0.0;
				if out>=2;
					println(@sprintf( "   .%d\t%3.2e\t%3.2e\t\t\t%3.2e",
						lsIter, F,       R,       Jt/J0))
				end
			end
			muLS /=2; lsIter += 1
			if lsIter > 6
			    outerFlag = -2
				break
			end
		end
		## End Line search

		## Check for termination
		stepNorm = norm(mt-mc,Inf)
		mc = mt
		Jc = Jt

		sig, dsig = pInv.modelfun(mc)

		Active = (mc .<=low) .| (mc.>=high)  # Compute active set

		#  Check stopping criteria for outer iteration.
		updateHis!(iter,His,Jc,-1.0,F,Dc,R,alpha[1],count(!iszero, Active),stepNorm,lsIter,tMis,tReg)

		dumpResults(mc,Dc,iter,pInv,pMis);
		if stepNorm < stepTol
			outerFlag = 1
			break
		elseif iter >= maxIter
			break
		end
		# Evaluate gradient

		if isempty(indCredit)
			His.timeGradMisfit[iter+1]+= @elapsed dF = computeGradMisfit(sig,Dc,pMis)
		else
			His.timeGradMisfit[iter+1]+= @elapsed dF = computeGradMisfit(sig,Dcp,pMis,indDebit)
		end

		dF = dsig'*dF
		gc = dF + dR

		His.dJ[iter+1] = norm(projGrad(gc,mc,low,high))

	end # while outer_flag == 0

	if out>=1
		if outerFlag==-1
			println("projGN iterated maxIter=$maxIter times but reached only stepNorm of $(stepNorm) instead $(stepTol)." )
		elseif outerFlag==-2
			println("projGN stopped at iteration $iter with a line search fail.")
		elseif outerFlag==1
			println("projGN reached desired accuracy at iteration $iter.")
		end
	end

	return mc,Dc,outerFlag,His
end  # Optimization code

projGNCG = projGN


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ThrottledRequestsInput <: SwaggerModel
    blobContainerSasUri::Any # spec type: Union{ Nothing, String } # spec name: blobContainerSasUri
    fromTime::Any # spec type: Union{ Nothing, DateTime } # spec name: fromTime
    toTime::Any # spec type: Union{ Nothing, DateTime } # spec name: toTime
    groupByThrottlePolicy::Any # spec type: Union{ Nothing, Bool } # spec name: groupByThrottlePolicy
    groupByOperationName::Any # spec type: Union{ Nothing, Bool } # spec name: groupByOperationName
    groupByResourceName::Any # spec type: Union{ Nothing, Bool } # spec name: groupByResourceName

    function ThrottledRequestsInput(;blobContainerSasUri=nothing, fromTime=nothing, toTime=nothing, groupByThrottlePolicy=nothing, groupByOperationName=nothing, groupByResourceName=nothing)
        o = new()
        validate_property(ThrottledRequestsInput, Symbol("blobContainerSasUri"), blobContainerSasUri)
        setfield!(o, Symbol("blobContainerSasUri"), blobContainerSasUri)
        validate_property(ThrottledRequestsInput, Symbol("fromTime"), fromTime)
        setfield!(o, Symbol("fromTime"), fromTime)
        validate_property(ThrottledRequestsInput, Symbol("toTime"), toTime)
        setfield!(o, Symbol("toTime"), toTime)
        validate_property(ThrottledRequestsInput, Symbol("groupByThrottlePolicy"), groupByThrottlePolicy)
        setfield!(o, Symbol("groupByThrottlePolicy"), groupByThrottlePolicy)
        validate_property(ThrottledRequestsInput, Symbol("groupByOperationName"), groupByOperationName)
        setfield!(o, Symbol("groupByOperationName"), groupByOperationName)
        validate_property(ThrottledRequestsInput, Symbol("groupByResourceName"), groupByResourceName)
        setfield!(o, Symbol("groupByResourceName"), groupByResourceName)
        o
    end
end # type ThrottledRequestsInput

const _property_map_ThrottledRequestsInput = Dict{Symbol,Symbol}(Symbol("blobContainerSasUri")=>Symbol("blobContainerSasUri"), Symbol("fromTime")=>Symbol("fromTime"), Symbol("toTime")=>Symbol("toTime"), Symbol("groupByThrottlePolicy")=>Symbol("groupByThrottlePolicy"), Symbol("groupByOperationName")=>Symbol("groupByOperationName"), Symbol("groupByResourceName")=>Symbol("groupByResourceName"))
const _property_types_ThrottledRequestsInput = Dict{Symbol,String}(Symbol("blobContainerSasUri")=>"String", Symbol("fromTime")=>"DateTime", Symbol("toTime")=>"DateTime", Symbol("groupByThrottlePolicy")=>"Bool", Symbol("groupByOperationName")=>"Bool", Symbol("groupByResourceName")=>"Bool")
Base.propertynames(::Type{ ThrottledRequestsInput }) = collect(keys(_property_map_ThrottledRequestsInput))
Swagger.property_type(::Type{ ThrottledRequestsInput }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ThrottledRequestsInput[name]))}
Swagger.field_name(::Type{ ThrottledRequestsInput }, property_name::Symbol) =  _property_map_ThrottledRequestsInput[property_name]

function check_required(o::ThrottledRequestsInput)
    (getproperty(o, Symbol("blobContainerSasUri")) === nothing) && (return false)
    (getproperty(o, Symbol("fromTime")) === nothing) && (return false)
    (getproperty(o, Symbol("toTime")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ThrottledRequestsInput }, name::Symbol, val)
end


text:
export NFFTAlg

using NFFT

function NFFTAlg(obs::EHTObservation; kwargs...)
    u, v = getuv(arrayconfig(obs))
    return NFFTAlg(u, v; kwargs...)
end

function NFFTAlg(u::AbstractArray, v::AbstractArray; padfac=1, m=10)
    uv = Matrix{eltype(u)}(undef, 2, length(u))
    uv[1,:] .= u
    uv[2,:] .= v
    return ObservedNUFT(NFFTAlg(;padfac, m), uv)
end

function NFFTAlg(ac::ArrayConfiguration; padfac=1, m=10)
    u, v = getuv(ac)
    return NFFTAlg(u, v; padfac, m)
end

function padimage(alg::NFFTAlg, img)
    padfac = alg.padfac
    # if no padding exit now
    (padfac == 1) && return img

    ny,nx = size(img)
    nnx = nextpow(2, padfac*nx)
    nny = nextpow(2, padfac*ny)
    nsx = nnx÷2-nx÷2
    nsy = nny÷2-ny÷2
    return PaddedView(zero(eltype(img)), img,
                      (1:nnx, 1:nny),
                      (nsx+1:nsx+nx, nsy+1:nsy+ny)
                     )
end

function plan_nuft(alg::ObservedNUFT{<:NFFTAlg}, img, dx, dy)
    uv2 = similar(alg.uv)
    uv2[1,:] .= alg.uv[1,:]*dx
    uv2[2,:] .= alg.uv[2,:]*dy
    plan = plan_nfft(uv2, size(img); precompute=NFFT.POLYNOMIAL)
    return plan
end

function make_phases(alg::ObservedNUFT{<:NFFTAlg}, img)
    dx, dy = pixelsizes(img)
    u = @view alg.uv[1,:]
    v = @view alg.uv[2,:]
    return cispi.((u.*dx .+ v.*dy)).*visibility_point.(Ref(img.pulse), u.*dx, v.*dy)
end

@inline function create_cache(alg::ObservedNUFT{<:NFFTAlg}, plan, phases, img)
    timg = IntensityMap(transpose(img.im), img.fovx, img.fovy, img.pulse)
    return NUFTCache(alg, plan, phases, img.pulse, timg)
end

function _frule_vis(m::ModelImage{M,<:IntensityMap{<:ForwardDiff.Dual{T,V,P}},<:NUFTCache{O}}) where {M,T,V,P,A<:NFFTAlg,O<:ObservedNUFT{A}}
    S = typeof(ForwardDiff.value(first(m.cache.img)))
    p = m.cache.plan
    # Compute the fft
    buffer = similar(m.cache.img, S)
    buffer .= ForwardDiff.value.(m.cache.img)
    xtil = p*(buffer .+ 0.0im)
    out = similar(buffer, Complex{ForwardDiff.Dual{T,V,P}})
    # Now take the deriv of nuft
    ndxs = ForwardDiff.npartials(first(m.cache.img))
    dxtils = ntuple(ndxs) do n
        buffer .= ForwardDiff.partials.(m.cache.img, n)
        p * (buffer .+ 0.0im)
    end
    out = similar(xtil, Complex{ForwardDiff.Dual{T,V,P}})
    for i in eachindex(out)
        dual = getindex.(dxtils, i)
        prim = xtil[i]
        red = ForwardDiff.Dual{T,V,P}(real(prim), ForwardDiff.Partials(real.(dual)))
        imd = ForwardDiff.Dual{T,V,P}(imag(prim), ForwardDiff.Partials(imag.(dual)))
        out[i] = Complex(red, imd)
    end
    return out
end

function _visibilities(m::ModelImage{M,<:IntensityMap{<:ForwardDiff.Dual{T,V,P}},<:NUFTCache{O}},
    u::AbstractArray,
    v::AbstractArray) where {M,T,V,P,A<:NFFTAlg,O<:ObservedNUFT{A}}
    checkuv(m.cache.alg.uv, u, v)
    # Now reconstruct everything

    vis = _frule_vis(m)
    conj.(vis).*m.cache.phases
#return vis
end


text:
module BIFMNodeTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "BIFMNode" begin

    @testset "Creation" begin
        node = make_node(BIFM, meta = BIFMMeta(ones(2,2), 2*ones(2,2), 3*ones(2,2)))

        @test functionalform(node)          === BIFM
        @test sdtype(node)                  === Deterministic()
        @test name.(interfaces(node))       === (:out, :in, :zprev, :znext)
        @test factorisation(node)           === ((1, 2, 3, 4), )
        @test metadata(node).A              ==  ones(2,2)
        @test metadata(node).B              ==  2*ones(2,2)
        @test metadata(node).C              ==  3*ones(2,2)

        node = make_node(BIFM, meta = 1)

        @test functionalform(node)          === BIFM
        @test sdtype(node)                  === Deterministic()
        @test name.(interfaces(node))       === (:out, :in, :zprev, :znext)
        @test factorisation(node)           === ((1, 2, 3, 4), )
        @test metadata(node)                === 1
    end
end
end


text:
# easier creation of atoms

function (p::Predicate)(args::Vararg{Union{Term,Int64,Float64}})
    Literal(p, [x for x in args])
end

function (f::Functor)(args::Vararg{Union{Term,Int64,Float64}})
    Structure(f, [x for x in args])
end

# \+(literal::Literal) = Negation(literal)

Base.:!(literal::Literal) = Negation(literal)

"""
    Easier syntax for constructing clauses
"""
function Base.:&(left::Union{Literal,Negation,Proposition}, right::Union{Negation,Literal,Proposition})
    Conj([left, right])
end

function Base.:&(left::Conj, right::Union{Literal,Negation,Proposition})
    Conj(vcat(left, right))
end

function Base.:(<=)(left::Union{Literal,Proposition}, right::Conj)
    Clause(left, right)
end

function Base.:(<=)(left::Union{Literal,Proposition}, right::Union{Literal,Negation})
    Clause(left, Conj([right]))
end


text:
fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
anyerrors = false

my_tests = [
    "utils/maxdiff.jl",
    "automatic/dualnumbers.jl",
    "symbolic/symbolic.jl",
    "symbolic/deparse.jl",
    "api/grad.jl",
    "api/hess.jl",
    "api/checker.jl",
    "api/integrate.jl",
]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
using Gtk
using Test
using DoseCalculators

function amoxicillin_po(; age, weight, height)
    if age <= 3/12
        return weight * 30.0
    elseif age <= 16 && weight <= 40
        return weight * 37.0
    else
        return 1500.0
    end
end

DoseCalculators._apps_should_persist[1] = false

dose_calculator_app(amoxicillin_po, "Amoxacillin PO Dosing", "Amoxacillin dosage(mg)")

@test DoseCalculators._apps[end] isa GtkWindow


text:
using ExplainableAI: augment_batch_dim, augment_indices, reduce_augmentation
using ExplainableAI: interpolate_batch

# augment_batch_dim
A = [1 2; 3 4]
B = @inferred augment_batch_dim(A, 3)
@test B == [
    1 1 1 2 2 2
    3 3 3 4 4 4
]
B = @inferred augment_batch_dim(A, 4)
@test B == [
    1 1 1 1 2 2 2 2
    3 3 3 3 4 4 4 4
]
A = reshape(1:8, 2, 2, 2)
B = @inferred augment_batch_dim(A, 3)
@test B[:, :, 1] == A[:, :, 1]
@test B[:, :, 2] == A[:, :, 1]
@test B[:, :, 3] == A[:, :, 1]
@test B[:, :, 4] == A[:, :, 2]
@test B[:, :, 5] == A[:, :, 2]
@test B[:, :, 6] == A[:, :, 2]

# augment_batch_dim
inds = [CartesianIndex(5, 1), CartesianIndex(3, 2)]
augmented_inds = @inferred augment_indices(inds, 3)
@test augmented_inds == [
    CartesianIndex(5, 1)
    CartesianIndex(5, 2)
    CartesianIndex(5, 3)
    CartesianIndex(3, 4)
    CartesianIndex(3, 5)
    CartesianIndex(3, 6)
]

# reduce_augmentation
A = Float32.(reshape(1:10, 1, 1, 10))
R = @inferred reduce_augmentation(A, 5)
@test R == reshape([sum(1:5), sum(6:10)] / 5, 1, 1, :)
A = Float64.(reshape(1:10, 1, 1, 1, 1, 10))
R = @inferred reduce_augmentation(A, 2)
@test R == reshape([3, 7, 11, 15, 19] / 2, 1, 1, 1, 1, :)

x = Float16.(reshape(1:4, 2, 2))
x0 = zero(x)
A = @inferred interpolate_batch(x, x0, 5)
@test A ≈ [
    0.0 0.25 0.5 0.75 1.0 0.0 0.75 1.5 2.25 3.0
    0.0 0.5 1.0 1.5 2.0 0.0 1.0 2.0 3.0 4.0
]


text:
p = pwd()
cd(Pkg.dir("RNGTest/deps/"))
run(`make`)
cd(p)

max_stars_repo_name: josejuanmijares/RNGTest

max_issues_repo_name: josejuanmijares/RNGTest

max_forks_repo_name: josejuanmijares/RNGTest


text:
# function require_julia_version(varargs...)::VersionNumber
#     current_julia_version = convert(VersionNumber, Base.VERSION)
#     version_meets_requirements = does_given_version_meet_requirements(
#         current_julia_version,
#         collect(varargs),
#         )
#     if !version_meets_requirements
#         error(
#             string(
#                 "Current Julia version (",
#                 current_julia_version,
#                 ") does not match any of the ",
#                 "user-specified version intervals.",
#                 )
#             )
#     end
#     return current_julia_version
# end

# function require_predictmd_version(varargs...)::VersionNumber
#     current_predictmd_version = convert(VersionNumber, version())
#     version_meets_requirements = does_given_version_meet_requirements(
#         current_predictmd_version,
#         collect(varargs),
#         )
#     if !version_meets_requirements
#         error(
#             string(
#                 "Current PredictMD version (",
#                 current_predictmd_version,
#                 ") does not match any of the ",
#                 "user-specified version intervals.",
#                 )
#             )
#     end
#     return current_predictmd_version
# end

# function does_given_version_meet_requirements(
#         given_version::VersionNumber,
#         version_requirements::AbstractVector,
#         )::Bool
#     num_version_requirements = length(version_requirements)
#     if num_version_requirements == 0
#         answer = true
#     else
#         if iseven(num_version_requirements)
#             num_intervals = Int((num_version_requirements)/(2))
#             answer_for_each_interval = Vector{Bool}(undef, num_intervals)
#             for interval = (1):(num_intervals)
#                 answer_for_each_interval[interval] =
#                     does_given_version_meet_requirements(
#                         given_version,
#                         version_requirements[2*interval - 1],
#                         version_requirements[2*interval],
#                         )
#             end
#         else
#             num_intervals = Int((num_version_requirements+1)/(2))
#             answer_for_each_interval = Vector{Bool}(undef, num_intervals)
#             for interval = (1):(num_intervals - 1)
#                 answer_for_each_interval[interval] =
#                     does_given_version_meet_requirements(
#                         given_version,
#                         version_requirements[2*interval - 1],
#                         version_requirements[2*interval],
#                         )
#             end
#             answer_for_each_interval[num_intervals] =
#                 does_given_version_meet_requirements(
#                     given_version,
#                     version_requirements[2*num_intervals - 1],
#                     )
#         end
#         answer = any(answer_for_each_interval)
#     end
#     return answer
# end

# function does_given_version_meet_requirements(
#         given_version::VersionNumber,
#         min_version,
#         )::Bool
#     min_version = VersionNumber(min_version)
#     answer = min_version <= given_version
#     return answer
# end

# function does_given_version_meet_requirements(
#         given_version::VersionNumber,
#         min_version,
#         max_version,
#         )::Bool
#     min_version = VersionNumber(min_version)
#     max_version = VersionNumber(max_version)
#     answer = min_version <= given_version < max_version
#     return answer
# end


text:
md"""
# Candidate Utilities

As seen in [Fitness Functions](@ref), fitness strategies require an `AbstractCandidate` to compute fitness. To be used by NaiveGAflux,
an `AbstractCandidate` needs to
1. Provide the data needed by the fitness strategy, most commonly the model but also things like lossfunctions and optimizers
2. Be able to create a new version of itself given a function which maps its fields to new fields.

Capability 1. is generally performed through functions of the format `someproperty(candidate; default)` where in general
`someproperty(::AbstractCandidate; default=nothing) = default`. The following such functions are currently implemented by NaiveGAflux:

* [`model(c; default)`](@ref model)  : Return a model
* [`opt(c; default)`](@ref opt)    : Return an optimizer
* [`lossfun(c; default)`](@ref lossfun) : Return a lossfunction

All such functions are obviously not used by all fitness strategies and some are used more often than others. Whether an
`AbstractCandidate` returns something other than `default` generally depends on whether it is a hyperparameter which is
being searched for or not. For example, the very simple [`CandidateModel`](@ref) has only a `model` while
[`CandidateOptModel`](@ref) has both a model and an own optimizer which may be mutated/crossedover when evolving.

Capability 2. is what is used then evolving a candidate into a new version of itself. The function to implement for new
`AbstractCandidate` types is `newcand(c::MyCandidate, mapfields)` which in most cases has the implementation
`newcand(c::MyCandidate, mapfield) = MyCandidate(map(mapfield, getproperty.(c, fieldnames(MyCandidate)))...)`.
Furthermore, candidates must also be functors from [Functors.jl](https://github.com/FluxML/Functors.jl) to
support things like GPU<->CPU movement.

Example with a new candidate type and a new fitness strategy for said type:
"""

@testset "Candidate handling" begin #src
import Functors
struct ExampleCandidate <: AbstractCandidate
    a::Int
    b::Int
end
aval(c::ExampleCandidate; default=nothing) = c.a
bval(c::ExampleCandidate; default=nothing) = c.b

Functors.@functor ExampleCandidate

struct ExampleFitness <: AbstractFitness end
NaiveGAflux._fitness(::ExampleFitness, c::AbstractCandidate) = aval(c; default=10) - bval(c; default=5)

# Ok, this is alot of work for quite little in this dummy example.
@test fitness(ExampleFitness(), ExampleCandidate(4, 3)) === 1

ctime, examplemetric = fitness(TimeFitness(ExampleFitness()), ExampleCandidate(3,1))
@test examplemetric === 2
@test ctime > 0
end #src


text:
function convolve(A::SubArray{T,n}, B::Array{T,n}, dims) where {T,n}
    common_size = tuple(map(max, size(A), size(B))...)
    pA = zeros(Complex{T}, common_size)
    pB = zeros(Complex{T}, common_size)
    rangesA = [1:x for x in size(A)]
    rangesB = [1:x for x in size(B)]
    pA[rangesA...] = A
    pB[rangesB...] = B

    real(ifft!(fft!(pA, dims) .* fft!(pB, dims), dims))
end

function convolve_ComplexFloat64(A, B, dims)
    common_size = tuple(map(max, size(A), size(B))...)
    pA = zeros(Complex{Float64}, common_size)
    pB = zeros(Complex{Float64}, common_size)
    rangesA = [1:x for x in size(A)]
    rangesB = [1:x for x in size(B)]
    pA[rangesA...] = A
    pB[rangesB...] = B

    real(ifft!(fft!(pA, dims) .* fft!(pB, dims), dims))
end

# got rid of dims argument
# using real to complex fft and ifft
function convolve_Float64(A, B)
    common_size = tuple(map(max, size(A), size(B))...)
    pA = zeros(Float64, common_size)
    pB = zeros(Float64, common_size)
    rangesA = [1:x for x in size(A)]
    rangesB = [1:x for x in size(B)]
    pA[rangesA...] = A
    pB[rangesB...] = B

    irfft(rfft(pA) .* rfft(pB), common_size[1])
end

function valid_convolve(A, B)
    ranges = [min(a, b):max(a, b) for (a, b) in zip(size(A), size(B))]
    convolve_Float64(A, B)[ranges...]
end

function cumsum2(A::Array{T,ndim}) where {T,ndim}
    # cumulative sum in two dimensions
    if ndim != 2
        throw(ArgumentError("input must be two-dimensional array"))
    end
    B = similar(A)
    # first row and column are 1D cumulative sums
    B[:, 1] = cumsum(A[:, 1], 1)
    B[1, :] = cumsum(A[1, :], 2)      # B[1,1] is redundantly computed twice
    # compute rest of matrix from recursion
    (m, n) = size(A)
    if m > 1 && n > 1
        for j = 2:n, i = 2:m
            B[i, j] = A[i, j] + B[i-1, j] + B[i, j-1] - B[i-1, j-1]
        end
    end
    B
end

function optimize_normxcorr2(img)
    p = plan_rfft(img, flags = FFTW.MEASURE)
    q = plan_irfft(rfft(img), flags = FFTW.MEASURE, size(img, 1))
    return Void
end

# extension:
# Params_session.jl: optimize_all_cores(params::Params)
function optimize_all_cores(img_d::Int64)
    img = rand(img_d, img_d)
    @sync begin
        for p = 1:num_procs
            @async begin
                remotecall_fetch(p, optimize_normxcorr2, img)
            end
        end
    end
    return Void
end

"""
`NORMXCORR2` - normalized cross correlation
slide template across img, compute Pearson correlation coefficient for each
template location result has same size as MATLAB-style 'valid' convolution
efficient algorithm of J. P. Lewis
http://scribblethink.org/Work/nvisionInterface/nip.html

` correlogram = normxcorr2(template,img)`

return NaN for image patch or template with zero variance
Pearson correlation coefficient is undefined in this case

need some argument checking
e.g. sizes of template should be less than those of img
this works for arrays.  extend to Image defined in Holy's package?
"""
function normxcorr2(template, img)
    # sufficient to subtract mean from just one variable
    dt = template - mean(template)
    templatevariance = sum(dt.^2)
    numerator = valid_convolve(img, dt[end:-1:1, end:-1:1])

    ##### local statistics of img
    # zero pad image in first row and column
    # so that cumulative sums will have zeros in the same place
    (m1, m2) = size(img)
    imgpad = zeros(m1 + 1, m2 + 1)
    imgpad[2:end, 2:end] = img
    # define four combinations of Small and Large ranges
    (n1, n2) = size(template)
    if templatevariance == 0
        return zeros(m1 - n1 + 1, m2 - n2 + 1) * NaN
    end
    LL = Any[1+(n1:m1), 1+(n2:m2)]
    SL = LL - [n1; 0]
    LS = LL - [0; n2]
    SS = LL - [n1; n2]
    # sum of img and its square in template-sized neighborhoods
    s = cumsum2(imgpad)
    localsum = s[LL...] - s[SL...] - s[LS...] + s[SS...]
    s2 = cumsum2(imgpad.^2)
    localsum2 = s2[LL...] - s2[SL...] - s2[LS...] + s2[SS...]
    localvariance = localsum2 - localsum.^2 / prod(size(template))
    # localvariance is zero for image patches that are constant
    # leading to undefined Pearson correlation coefficient
    # should only be negative due to roundoff error
    localvariance[localvariance.<=0] *= NaN
    denominator = sqrt(localvariance * templatevariance)

    numerator ./ denominator
end


text:
# Gurobi environment and other supporting facilities

mutable struct Env
    ptr_env::Ptr{Cvoid}

    function Env()
        a = Ref{Ptr{Cvoid}}()
        ret = @grb_ccall(loadenv, Cint, (Ptr{Ptr{Cvoid}}, Ptr{UInt8}),
            a, C_NULL)
        if ret != 0
            if ret == 10009
                error("Invalid Gurobi license")
            else
                error("Failed to create environment (error $ret).")
            end
        end
        env = new(a[])
        # finalizer(env, free_env)  ## temporary disable: which tends to sometimes caused warnings
        env
    end
end

Base.unsafe_convert(ty::Type{Ptr{Cvoid}}, env::Env) = env.ptr_env::Ptr{Cvoid}

function is_valid(env::Env)
    env.ptr_env != C_NULL
end

function free_env(env::Env)
    if env.ptr_env != C_NULL
        @grb_ccall(freeenv, Nothing, (Ptr{Cvoid},), env.ptr_env)
        env.ptr_env = C_NULL
    end
end

function get_error_msg(env::Env)
    @assert env.ptr_env != C_NULL
    sz = @grb_ccall(geterrormsg, Ptr{UInt8}, (Ptr{Cvoid},), env.ptr_env)
    unsafe_string(sz)
end

# error

mutable struct GurobiError <: Exception
    code::Int
    msg::String

    function GurobiError(env::Env, code::Integer)
        new(convert(Int, code), get_error_msg(env))
    end
end


text:
functor(T, x) = (), _ -> x
functor(x) = functor(typeof(x), x)

functor(::Type{<:Tuple}, x) = x, y -> y
functor(::Type{<:NamedTuple}, x) = x, y -> y

functor(::Type{<:AbstractArray}, x) = x, y -> y
functor(::Type{<:AbstractArray{<:Number}}, x) = (), _ -> x

function makefunctor(m::Module, T, fs = fieldnames(T))
  @eval m begin
    $Functors.functor(::Type{<:$T}, x) = ($([:($f=x.$f) for f in fs]...),), y -> $T(y...)
  end
end

function functorm(T, fs = nothing)
  fs == nothing || isexpr(fs, :tuple) || error("@functor T (a, b)")
  fs = fs == nothing ? [] : [:($(map(QuoteNode, fs.args)...),)]
  :(makefunctor(@__MODULE__, $(esc(T)), $(fs...)))
end

macro functor(args...)
  functorm(args...)
end

isleaf(x) = functor(x)[1] === ()

function fmap1(f, x)
  func, re = functor(x)
  re(map(f, func))
end

function fmap(f, x; cache = IdDict())
  haskey(cache, x) && return cache[x]
  cache[x] = isleaf(x) ? f(x) : fmap1(x -> fmap(f, x, cache = cache), x)
end


text:
using LinearAlgebra
using YAML

parameters = YAML.load(open("cfg_default.yml"))

mutable struct Weights

    e::Array{Array{Float64,2}}
    ie::Array{Any}
    ei::Array{Any}
    de::Array{Array{Float64,2}}

    function Weights(arch::Array{Int64}, param::Dict = parameters)
        # exc
        e = Array{Float64,2}[]
        de = Array{Float64,2}[]
        for l in 1:length(arch)-1
            push!(e, clamp.(randn(arch[l],arch[l+1]) .* param["std"] .+ param["m"], 0, param["wmax"]))
            push!(de, zeros(arch[l],arch[l+1]))
        end
        # inh
        ie = Array{Any}[]
        ei = Array{Any}[]
        for l in 1:length(arch)
            if l>1 && l!=length(arch)
                push!(ie, param["wie"] .* (ones(arch[l],arch[l])-I(arch[l])))
                push!(ei, param["wei"] .* I(arch[l]))
            else
                push!(ie, [[] for _ in 1:arch[l]])
                push!(ei, [[] for _ in 1:arch[l]])
            end
        end
        new(e,ie,ei,de)
    end
end

mutable struct Ind

    exc::Array{Int64}
    inh::Array{Int64}
    el::Array{Array{Int64}}
    il::Array{Array{Int64}}

    function Ind(arch)
        n_exc = sum(arch)
        n_inh = sum(arch[2:end-1])
        n_neurons = n_exc + n_inh
        exc = collect(1:n_exc)
        inh = collect(n_exc+1:n_neurons)

        architr = [0; arch]
        el = Array{Int64}[]
        il = Array{Int64}[]

        for l in 2:length(architr)
            push!(el, collect(sum(architr[1:l-1])+1:sum(architr[1:l])))
            if l > 2 && l!=length(architr)
                push!(il, collect(sum(architr[2:l-1])+1:sum(architr[2:l])) .+ n_exc .- arch[1])
            else
                push!(il, [])
            end
        end
        new(exc,inh,el,il)
    end
end

mutable struct Network

    n_exc::Int64
    n_inh::Int64
    n_neurons::Int64

    i::Ind

    v::Array{Float64}
    u::Array{Float64}
    d::Array{Float64}
    a::Array{Float64}

    s::Weights
    trace_e::Array{Float64}
    ho::Array{Float64}

    I::Array{Float64}
    da::Float64
    arch::Array{Int64}
    param::Dict

    function Network(arch::Array{Int64}, param::Dict = parameters)
        n_exc = sum(arch)
        n_inh = sum(arch[2:end-1])
        n_neurons = n_exc + n_inh

        i = Ind(arch)

        v = param["c"]*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? param["de"] : param["di"] for i in 1:n_neurons]
        a = [(i<=n_exc) ? param["ae"] : param["ai"] for i in 1:n_neurons]

        s = Weights(arch,param)
        trace_e = zeros(n_exc)
        ho = zeros(n_neurons)

        I = zeros(n_neurons)
        da = 0.0
        arch = arch

        new(n_exc,n_inh,n_neurons,i,v,u,d,a,s,trace_e,ho,I,da,arch,param)
    end
end

# %% functions

# direct spike of the neuron
function input!(net::Network, input_spikes::Array{Int64})
    net.v[input_spikes] .= 40.0
end

# superthreshold current processed by the neuron
function input!(net::Network, input_spikes::Array{Tuple{Int64,Float64}})
    #net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    net.I = zeros(net.n_neurons)
    max_current = net.param["imax"]
    min_current = net.param["imin"]
    indices = [input_spikes[i][1] for i in eachindex(input_spikes)]
    intensity = [input_spikes[i][2] for i in eachindex(input_spikes)]
    net.I[indices] .+= (max_current-min_current) .* intensity .+ min_current
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
end

function spike!(net::Network)
    # spikes
    net.I = zeros(net.n_neurons)
    spiked = findall(x->x>net.param["vthresh"],net.v .- net.ho)
    spiked_e = filter(x->(x in net.i.exc), spiked)
    spiked_ho = filter(x-> x > net.arch[1], spiked_e)
    net.ho[spiked_ho] .+= net.param["theta"]
    net.v[spiked] .= net.param["c"]
    net.u[spiked] .+= net.d[spiked]
    # increment current
    for neuron in spiked
        for l in 1:length(net.arch)-1
            if neuron in net.i.el[l]
                neur = neuron - (net.i.el[l][1] - 1)
                net.I[net.i.el[l+1]] .+= net.s.e[l][neur,:]
                net.I[net.i.il[l]] .+= net.s.ei[l][neur,:]
            end
            if neuron in net.i.il[l]
                neur = neuron - (net.i.il[l][1] - 1)
                net.I[net.i.el[l]] .+= net.s.ie[l][neur,:]
            end
        end
    end

    #time step (Izhikevic model)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.u .= @. net.u+net.a*(0.2*net.v-net.u)
    net.trace_e[spiked_e] .= 0.1

    # # LTP & LTD (STDP)
    # for neuron in spiked_e
    #     for l in 1:length(net.arch)
    #         if neuron in net.i.el[l]
    #             neur = neuron - (net.i.el[l][1] - 1)
    #             if l > 1
    #                 net.s.de[l-1][:,neur] .+= net.param["apost"] .* net.trace_e[net.i.el[l-1]]
    #             end
    #             if l < length(net.arch)
    #                 net.s.de[l][neur,:] .-= net.param["apre"] .* net.trace_e[net.i.el[l+1]]
    #             end
    #         end
    #     end
    # end

    # time decay
    net.ho .*= net.param["ttheta"]
    net.trace_e .*= 0.95
    net.da *= net.param["tda"]
    spiked
end

function learn!(net::Network)
    #increase the synaptic weigths
    for l in 1:length(net.arch)-1
        net.s.e[l] .+= (0.002+net.da) .* net.s.de[l]
        net.s.e[l] .= clamp.(net.s.e[l], 0, net.param["wmax"])
        net.s.de[l] .*= 0.99
    end
end

# %% main functions

function step!(net::Network, train::Bool=false)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Int64}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

# %% visualization function

function weight_viz(net::Network, method::String = "histogram")
    weights = net.s.e
    if method == "histogram"
        h = []
        for i in 1:length(net.arch)-1
            vec = filter(x-> x>=0 ,weights[i])
            h_tmp = histogram(vec,
                          legend = false,
                          xlabel = "Synaptic strength layer $(i)",
                          ylabel = "Synapse number")
            push!(h, h_tmp)
        end
        display(plot(h..., layout = length(net.arch)-1))

    elseif method == "heatmap"
        h = []
        for i in 1:length(net.arch)-1
            h_tmp = heatmap(weights[i],
                            xlabel = "connection matrix layer $(i)")
            push!(h, h_tmp)
        end
        display(plot(h..., layout = length(net.arch)-1))

    end
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
#########
# TRACE #
#########

function update!(tr::OptimizationTrace{T,O}, state::S,
                 iteration::Integer, val::T, dt::Dict,
                 store_trace::Bool, show_trace::Bool, show_every::Int = 1,
                 callback = nothing) where {S<:AbstractOptimizerState, T, O}
    os = OptimizationTraceRecord{T,O}(iteration, val, dt)
    store_trace && push!(tr, os)
    if show_trace
        if iteration % show_every == 0
            show(os)
            print("\n")
            flush(stdout)
        end
    end
    if callback !== nothing && (iteration % show_every == 0)
        if store_trace
            stopped = callback(tr)
        else
            stopped = callback(os)
        end
    else
        stopped = false
    end
    stopped
end

function trace!(tr, iteration, objfun, state, population, method, options, curr_time=time())
    dt = Dict{String,Any}()
    dt["time"] = curr_time
    # set additional trace value
    trace!(dt, objfun, state, population, method, options)
    update!(tr,
            state,
            iteration,
            value(state),
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

"""
    trace!(record::Dict{String,Any}, objfun, state, population, method, options)

Update the trace `record`. This function allows to supplement an additional information into the optimization algorithm trace by modifying a trace `record`. It can be overridden by specifying particular parameter types.
"""
trace!(record::Dict{String,Any}, objfun, state, population, method, options) = ()

##############
# POPULATION #
##############

"""
    initial_population(method, individual::AbstractVector)

Initialize population by replicating the `individual` vector.
"""
initial_population(method::M, individual::I; kwargs...) where {M<:AbstractOptimizer, I<:AbstractVector} =
    [copy(individual) for i in 1:population_size(method)]

"""
    initial_population(method, individuals::AbstractVector{<:AbstractVector})

Initialize population from the collection of `individuals` vectors.
"""
function initial_population(method::M, individuals::AbstractVector{I}; kwargs...) where {M<:AbstractOptimizer, I<:AbstractVector}
    n = population_size(method)
    @assert length(individuals) ==  n "Size of initial population must be $n"
    individuals
end

"""
    initial_population(method, individual::Function)

Initialize population from the `individual` function which returns an individual object.
"""
initial_population(method::M, individualFunc::Function; kwargs...) where {M<:AbstractOptimizer} =
    [individualFunc() for i in 1:population_size(method)]

"""
    initial_population(method, individual::AbstractMatrix)

Initialize population by replicating the `individual` matrix.
"""
function initial_population(method::M, individual::I; kwargs...) where {M<:AbstractOptimizer, I<:AbstractMatrix}
    [copy(individual) for i in 1:population_size(method)]
end

"""
    initial_population(method, bounds::ConstraintBounds)

Initialize a random population within the individual `bounds`.
"""
function initial_population(method::M, bounds::ConstraintBounds;
                            rng::AbstractRNG=default_rng()
                           ) where {M<:AbstractOptimizer}
    n = population_size(method)
    cn = nconstraints_x(bounds)
    indv = rand(rng, cn, n)
    if length(bounds.eqx) > 0
        indv[bounds.eqx,:] .= bounds.valx
    end
    T = eltype(bounds)
    rngs = Dict(i=>zeros(T,2) for i in bounds.ineqx)
    for (j,v,s) in zip(bounds.ineqx, bounds.bx, bounds.σx)
        rngs[j][1] -= s*v
        if s > 0
            rngs[j][2] += v
        end
    end
    for i in 1:n
        for (j,(r,l)) in rngs
            indv[j,i] *= r
            indv[j,i] += l
        end
    end
    initial_population(method,  [collect(i) for i in  eachcol(indv)])
end

########
# MISC #
########

default_values(x::AbstractArray{T}) where {T} = fill!(similar(x), zero(T))
default_values(x::AbstractArray{T}) where {T<:AbstractFloat} = fill!(similar(x), T(NaN))

function funargnum(f)
    fobj = try
        first(methods(f))
    catch
        @error "Cannot find function $f"
    end
    fobj.nargs
end


text:
include("HtmlParser.jl")

function htmljson_serializer(json::HtmlJSON)::String
    html = string('<', json.tag)
    for (attr, content) in json.attributes
        html = string(html, " ", attr)
        type = typeof(content)
        if type == String
            if content != ""
                html = string(html, "=\"", content, '\"')
            end
        else
            html = string(html, '=', content)
        end
    end
    if length(json.children) == 0
        html = string(html, "/>")
    else
        html = string(html, '>')
        for child in json.children
            type = typeof(child)
            if type == HtmlJSON
                html = string(html, htmljson_serializer(child))
            elseif type == String
                html = string(html, child)
            end
        end
        html = string(html, "</", json.tag, '>')
    end
    return html
end

max_stars_repo_name: haverzard/Husk

max_issues_repo_name: haverzard/Husk

max_forks_repo_name: haverzard/Husk


text:
using BayesLaplace
using Test

@testset "BayesLaplace.jl" begin
    @testset "Curvature" begin
        include("curvature.jl")
    end

    @testset "Laplace" begin
        include("laplace.jl")
    end
end


text:
# P(n) = (1/2)*n*(sin(2*pi/n))

P(n) = pi

function c(n = BigFloat)
    if(n==2)
        return 0
    else
        return BigFloat(sqrt((1/2)*(1+c(n-1))))
    end
end

function s(n = BigFloat)
    if(n==2)
        return 1
    else
        return BigFloat(BigFloat(sqrt((1/2)*(1-c(n-1)))))    end
end

function s2(n = BigFloat)
    if (n==2)
        return 1
    else
        return (BigFloat(s2(2*n)))/(2*(BigFloat(c(n))))
    end
end

A(n = BigFloat) =  BigFloat((2^(n-1))*s(n))

println("ZADANIE 5")

setprecision(128) do
    for i in 2:128
        print(i)
        print("  ")
        print(P(2^i))
        print("  ")
        print(A(i))
        print("  ")
        println(abs(P(2^i) - A(i)))
    end
end

max_stars_repo_path: II rok/I semestr/ANM/Cwiczenia 02/Zadanie5.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
@testset "ByteLevel" begin
  bbpe = ByteLevelBPE(Base.download("https://s3.amazonaws.com/models.huggingface.co/bert/gpt2-merges.txt"))
  weird_sentence_answer = ["ðĿ", "ĺ", "Ĳ", "Ġ", "ðĿ", "Ĺ", "Ĩ", "ðĿ", "Ĺ", "Ĥ", "ðĿ", "Ĺ", "´", "ðĿ", "ĺ", "©", "ðĿ", "Ļ", "©", "Ġ", "ðĿ", "ļ", "ķ", "ðĿ", "ļ", "ĺ", "ðĿ", "ļ", "ĺ", "k", "Ġ", "ðĿ", "Ļ", "ĩ", "ðĿ", "ĺ", "ª", "ðĿ", "Ĺ", "¸", "ðĿ", "ĸ", "¾", "Ġ", "ðĿ", "ĸ", "Ĩ", "Ġ", "ðĿ", "ĸ", "ĵ", "ðĿ", "ĸ", "Ķ", "ðĿ", "ĸ", "Ĺ", "ðĿ", "ķ", "ŀ", "ðĿ", "ķ", "Ĵ", "ðĿ", "ķ", "Ŀ", "Ġ", "ðĿ", "Ķ", "°", "ðĿ", "ĵ", "®", "ðĿ", "ĵ", "ĥ", "ðĿ", "Ĵ", "ķ", "ðĿ", "ĳ", "Ĵ", "ðĿ", "Ĳ", "§", "ce"]
  @test bbpe("𝘐 𝗆𝗂𝗴𝘩𝙩 𝚕𝚘𝚘k 𝙇𝘪𝗸𝖾 𝖆 𝖓𝖔𝖗𝕞𝕒𝕝 𝔰𝓮𝓃𝒕𝑒𝐧ce") == weird_sentence_answer

  bbpe1 = GenericBPE{String}(; input_transform = BytePairEncoding.gpt2_tokenizer, codemap = BytePairEncoding.default_codemap(), merging_rank=Dict{Tuple{String, String}, Int}(("\u0120", "l")=>1, ("\u0120l", "o")=>2, ("\u0120lo", "w")=>3, ("e", "r")=>4))
  @test bbpe1(" lower newer") == ["\u0120low", "er", "\u0120", "n", "e", "w", "er"]
end


text:
module Resource
# This module handles all the handling, parsing and transforming HTTP commands,
# and brokers between HTTP requests and Service that fullfils them.
using HTTP, JSON3, Dates, LoggingExtras, Sockets, ThreadPools
using ..DFControl, ..Service
using ..Utils
using ..Servers
using ..Jobs
using ..Calculations
using ..Structures
using ..FileIO

const ROUTER = HTTP.Router()

include("job.jl")
include("fileio.jl")

# GENERAL
kill_server(req) = exit()
HTTP.@register(ROUTER, "PUT", "/kill_server", kill_server)

get_server_config(req) = Service.server_config()
HTTP.@register(ROUTER, "GET", "/server_config", get_server_config)

get_ispath(req) = ispath(job_path(req))
HTTP.@register(ROUTER, "GET", "/get_ispath/*", get_ispath)

get_readdir(req) = readdir(job_path(req))
HTTP.@register(ROUTER, "GET", "/readdir/*", get_readdir)

# PSEUDOS

function pseudos(req)
    fuzzy = String(req.body)
    fuzzy = fuzzy == "\"\"" ? "" : fuzzy
    return Service.pseudos(splitpath(req.target)[end], fuzzy)
end
HTTP.@register(ROUTER, "GET", "/pseudos/*", pseudos)

pseudo_sets(req) = Service.pseudo_sets()
HTTP.@register(ROUTER, "GET", "/pseudos", pseudo_sets)

configure_pseudoset(req) = Service.configure_pseudoset(JSON3.read(req.body,String), job_path(req))
HTTP.@register(ROUTER, "POST", "/configure_pseudoset/*", configure_pseudoset)

rm_pseudos!(req) = Service.rm_pseudos!(JSON3.read(req.body, String))
HTTP.@register(ROUTER, "PUT", "/rm_pseudos", rm_pseudos!)
# EXECS
verify_exec(req) = Service.verify_exec(JSON3.read(req.body, Exec))
HTTP.@register(ROUTER, "GET", "/verify_exec", verify_exec)

known_execs(req) = (d = JSON3.read(req.body); Service.known_execs(d["exec"],d["dir"]))
HTTP.@register(ROUTER, "GET", "/known_execs/", known_execs)

get_exec(req) = Service.load_exec(splitpath(req.target)[end])
HTTP.@register(ROUTER, "GET", "/exec/*", get_exec)

save_exec(req) = Service.save(JSON3.read(req.body, Exec))
HTTP.@register(ROUTER, "POST", "/exec/", save_exec)

add_environment(req) = Service.add_environment(JSON3.read(req.body,Jobs.Environment), splitpath(req.target)[end])
HTTP.@register(ROUTER, "POST", "/environment/*", add_environment)

get_environment(req) = Service.get_environment(splitpath(req.target)[end])
HTTP.@register(ROUTER, "GET", "/environment/*", get_environment)

rm_environment!(req) = Service.rm_environment!(splitpath(req.target)[end])
HTTP.@register(ROUTER, "PUT", "/environment/*", rm_environment!)

# RUNNING
function requestHandler(req)
    start = Dates.now()
    @info (timestamp = start, event = "ServiceRequestBegin", tid = Threads.threadid(),
           method = req.method, target = req.target)
    local resp
    try
        obj = HTTP.handle(ROUTER, req)
        if obj === nothing
            resp = HTTP.Response(204)
        else
            resp = HTTP.Response(200, JSON3.write(obj))
        end
    catch e
        s = IOBuffer()
        showerror(s, e, catch_backtrace(); backtrace = true)
        errormsg = String(resize!(s.data, s.size))
        @error errormsg
        resp = HTTP.Response(500, errormsg)
    end
    stop = Dates.now()
    @info (timestamp = stop, event = "ServiceRequestEnd", tid = Threads.threadid(),
           method = req.method, target = req.target, duration = Dates.value(stop - start),
           status = resp.status, bodysize = length(resp.body))
    return resp
end

function run()
    s = Server("localhost")
    cd(s.default_jobdir)
    Service.global_logger(Service.daemon_logger())
    port, server = listenany(ip"0.0.0.0", 8080)
    s.port = port
    Servers.save(s)
    @tspawnat 2 HTTP.serve(requestHandler, "0.0.0.0", port, server=server)
    with_logger(Service.daemon_logger()) do
        Service.main_loop(s)
    end
    close(server)
    return
    end
end


text:
# ---
# title: 1641. Count Sorted Vowel Strings
# id: problem1641
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/count-sorted-vowel-strings/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _the number of strings of length_`n` _that
# consist only of vowels (_`a` _,_`e` _,_`i` _,_`o` _,_`u` _) and are
# **lexicographically sorted**._
#
# A string `s` is **lexicographically sorted** if for all valid `i`, `s[i]` is
# the same as or comes before `s[i+1]` in the alphabet.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: 5
#     Explanation: The 5 sorted strings that consist of vowels only are ["a","e","i","o","u"].
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: 15
#     Explanation: The 15 sorted strings that consist of vowels only are
#     ["aa","ae","ai","ao","au","ee","ei","eo","eu","ii","io","iu","oo","ou","uu"].
#     Note that "ea" is not a valid string since 'e' comes after 'a' in the alphabet.
#
#
# **Example 3:**
#
#
#
#     Input: n = 33
#     Output: 66045
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 50`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# File for Activation Function

# Sigmoid Functiom

Sigmoid(x::Float64) = 1/(1+exp(-x))

# ReLU Function

ReLU(x::Float64) = max(x,0)

# Identity Function

identity(x::Float64) = x


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

using Markdown
using InteractiveUtils

# ╔═╡ 194baef2-0417-11ec-05ab-4527ef614024
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 38c9ced5-dcf8-4e03-ac07-7c435687861b
begin
	using DoublePendulum
	using Plots
end

# ╔═╡ 396b5d7a-a7a4-4f22-a87e-39b405e8d62a
md"
# Double Pendulum

JingYu Ning
"

# ╔═╡ 2a606ecf-acf0-41ad-9290-7569dbb22b5a
md"
The data is provided by [IBM](https://developer.ibm.com/exchanges/data/all/double-pendulum-chaotic/)

> In this dataset, videos of the double pendulum were taken using a high-speed Phantom Miro EX2 camera. To make the extraction of the arm positions easier, a matte black background was used, and the three datums were marked with red, green and blue fiducial markers. The camera was placed at 2 meters from the pendulum, with the axis of the objective aligned with the first pendulum datum. The pendulum was launched by hand, and the camera was motion triggered. The dataset was generated on the basis of 21 individual runs of the pendulum. Each of the recorded sequences lasted around 40s and consisted of around 17500 frames.
"

# ╔═╡ 5268feee-bda2-4612-9d4c-a1db424a11c7
data_x, data_y = DoublePendulum.preprocess(DoublePendulum.get_data(i=20));

# ╔═╡ 4d0b08a4-8a54-41fd-997f-ad54d4c984cd
m = DoublePendulum.get_model();

# ╔═╡ ad6302b2-3d62-4a3f-b8bf-f69bab80c7a4
ground_truth_data = cat(
	[data_x[:, :, :, i]
	for i in 1:size(data_x, 3):size(data_x)[end]]..., dims=3
)[1, :, :];

# ╔═╡ 794374ce-6674-481d-8a3b-04db0f32d233
begin
	n = 5
	inferenced_data = data_x[:, :, :, 1:1]
	for i in 1:n
		inferenced_data = cat(
			inferenced_data,
			m(inferenced_data[:, :, :, i:i]),
			dims=4
		)
	end

	inferenced_data = cat(
		[inferenced_data[:, :, :, i] for i in 1:n]..., dims=3
	)[1, :, :]
end;

# ╔═╡ 9c8b3f8a-1b85-4c32-a416-ead51b244b94
begin
	c = [
		RGB([239, 71, 111]/255...),
		RGB([6, 214, 160]/255...),
		RGB([17, 138, 178]/255...)
	]
	xi, yi = [2, 4, 6], [1, 3, 5]

	anim = @animate for i in 1:size(inferenced_data)[end]
		i_data = [0, 0, inferenced_data[:, i]...]
		g_data = [0, 0, ground_truth_data[:, i]...]

		scatter(
			legend=false, ticks=false,
			xlim=(-1500, 1500), ylim=(-1500, 1500), size=(400, 350)
		)
		plot!(i_data[xi], i_data[yi], color=:black)
		scatter!(i_data[xi], i_data[yi], color=c, markersize=8)
		plot!(g_data[xi], g_data[yi], color=:gray)
		scatter!(g_data[xi], g_data[yi], color=c, markersize=4)

		if i ≤ 30
			annotate!(-1400, -1400, text("t=$i", :left, color=:black))
		else
			annotate!(-1400, -1400, text("t=$i", :left, color=:red))
		end
	end

	gif(anim)
end

# ╔═╡ Cell order:
# ╟─396b5d7a-a7a4-4f22-a87e-39b405e8d62a
# ╟─2a606ecf-acf0-41ad-9290-7569dbb22b5a
# ╟─194baef2-0417-11ec-05ab-4527ef614024
# ╠═38c9ced5-dcf8-4e03-ac07-7c435687861b
# ╠═5268feee-bda2-4612-9d4c-a1db424a11c7
# ╠═4d0b08a4-8a54-41fd-997f-ad54d4c984cd
# ╠═ad6302b2-3d62-4a3f-b8bf-f69bab80c7a4
# ╠═794374ce-6674-481d-8a3b-04db0f32d233
# ╟─9c8b3f8a-1b85-4c32-a416-ead51b244b94


text:
include("../examples/lorenz.jl")
using JLD
function p_y_g_x(a, σ_o)
        term =  log(1/sqrt(2π)/σ_o)
        return sum(term .- 0.5*a.*a/σ_o/σ_o)
end
function obs(x, σ_o)
    return x .+ σ_o*randn(d)
end
function obs1(x, σ_o)
	return [x[3] + σ_o*randn()]
end

function resample(x, w)
    # multinomial resampling
    cdfw = cumsum(w)
    new_pts = similar(x)
	N_p = size(w)[1]
    for j = 1:N_p
        r = rand()
        for i = 1:N_p
            if cdfw[i] >= r
                new_pts[:,j] .= x[:,i]
                break
            end
        end
    end
    return new_pts
end

function sir(y,x_ip,Δ,obs_fun,σ_o,σ_d,N_thr=10)
    K = size(y)[2]
    N_p = size(x_ip)[2]

    N_y = K*Δ
	@show N_y
	w_trj = zeros(N_p, N_y)
    x_trj = zeros(d, N_p, N_y)

    w = ones(N_p)./N_p
    x = copy(x_ip)
	count = 0
	j = 1
    for k = 1:N_y
        if k % Δ == 0
		    for i = 1: N_p
                logwi = log(w[i]) + p_y_g_x(y[:,j] .- obs_fun(x[:,i],σ_o), σ_o)
                w[i] = exp(logwi)
			end
            w ./= sum(w)
            N_eff = 1.0/sum(w.*w)
            if (N_eff < N_thr)
				count = count + 1
			    x .= resample(x, w)
		    end
			j = j + 1
		end
        for i = 1:N_p
            x[:,i] .= next(x[:,i],σ_d)
        end
        x_trj[:,:,k] .= x
        w_trj[:,k] .= w
    end
	@show "Resampling was triggered ", count, " times out of ", N_y
    return x_trj, w_trj
end
"""
    assimilate(K, Np, σ_o, σ_d,
				   Δ, Nth)

Perform `K` data assimilation steps using an SIR algorithm with `Np` particles. Other inputs:

    1. `σ_o`: std of Gaussian observation noise
	2. `σ_d`: std of Gaussian dynamics noise, which is added to every component and at every timestep.
	3. `Δ`: inter-observation number of timesteps, e.g., if `Δ = 5,` observations are assumed available at timestep `0,5,10,...,5(K-1)`.
	4. `Nth`: number of particles below which to resample
	5. `obs`: observation map, e.g., `obs(x) = x[3]`

Outputs:

    1. `x`: orbit of particles. size: `dX x Np x (Δ K)`
	2. `w`: orbit of weights. size: `Np x (Δ K)`
	3. `y`: synthetic observations. size: `dY x K`
	4. `x_true`: ``true'' orbit that generated the observations. size: `dX x (Δ K)`

# Examples
```julia-repl
julia> x, w, y, x_true = assimilate(500, 1000, 0.1, 0.1, 1, 10, obs)
```
"""
function assimilate(K, Np, σ_o, σ_d,
				   Δ, Nth, obsfun)
    x = rand(d,Np)
	Ny = K*Δ

	ytest = obsfun(rand(d),σ_o)
	dy = size(ytest)[1]
    x_true = ones(d,Ny)
    x0_true = rand(d)
    Nrunup = 2000
    for i = 1:Nrunup
        x0_true = next(x0_true, 0.0)
    end
    for i = 1:Nrunup
        for j = 1:Np
            x[:,j] = next(x[:,j],σ_d)
        end
    end
    x_true[:,1] .= x0_true
    y = zeros(dy, K)
	y[:,1] .= obsfun(x0_true, σ_o)
	k = 1
    for i = 2:Ny
        x_true[:,i] = next(x_true[:,i-1], σ_d)
		if i % Δ == 0
        	y[:,k] .= obsfun(x_true[:,i], σ_o)
			k = k+1
		end
    end
    # store trajectory of w and x.
    x_trj, w_trj = sir(y, x, Δ, obsfun, σ_o, σ_d, Nth)
    return x_trj, w_trj, y, x_true
end


text:
include("utils/main.jl")

server_start()

@testset "Space before Header colon" begin
    r = HTTP.request("GET", "http://127.0.0.1:8000/", ["a " => "b"])
    parse_and_test_request(r, 200, ["a" => "b"], 230, false, true)
end

@testset "Space after Header colon" begin
    r = HTTP.request("GET", "http://127.0.0.1:8000/", ["a" => " b"])
    parse_and_test_request(r, 200, ["a" => "b"], 230, false, true)
end

server_close()


text:
using Statistics:mean, median

# TODO: Add AOM/MOA/LSCP combination strategies

"""
    combine_mean(scores_mat)

Combination method to merge outlier scores from multiple detectors using the mean value of scores.

Parameters
----------
    scores_mat::AbstractMatrix{T}
A matrix, with each row representing the scores for a specific instance and each column representing a detector.

Returns
----------
    combined_scores::AbstractVector{T}
    The combined scores, i.e. column-wise mean.

Examples
----------
    scores = [1 2; 3 4; 5 6]
    combine_mean(scores) # [1.5, 3.5, 5.5]
"""
combine_mean(scores_mat::AbstractMatrix{<:Real}) = dropdims(mean(scores_mat, dims=2), dims=2)
combine_mean(scores::Scores...) = combine_mean(hcat(scores...))

"""
    combine_median(scores_mat)

Combination method to merge outlier scores from multiple detectors using the median value of scores.

Parameters
----------
    scores_mat::AbstractMatrix{T}
A matrix, with each row representing the scores for a specific instance and each column representing a detector.

Returns
----------
    combined_scores::AbstractVector{T}
The combined scores, i.e. column-wise median.

Examples
----------
    scores = [1 2; 3 4; 5 6]
    combine_median(scores) # [1.5, 3.5, 5.5]
"""
combine_median(scores_mat::AbstractMatrix{<:Real}) = dropdims(median(scores_mat, dims=2), dims=2)
combine_median(scores::Scores...) = combine_median(hcat(scores...))

"""
    combine_max(scores_mat)

Combination method to merge outlier scores from multiple detectors using the maximum value of scores.

Parameters
----------
    scores_mat::AbstractMatrix{T}
A matrix, with each row representing the scores for a specific instance and each column representing a detector.

Returns
----------
    combined_scores::AbstractVector{T}
The combined scores, i.e. column-wise maximum.

Examples
----------
    scores = [1 2; 3 4; 5 6]
    combine_max(scores) # [2, 4, 6]
"""
combine_max(scores_mat::AbstractMatrix{<:Real}) = dropdims(maximum(scores_mat, dims=2), dims=2)
combine_max(scores::Scores...) = combine_max(hcat(scores...))