text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module ZygoteGMM
include("../../shared/load.jl")
using ADPerfTest
using GMMData
using Zygote
using SpecialFunctions
using LinearAlgebra

export get_gmm_test

"Computes logsumexp. Input should be 1 dimensional"
function logsumexp(x)
    mx = maximum(x)
    log(sum(exp.(x .- mx))) + mx
end

function ltri_unpack(D, LT)
    d = length(D)
    make_col(r::Int, L) = vcat(zeros(r - 1), D[r], reshape([L[i] for i=1:d-r], d - r))
    col_start(r::Int) = (r - 1) * (2d - r) ÷ 2
    inds(r) = col_start(r) .+ (1:d-r)
    hcat([make_col(r, LT[inds(r)]) for r=1:d]...)
end

function get_Q(d, icf)
    ltri_unpack((icf[1:d]), icf[d+1:end])
end

function log_gamma_distrib(a, p)
    out = 0.25 * p * (p - 1) * 1.1447298858494002 #convert(Float64, log(pi))
    out += sum(j -> loggamma(a + 0.5 * (1 - j)), 1:p)
    out
end

Zygote.@adjoint function loggamma(x)
    loggamma(x), Δ -> (Δ * digamma(x),)
end

function log_wishart_prior(wishart::Wishart, sum_qs, Qs, k)
    p = size(Qs, 1)
    n = p + wishart.m + 1
    C = n * p * (log(wishart.gamma) - 0.5 * log(2)) - log_gamma_distrib(0.5 * n, p)

    frobenius = sum(abs2, Qs)
	0.5 * wishart.gamma^2 * frobenius - wishart.m * sum(sum_qs) - k * C
end

function diagsums(Qs)
    mapslices(slice -> sum(diag(slice)), Qs; dims=[1, 2])
end

Zygote.@adjoint function diagsums(Qs)
    diagsums(Qs),
    function (Δ)
        Δ′ = zero(Qs)
        for (i, δ) in enumerate(Δ)
            for j in 1:size(Qs, 1)
                Δ′[j, j, i] = δ
            end
        end
        (Δ′,)
    end
end

function expdiags(Qs)
    mapslices(Qs; dims=[1, 2]) do slice
        slice[diagind(slice)] .= exp.(slice[diagind(slice)])
        slice
    end
end

Zygote.@adjoint function expdiags(Qs)
    expdiags(Qs),
    function (Δ)
        Δ′ = zero(Qs)
        Δ′ .= Δ
        for i in 1:size(Qs, 3)
            for j in 1:size(Qs, 1)
                Δ′[j, j, i] *= exp(Qs[j, j, i])
            end
        end
        (Δ′,)
    end
end

function unzip(tuples)
    map(1:length(first(tuples))) do i
        map(tuple -> tuple[i], tuples)
    end
end

Zygote.@adjoint function map(f, args...)
    ys_and_backs = map((args...) -> Zygote._forward(__context__, f, args...), args...)
    ys, backs = unzip(ys_and_backs)
    ys, function (Δ)
        Δf_and_args_zipped = map((f, δ) -> f(δ), backs, Δ)
        Δf_and_args = unzip(Δf_and_args_zipped)
        Δf = reduce(Zygote.accum, Δf_and_args[1])
        (Δf, Δf_and_args[2:end]...)
    end
end

Base.:*(::Float64, ::Nothing) = nothing

function gmm_objective(alphas, means, Qs, x, wishart::Wishart)
    d = size(x, 1)
    n = size(x, 2)
    k = size(means, 2)
    CONSTANT = -n * d * 0.5 * log(2 * pi)
    sum_qs = reshape(diagsums(Qs), 1, size(Qs, 3))
    Qs = expdiags(Qs)

    slse = 0.
    for ix=1:n
        formula(ik) = -0.5 * sum(abs2, Qs[:, :, ik] * (x[:,ix] .- means[:, ik])) + alphas[ik] + sum_qs[ik]
        terms = map(formula, 1:k)
        slse += logsumexp(terms)
    end

    CONSTANT + slse - n * logsumexp(alphas) + log_wishart_prior(wishart, sum_qs, Qs, k)
end

function zygote_J_to_packed_J(J, k, d)
    alphas = reshape(J[1], :)
    means = reshape(J[2], :)
    icf_unpacked = map(1:k) do Q_idx
        Q = J[3][:, :, Q_idx]
        lt_cols = map(1:d-1) do col
            Q[col+1:d, col]
        end
        vcat(diag(Q), lt_cols...)
    end
    icf = collect(Iterators.flatten(icf_unpacked))
    packed_J = vcat(alphas, means, icf)
    packed_J
  end

mutable struct ZygoteGMMContext
    input::Union{GMMInput, Nothing}
    Qs
    wrapper_gmm_objective::Union{Function, Nothing}
    objective::Float64
    zygote_J
end

function zygote_gmm_prepare!(ctx::ZygoteGMMContext, input::GMMInput)
    # Running Zygote.gradient with test input, because the first invocation
    # of gradient on a given function is very long.
    # Using test input ensures that all computations related to the actual input
    # are done in calculate_jacobian!
    # Test input taken from ../../../../data/gmm/test.txt
    testinput = GMMInput([-0.649014 1.181166 -0.758453], [0.092339 0.345561 0.538817; 0.18626 0.396767 0.419195], [0.586443 -1.509405 0.166813; -0.851887 0.875874 -1.965419; 0.800321 -0.24279 -1.270071], [1.175171 2.02916]', Wishart(1.0, 0))
    testd = size(testinput.x, 1)
    testk = size(testinput.means, 2)
    testQs = cat([get_Q(testd, testinput.icfs[:, ik]) for ik in 1:testk]...; dims=[3])
    test_wrapper_gmm_objective = (alphas, means, Qs) -> gmm_objective(alphas, means, Qs, testinput.x, testinput.wishart)
    Zygote.gradient(test_wrapper_gmm_objective, testinput.alphas, testinput.means, testQs)

    ctx.input = input
    d = size(input.x, 1)
    k = size(input.means, 2)
    ctx.Qs = cat([get_Q(d, input.icfs[:, ik]) for ik in 1:k]...; dims=[3])
    ctx.wrapper_gmm_objective = (alphas, means, Qs) -> gmm_objective(alphas, means, Qs, input.x, input.wishart)
end

function zygote_gmm_calculate_objective!(ctx::ZygoteGMMContext, times)
    for i in 1:times
        ctx.objective = gmm_objective(ctx.input.alphas, ctx.input.means, ctx.Qs, ctx.input.x, ctx.input.wishart)
    end
end

function zygote_gmm_calculate_jacobian!(ctx::ZygoteGMMContext, times)
    for i in 1:times
        ctx.zygote_J = Zygote.gradient(ctx.wrapper_gmm_objective, ctx.input.alphas, ctx.input.means, ctx.Qs)
    end
end

function zygote_gmm_output!(out::GMMOutput, ctx::ZygoteGMMContext)
    out.objective = ctx.objective
    d = size(ctx.input.x, 1)
    k = size(ctx.input.means, 2)
    out.gradient = ctx.zygote_J === nothing ? [] : zygote_J_to_packed_J(ctx.zygote_J, k, d)
end

get_gmm_test() = Test{GMMInput, GMMOutput}(
    ZygoteGMMContext(nothing, nothing, nothing, 0.0, nothing),
    zygote_gmm_prepare!,
    zygote_gmm_calculate_objective!,
    zygote_gmm_calculate_jacobian!,
    zygote_gmm_output!
)

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
module FameSVD
using LinearAlgebra

export fsvd

function fsvd(A::Array{R, 2}) where R<:AbstractFloat
  dim = size(A)

  V = Array{R, 2}(undef, dim[2], dim[2])
  U = Array{R, 2}(undef, dim[1], dim[2])
  S = Array{R, 1}(undef, dim[2])

  V = BLAS.syrk('U', 'T', 1.0, A)
  (S, V) = LAPACK.syevr!('V', 'A', 'U', V, 0., 0., 0, 0, 0.)

  reverse!(S)
  @inbounds for i = 1:dim[2]
    @fastmath S[i] = sqrt(S[i])
  end

  V = reverse(V; dims = 2)
  U = BLAS.gemm('N', 'N', A, V)

  # TODO check if the loop order is optimal
  @inbounds for i = 1:dim[2]
    @inbounds for j = 1:dim[1]
      U[j, i] /= S[i]
    end
  end

  return LinearAlgebra.SVD(U, S, V')
end

end # module

ext: jl


text:
using Test
using Pat

@testset "Basic Two Body Functionality" begin
    sma = 10_000.0
    ecc = 0.4
    ta = π / 4

    @test Pat.parameter(sma, ecc) == 8_400

    @testset "Apsis Radii" begin
        @test Pat.periapsis_radius(sma, ecc) == 6_000
        @test Pat.apoapsis_radius(sma, ecc) == 14_000
    end

    @test Pat.radius(sma, ecc, ta) ≈ 6.547957842623e3

    @testset "Eccentricity" begin
        @test Pat.isclosed(0.0) == true
        @test Pat.isclosed(0.3) == true
        @test Pat.isclosed(0.9) == true
        @test Pat.isclosed(1.0) == false
        @test Pat.isclosed(1.1) == false

        @test Pat.iscircular(0.0) == true
        @test Pat.iscircular(0.3) == false
        @test Pat.iscircular(0.9) == false
        @test Pat.iscircular(1.0) == false
        @test Pat.iscircular(1.1) == false

        @test Pat.iselliptical(0.0) == true
        @test Pat.iselliptical(0.3) == true
        @test Pat.iselliptical(0.9) == true
        @test Pat.iselliptical(1.0) == false
        @test Pat.iselliptical(1.1) == false

        @test Pat.isparabolic(0.0) == false
        @test Pat.isparabolic(0.3) == false
        @test Pat.isparabolic(0.9) == false
        @test Pat.isparabolic(1.0) == true
        @test Pat.isparabolic(1.1) == false

        @test Pat.ishyperbolic(0.0) == false
        @test Pat.ishyperbolic(0.3) == false
        @test Pat.ishyperbolic(0.9) == false
        @test Pat.ishyperbolic(1.0) == false
        @test Pat.ishyperbolic(1.1) == true
    end

    @test Pat.circular_velocity(12_000, 100) ≈ 1.095445115010e1
    @test Pat.mean_motion(12_000, 100) ≈ 1.095445115010e-1
    @test Pat.period(12_000, 100) ≈ 5.735737209545e1
    @test Pat.angular_momentum(12_000, sma, ecc) ≈ 1.003992031841e4
    @test Pat.flight_path_angle(sma, ecc, ta) ≈ 2.170092525064e-1

    cb = Pat.CelestialBody(12_000)
    @test Pat.circular_velocity(cb, 100) ≈ 1.095445115010e1
    @test Pat.mean_motion(cb, 100) ≈ 1.095445115010e-1
    @test Pat.period(cb, 100) ≈ 5.735737209545e1
    @test Pat.angular_momentum(cb, sma, ecc) ≈ 1.003992031841e4
end

@testset "Keplerian Elements" begin
    sma = 10_000
    ecc = 0.4
    inc = 0.5
    aop = 0.6
    raan = 0.7
    ta = 0.8

    @testset "Single Argument Constructor" begin
        @test Pat.semi_major_axis(Pat.KeplerianElements(sma)) == sma
    end

    @testset "Default Inner Constructor" begin
        ke = Pat.KeplerianElements(sma, ecc, inc, aop, raan, ta)
        @test Pat.semi_major_axis(ke) == sma
        @test Pat.eccentricity(ke) == ecc
        @test Pat.inclination(ke) == inc
        @test Pat.argument_of_periapsis(ke) == aop
        @test Pat.right_ascension(ke) == raan
        @test Pat.true_anomaly(ke) == ta
    end

    @testset "Invalid Construction" begin
        @test_throws DomainError Pat.KeplerianElements(0.0)
        @test_throws DomainError Pat.KeplerianElements(1.0; ecc=-0.4)
    end

    @testset "Derived Values" begin
        sma = 10_000.0
        ecc = 0.4
        ta = π / 4

        ke = Pat.KeplerianElements(sma; ecc=ecc, ta=ta)
        k2 = Pat.KeplerianElements(100)

        @test Pat.parameter(sma, ecc) == 8_400

        @testset "Apsis Radii" begin
            @test Pat.periapsis_radius(ke) == 6_000
            @test Pat.apoapsis_radius(ke) == 14_000
        end

        @test Pat.radius(ke) ≈ 6.547957842623e3
        @test Pat.circular_velocity(12_000, k2) ≈ 1.095445115010e1
        @test Pat.mean_motion(12_000, k2) ≈ 1.095445115010e-1
        @test Pat.period(12_000, k2) ≈ 5.735737209545e1
        @test Pat.angular_momentum(12_000, ke) ≈ 1.003992031841e4
        @test Pat.flight_path_angle(ke) ≈ 2.170092525064e-1

        cb = Pat.CelestialBody(12_000)
        @test Pat.circular_velocity(cb, k2) ≈ 1.095445115010e1
        @test Pat.mean_motion(cb, k2) ≈ 1.095445115010e-1
        @test Pat.period(cb, k2) ≈ 5.735737209545e1
        @test Pat.angular_momentum(cb, ke) ≈ 1.003992031841e4
    end
end


text:
__precompile__()

module IRLS

# Special Matrices
using Base.LinAlg: Diagonal, Hermitian, LowerTriangular, UpperTriangular
# Factorizations
using Base.LinAlg: qrfact, cholfact!
# Vector functions
using Base.LinAlg: norm

using StatsFuns: logistic
using Distributions: Distribution, Normal, Bernoulli, Binomial, Poisson

"""
    AbstractLink

Abstract type for link functions.
"""
abstract type AbstractLink end
"""
    IdentityLink <: AbstractLink

Return a contrete type of the identity link function.
"""
struct IdentityLink <: AbstractLink end
"""
    LogitLink <: AbstractLink

Return a contrete type of the logit link function.
"""
struct LogitLink <: AbstractLink end
"""
    LogLink <: AbstractLink

Return a contrete type of the log link function.
"""
struct LogLink <: AbstractLink end

"""
    canonicallink(::Distribution)

Return the canonical link function for the distribution.
"""
canonicallink(obj::Distribution) = error("canonicallink is not defined for $(typeof(obj))")
canonicallink(obj::Normal) = IdentityLink()
canonicallink(obj::Union{Bernoulli,Binomial}) = LogitLink()
canonicallink(obj::Poisson) = LogLink()

"""
    μμ′σ²(::Distribution, ::AbstractLink, ::AbstractVector{<:Real})
    μμ′σ²(::Distribution, ::AbstractLink, ::Real)

Return the value of the inverse link function, its derivative and the variance
of the distribution evaluated at the given value.
"""
μμ′σ²(::Distribution, ::AbstractLink, η::Real) = error("μμ′σ² is not defined for that combination of arguments.")

function μμ′σ²(distribution::Distribution, link::AbstractLink, η::AbstractVector{<:Real})
    m    = length(η)
    μ    = Vector{Float64}(m)
    μ′ = Vector{Float64}(m)
    σ²   = Vector{Float64}(m)
    @inbounds for idx ∈ eachindex(μ, μ′, σ², η)
        μ[idx], μ′[idx], σ²[idx] = μμ′σ²(distribution, link, η[idx])
    end
    return (μ, μ′, σ²)
end
function μμ′σ²(::Normal, ::IdentityLink, η::Real)
    return (η, one(Float64), one(Float64))
end
function μμ′σ²(::Union{Bernoulli,Binomial}, ::LogitLink, η::Real)
    g  = logistic(η)
    g′ = g * (one(g) - g)
    σ²  = g′
    return (g, g′, σ²)
end
function μμ′σ²(::Poisson, ::LogLink, η::Real)
    g  = exp(η)
    return (g, g, g)
end

function irls(A::AbstractMatrix{<:Real},
              b::AbstractVector{<:Real},
              distribution::Distribution;
              link::AbstractLink = canonicallink(distribution),
              maxit::Integer = 25)
    F = qrfact(A)
    Q = Matrix(F.Q)
    m, n = size(A)
    x = fill(zero(Float64), n)
    η = fill(zero(Float64), m)
    w = Vector{Float64}(undef, m)
    for itt ∈ 1:maxit
        μ, μ′, σ² = μμ′σ²(distribution, link, η)
        z  = η + (b - μ) ./ μ′
        w  = μ′.^2 ./ σ²
        x₀ = x
        C  = cholfact!(Hermitian(Q' * Diagonal(w) * Q)).factors
        x  = LowerTriangular(C') \ (Q' * Diagonal(w) * z)
        x  = UpperTriangular(C) \ x
        η  = Q * x
        norm(x - x₀, 2) < 1e-8 && return (F \ η, w)
    end
end
export
    irls
end


text:
using ModernHopfieldNetworks
using Test

@test hopfieldlayers([2,3],[1,1]) == 6

my(2,3)

@testset "ModernHopfieldNetworks.jl" begin
     # hopfieldlayer([1, 2,3],[1,2,1])

end


text:
function assemble_dim(iop::IntegralOperator;compress=Matrix,location=:onsurface)
    X    = target_surface(iop)
    Y    = source_surface(iop)
    pde  = iop.kernel.pde
    T    = kernel_type(iop)
    if T === SingleLayer()
        singlelayer_dim(pde,X,Y;compress,location)
    elseif T === DoubleLayer()
        doublelayer_dim(pde,X,Y;compress,location)
    elseif T === AdjointDoubleLayer()
        adjointdoublelayer_dim(pde,X,Y;compress,location)
    elseif T === HyperSingular()
        hypersingular_dim(pde,X,Y;compress,location)
    else
        notimplemented()
    end
end

function single_doublelayer_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
    msg = "unrecognized value for kw `location`: received $location.
           Valid options are `:onsurface`, `:inside` and `:outside`."
    σ = location === :onsurface ? -0.5 : location === :inside ? 0 : location === :outside ? -1 : error(msg)
    Sop  = SingleLayerOperator(pde,X,Y)
    Dop  = DoubleLayerOperator(pde,X,Y)
    # convert to a possibly more efficient format
    @timeit_debug "assemble dense part" begin
        S = compress(Sop)
        D = compress(Dop)
    end
    @timeit_debug "compute near interaction list" begin
        dict_near = near_interaction_list(dofs(X),Y;atol=0)
    end
    # precompute dim quantities
    @timeit_debug "auxiliary dim quantities" begin
        basis,γ₁_basis = _basis_dim(Sop)  # list of functions γ₀ and γ₁ for each source
        γ₀B,γ₁B,R      = _auxiliary_quantities_dim(Sop,S,D,basis,γ₁_basis,σ)
    end
    # compute corrections
    @timeit_debug "compute dim correction" begin
        if pde isa Maxwell
            δS = _singular_weights_dim_maxwell(Sop,γ₀B,γ₁B,R,dict_near)
            δD = _singular_weights_dim_maxwell(Dop,γ₀B,γ₁B,R,dict_near)
        else
            δS = _singular_weights_dim(Sop,γ₀B,γ₁B,R,dict_near)
            δD = _singular_weights_dim(Dop,γ₀B,γ₁B,R,dict_near)
        end
    end
    # convert to DiscreteOp
    S_discrete = DiscreteOp(S) + DiscreteOp(δS)
    D_discrete = DiscreteOp(D) + DiscreteOp(δD)
    return S_discrete,D_discrete
end
singlelayer_dim(args...;kwargs...) = single_doublelayer_dim(args...;kwargs...)[1]
doublelayer_dim(args...;kwargs...) = single_doublelayer_dim(args...;kwargs...)[2]

function adjointdoublelayer_hypersingular_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
    msg = "unrecognized value for kw `location`: received $location.
    Valid options are `:onsurface`, `:inside` and `:outside`."
    σ = location === :onsurface ? -0.5 : location === :inside ? 0 : location === :outside ? -1 : error(msg)
    Kop  = AdjointDoubleLayerOperator(pde,X,Y)
    Hop  = HyperSingularOperator(pde,X,Y)
    # convert to a possibly more efficient compress
    K = compress(Kop)
    H = compress(Hop)
    dict_near = near_interaction_list(dofs(X),Y;atol=0)
    # precompute dim quantities
    basis,γ₁_basis = _basis_dim(Kop)
    γ₀B,γ₁B,R      = _auxiliary_quantities_dim(Kop,K,H,basis,γ₁_basis,σ)
    # compute corrections
    δK        = _singular_weights_dim(Kop,γ₀B,γ₁B,R,dict_near)
    δH        = _singular_weights_dim(Hop,γ₀B,γ₁B,R,dict_near)
    # convert to DiscreteOp
    K_discrete = DiscreteOp(K) + DiscreteOp(δK)
    H_discrete = DiscreteOp(H) + DiscreteOp(δH)
    return K_discrete,H_discrete
end
adjointdoublelayer_dim(args...;kwargs...)  = adjointdoublelayer_hypersingular_dim(args...;kwargs...)[1]
hypersingular_dim(args...;kwargs...)       = adjointdoublelayer_hypersingular_dim(args...;kwargs...)[2]

function singular_weights_dim(iop::IntegralOperator,compress=Matrix)
    X,Y,op = iop.X, iop.Y, iop.kernel.op
    σ = X == Y ? -0.5 : 0.0
    #
    dict_near = near_interaction_list(dofs(X),Y;atol=0)
    basis,γ₁_basis = _basis_dim(iop)
    Op1, Op2       = _auxiliary_operators_dim(iop,compress)
    γ₀B,γ₁B,R      = _auxiliary_quantities_dim(iop,Op1,Op2,basis,γ₁_basis,σ)
    Sp = _singular_weights_dim(iop,γ₀B,γ₁B,R,dict_near)
    return Sp # a sparse matrix
end

function _auxiliary_quantities_dim(iop,Op0,Op1,basis,γ₁_basis,σ)
    T      = eltype(iop)
    X,Y    = target_surface(iop), source_surface(iop)
    num_basis = length(basis)
    # compute matrix of basis evaluated on Y
    ynodes   = dofs(Y)
    γ₀B      = Matrix{T}(undef,length(ynodes),num_basis)
    γ₁B      = Matrix{T}(undef,length(ynodes),num_basis)
    for k in 1:num_basis
        for i in 1:length(ynodes)
            γ₀B[i,k] = basis[k](ynodes[i])
            γ₁B[i,k] = γ₁_basis[k](ynodes[i])
        end
    end
    # integrate the basis over Y
    xnodes      = dofs(X)
    num_targets = length(xnodes)
    num_sources = length(ynodes)
    R = Matrix{T}(undef,num_targets,num_basis)
    @timeit_debug "integrate dim basis" begin
        @threads for k in 1:num_basis
            # integrate basis over the surface
            # FIXME: use inplace multiplication through `mul!` instead of the line
            # below. The issue at the moment is that
            R[:,k] = Op0*γ₁B[:,k] - Op1*γ₀B[:,k]
            # @views mul!(R[:,k],Op0,γ₁B[:,k])
            # @views mul!(R[:,k:k],Op1,γ₀B[:,k:k],-1,1)
            for i in 1:num_targets
                # analytic correction for on-surface evaluation of Greens identity
                if kernel_type(iop) isa Union{SingleLayer,DoubleLayer}
                    R[i,k] += σ*basis[k](xnodes[i])
                elseif kernel_type(iop) isa Union{AdjointDoubleLayer,HyperSingular}
                    R[i,k] += σ*γ₁_basis[k](xnodes[i])
                end
            end
        end
    end
    return γ₀B, γ₁B, R
end

function _auxiliary_operators_dim(iop,compress)
    X,Y,op = iop.X, iop.Y, iop.kernel.op
    T = eltype(iop)
    # construct integral operators required for correction
    if kernel_type(iop) isa Union{SingleLayer,DoubleLayer}
        Op1 = IntegralOperator{T}(SingleLayerKernel(op),X,Y) |> compress
        Op2 = IntegralOperator{T}(DoubleLayerKernel(op),X,Y) |> compress
    elseif kernel_type(iop) isa Union{AdjointDoubleLayer,HyperSingular}
        Op1 = IntegralOperator{T}(AdjointDoubleLayerKernel(op),X,Y) |> compress
        Op2 = IntegralOperator{T}(HyperSingularKernel(op),X,Y) |> compress
    end
    return Op1,Op2
end

function _basis_dim(iop)
    op = pde(kernel(iop))
    xs = _source_gen(iop) # list of Lebedev sources
    basis     = [(source) -> SingleLayerKernel(op)(x,source) for x in xs]
    γ₁_basis  = [(source) -> transpose(DoubleLayerKernel(op)(x,source)) for x in xs]
    return basis,γ₁_basis
end

function _singular_weights_dim(iop::IntegralOperator,γ₀B,γ₁B,R,dict_near)
    X,Y = target_surface(iop), source_surface(iop)
    T   = eltype(iop)
    num_basis = size(γ₀B,2)
    a,b = combined_field_coefficients(iop)
    # we now have the residue R. For the correction we need the coefficients.
    Is = Int[]
    Js = Int[]
    Vs = T[]
    for (E,list_near) in dict_near
        el2qnodes = elt2dof(Y,E)
        num_qnodes, num_els   = size(el2qnodes)
        M                     = Matrix{T}(undef,2*num_qnodes,num_basis)
        @assert length(list_near) == num_els
        for n in 1:num_els
            j_glob                = @view el2qnodes[:,n]
            M[1:num_qnodes,:]     = @view γ₀B[j_glob,:]
            M[num_qnodes+1:end,:] = @view γ₁B[j_glob,:]
            # distinguish scalar and vectorial case
            if T <: Number
                F                     = qr(M)
            elseif T <: SMatrix
                M_mat = blockmatrix_to_matrix(M)
                F                     = qr!(M_mat)
            else
                error("unknown element type T=$T")
            end
            for (i,_) in list_near[n]
                if T <: Number
                    tmp = ((R[i:i,:])/F.R)*adjoint(F.Q)
                elseif T <: SMatrix
                    tmp_scalar  = (blockmatrix_to_matrix(R[i:i,:])/F.R)*adjoint(F.Q)
                    tmp  = matrix_to_blockmatrix(tmp_scalar,T)
                else
                    error("unknown element type T=$T")
                end
                w    = axpby!(a,view(tmp,1:num_qnodes),b,view(tmp,(num_qnodes+1):(2*num_qnodes)))
                append!(Is,fill(i,num_qnodes))
                append!(Js,j_glob)
                append!(Vs,w)
            end
        end
    end
    Sp = sparse(Is,Js,Vs,size(iop)...)
    return Sp
end

function _singular_weights_dim_maxwell(iop::IntegralOperator,γ₀B,γ₁B,R,dict_near)
    # initialize vectors for the sparse matrix, then dispatch to type-stable
    # method for each element type
    T   = eltype(iop)
    Is = Int[]
    Js = Int[]
    Vs = T[]
    for (E,list_near) in dict_near
        _singular_weights_dim_maxwell!(Is,Js,Vs,iop,γ₀B,γ₁B,R,E,list_near)
    end
    Sp = sparse(Is,Js,Vs,size(iop)...)
    return Sp
end

@noinline function _singular_weights_dim_maxwell!(Is,Js,Vs,iop,γ₀B,γ₁B,R,E,list_near)
    X,Y    = target_surface(iop), source_surface(iop)
    qnodes = dofs(Y)
    T   = eltype(iop)
    a,b = combined_field_coefficients(iop)
    el2qnodes = elt2dof(Y,E)
    num_qnodes, num_els   = size(el2qnodes)
    num_basis = size(γ₀B,2)
    T2 = SMatrix{2,3,ComplexF64,6}
    T3 = SMatrix{3,2,ComplexF64,6}
    M      = Matrix{T2}(undef,2*num_qnodes,num_basis)
    M_mat  = Matrix{ComplexF64}(undef,2*2*num_qnodes,3*num_basis)
    @assert length(list_near) == num_els
    for n in 1:num_els
        j_glob                = @view el2qnodes[:,n]
        for p in 1:num_basis
            for k in 1:num_qnodes
                # J,_ = jacobian(qnodes[j_glob[k]]) |> qr
                J = jacobian(qnodes[j_glob[k]])
                M[k,p]            = transpose(J)*γ₀B[j_glob[k],p]
                M[num_qnodes+k,p] = transpose(J)*γ₁B[j_glob[k],p]
            end
        end
        # M_mat                 = blockmatrix_to_matrix(M)
        blockmatrix_to_matrix!(M_mat,M)
        F                     = qr!(M_mat)
        for (i,_) in list_near[n]
            tmp_scalar  = (blockmatrix_to_matrix(R[i:i,:])/F.R)*adjoint(F.Q)
            tmp         = matrix_to_blockmatrix(tmp_scalar,T3)
            tmp         = axpby!(a,view(tmp,1:num_qnodes),b,view(tmp,(num_qnodes+1):(2*num_qnodes)))
            w           = Vector{T}(undef,num_qnodes)
            for k in 1:num_qnodes
                # J,_ = jacobian(qnodes[j_glob[k]]) |> qr
                J    = jacobian(qnodes[j_glob[k]])
                w[k] = tmp[k]*transpose(J)
            end
            append!(Is,fill(i,num_qnodes))
            append!(Js,j_glob)
            append!(Vs,w)
        end
    end
    return Is,Js,Vs
end

function _source_gen(iop::IntegralOperator,kfactor=5)
    Y      =  source_surface(iop)
    nquad  = 0
    for (E,tags) in elt2dof(Y)
        nquad = max(nquad,size(tags,1))
    end
    nbasis = 3*nquad
    # construct source basis
    return _source_gen(iop,nbasis;kfactor)
end

function _source_gen(iop,nsources;kfactor)
    N      = ambient_dimension(iop)
    Y      = source_surface(iop)
    pts    = dofs(Y)
    # create a bounding box
    bbox   = HyperRectangle(pts)
    xc     = center(bbox)
    d      = diameter(bbox)
    if N == 2
        xs = _circle_sources(;nsources,center=xc,radius=kfactor*d/2)
    elseif N == 3
        xs = _sphere_sources_lebedev(;nsources,center=xc,radius=kfactor*d/2)
    else
        error("dimension must be 2 or 3. Got $N")
    end
    return xs
end

function _sphere_sources_lebedev(;nsources, radius=10, center=SVector(0.,0.,0.))
    lpts = lebedev_points(nsources)
    Xs = SVector{3,Float64}[]
    for pt in lpts
        push!(Xs,radius*pt .+ center)
    end
    return Xs
end

function _circle_sources(;nsources, radius=10, center=SVector(0.,0.))
    par   = (s) -> center .+ radius .* SVector(cospi(2 * s[1]), sinpi(2 * s[1]))
    x     = [i/(nsources) for i in 0:nsources-1]
    Xs    = SVector{2,Float64}[]
    for pt in x
        push!(Xs,par(pt))
    end
    return Xs
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module BoundedPriorityQueues

export BoundedPriorityQueue, enqueue!

using DataStructures

mutable struct BoundedPriorityQueue{K,V}
    pq::PriorityQueue{K,V}
    N::Int64

    function BoundedPriorityQueue{K,V}(N::Int64, o::Base.Order.Ordering=Base.Order.Forward) where {K,V}
        #higher is kept
        new{K,V}(PriorityQueue{K,V}(o), N)
    end
end

function DataStructures.enqueue!(q::BoundedPriorityQueue{K,V}, k::K, v::V;
    make_copy::Bool=false) where {K,V}
    haskey(q.pq, k) && return -1 #keys must be unique, return -1 if collision
    if make_copy
        k = deepcopy(k)
    end
    n = length(q.pq)
    enqueue!(q.pq, k, v)
    n_add = n - length(q.pq) #number of items added
    while length(q.pq) > q.N
        dequeue!(q.pq)
    end
    n_add
end

Base.length(q::BoundedPriorityQueue) = length(q.pq)
Base.isempty(q::BoundedPriorityQueue) = isempty(q.pq)
Base.haskey(q::BoundedPriorityQueue) = haskey(q.pq)
Base.keys(q::BoundedPriorityQueue) = keys(q.pq)
Base.values(q::BoundedPriorityQueue) = values(q.pq)
function Base.empty!(q::BoundedPriorityQueue)
    while !isempty(q.pq)
        dequeue!(q.pq)
    end
end

"""
Ordered iterator
"""
function Base.iterate(q::BoundedPriorityQueue)
    pairs = collect(q.pq)
    sort!(pairs, by=x->x[2], rev=(q.pq.o==Base.Order.ForwardOrdering()))
    result = iterate(pairs)
    if result == nothing
        return nothing
    else
        item, state = result
        return item, (pairs=pairs, state=state)
    end
end
function Base.iterate(q::BoundedPriorityQueue, x)
    result = iterate(x.pairs, x.state)
    if result == nothing
        return nothing
    else
        item, state = result
        return item, (pairs=x.pairs, state=state)
    end
end

end #module


text:
@with_kw mutable struct AStarGrid2DEnvironment <: MAPFEnvironment
    dimx::Int64
    dimy::Int64
    obstacles::Set{Grid2DLocation}
    goal::Grid2DLocation
    state_graph::SimpleVListGraph{Grid2DLocation}    = SimpleVListGraph{Grid2DLocation}()
    state_to_idx::Dict{Grid2DLocation,Int64}         = Dict{Grid2DLocation,Int64}()
    goal_idx::Int64                                  = 0
end

function admissible_heuristic_a_star_eps(env::AStarGrid2DEnvironment, s::Grid2DLocation)
    return abs(s.x - env.goal.x) + abs(s.y - env.goal.y)
end

function focal_state_heuristic_a_star_eps(env::AStarGrid2DEnvironment, s::Grid2DLocation, gscore::Int64)
    return gscore
end

function focal_transition_heuristic_a_star_eps(env::AStarGrid2DEnvironment, s1::Grid2DLocation, s2::Grid2DLocation,
                                    gscore1::Int64, gscore2::Int64)
    return gscore2 - gscore1
end

function state_valid(env::AStarGrid2DEnvironment, s::Grid2DLocation)
    return s.x >=0 && s.x < env.dimx && s.y >=0 && s.y < env.dimy &&
            ~(s in env.obstacles)
end

function get_env_state_idx!(env::AStarGrid2DEnvironment, s::Grid2DLocation)

    idx = get(env.state_to_idx, s, 0)

    if idx == 0
        add_vertex!(env.state_graph, s)
        nbr_idx = num_vertices(env.state_graph)
        env.state_to_idx[s] = nbr_idx
        idx = nbr_idx
    end

    return idx
end

function get_mapf_state_from_idx(env::AStarGrid2DEnvironment, idx::Int64)
    return env.state_graph.vertices[idx]
end

function get_mapf_action(env::AStarGrid2DEnvironment, source::Int64, target::Int64)

    state1 = env.state_graph.vertices[source]
    state2 = env.state_graph.vertices[target]

    if state2.x > state1.x
        return Grid2DAction(Right::Action)
    elseif state2.x < state1.x
        return Grid2DAction(Left::Action)
    elseif state2.y > state1.y
        return Grid2DAction(Up::Action)
    else
        return Grid2DAction(Down::Action)
    end
    return nothing
end

function low_level_search!(env::AStarGrid2DEnvironment, s::Grid2DLocation, weight::Float64)

    source_idx = get_env_state_idx!(env, s)

    # Set all edge weights to one
    edge_wt_fn(u, v) = 1

    # Set the heuristics
    admissible_heuristic(s) = admissible_heuristic_a_star_eps(env, s)
    focal_state_heuristic(s, gs) = focal_state_heuristic_a_star_eps(env, s, gs)
    focal_transition_heuristic(s1, s2, gs1, gs2) = focal_transition_heuristic_a_star_eps(env, s1, s2, gs1, gs2)

    # Run the search
    a_star_epsilon_states = a_star_light_epsilon_shortest_path_implicit!(env.state_graph,
                                                                        edge_wt_fn, source_idx,
                                                                        AStarEpsilonGrid2DVisitor(env),
                                                                        weight,
                                                                        admissible_heuristic,
                                                                        focal_state_heuristic,
                                                                        focal_transition_heuristic,
                                                                        Int64)

    # Extract the solution as a vector of locations and dist tuples
    if ~(haskey(a_star_epsilon_states.dists, env.goal_idx))
        return nothing, 0
    end

    cost = a_star_epsilon_states.dists[env.goal_idx]

    goal_state = env.state_graph.vertices[env.goal_idx]
    solution = [goal_state]

    curr_idx = env.goal_idx
    while curr_idx != source_idx

        prev_idx = a_star_epsilon_states.parent_indices[curr_idx]
        pushfirst!(solution, env.state_graph.vertices[prev_idx])
        curr_idx = prev_idx

    end

    return solution, cost
end

mutable struct AStarEpsilonGrid2DVisitor <: AbstractDijkstraVisitor
    env::AStarGrid2DEnvironment
end

function Graphs.include_vertex!(vis::AStarEpsilonGrid2DVisitor, u::Grid2DLocation, v::Grid2DLocation,
                                d::N, nbrs::Vector{Int64}) where {N <: Number}

    env = vis.env
    if v == env.goal
        env.goal_idx = get_env_state_idx!(env, v)
        return false
    end

    empty!(nbrs)

    temp_locs = [Grid2DLocation((x = v.x-1, y = v.y)), Grid2DLocation((x = v.x+1, y = v.y)),
                 Grid2DLocation((x = v.x, y = v.y-1)), Grid2DLocation((x = v.x, y = v.y+1))]

    for temp_loc in temp_locs
        if state_valid(env, temp_loc)
            nbr_idx = get_env_state_idx!(env, temp_loc)
            push!(nbrs, nbr_idx)
        end
    end

    return true
end


text:
"""
# module Elements

# Examples

```jldoctest
julia>
```
"""
module Elements

import PeriodicTable

export AtomicSpecies,
    get_element,
    get_atomic_species

struct AtomicSpecies
    symbol::String
    mass::Float64
    pseudopotential::String
end

AtomicSpecies(symbol::String, mass::Float64) = pseudopotential::String -> AtomicSpecies(symbol, mass, pseudopotential)
AtomicSpecies(symbol::String) = mass::Float64 -> AtomicSpecies(symbol, mass)

get_element(index::Symbol) = PeriodicTable.elements.bysymbol[(Symbol ∘ uppercasefirst ∘ String)(index)]
get_element(index::Int) = PeriodicTable.elements.bynumber[index]
get_element(index::String) = PeriodicTable.elements.byname[lowercase(index)]

function get_atomic_species(index::Union{Symbol, Int, String})
    element = get_element(index)  # Separate kernel functions (aka, function barriers)
    AtomicSpecies(element.symbol, element.atomic_mass.val)
end

end


text:
export doubledescription
# Naive implementation of the double description Algorithm
# See JuliaPolyhedra/ConvexHull.jl for an efficient implementation

polytypefor(::Type{T}) where {T <: Integer} = Rational{BigInt}
polytypefor(::Type{Float32}) = Float64
polytypefor(::Type{T}) where {T} = T

polyvectortype(a) = a
# TODO sparse halfspaces does not mean sparse points
polyvectortype(::Type{<:SparseVector{T}}) where T = Vector{T}

dualtype(RepT::Type{<:Representation}) = dualtype(RepT, polyvectortype(vectortype(RepT)))
function dualfullspace(rep::Representation, d::FullDim, ::Type{T}) where T
    dualfullspace(rep, d, T, polyvectortype(similar_type(vectortype(rep), d, T)))
end
function dualfullspace(rep::Representation{T}) where T
    dualfullspace(rep, FullDim(rep), polytypefor(T))
end

"""
    doubledescription(h::HRepresentation)

Computes the V-representation of the polyhedron represented by `h` using the Double-Description algorithm [1, 2].

    doubledescription(V::VRepresentation)

Computes the H-representation of the polyhedron represented by `v` using the Double-Description algorithm [1, 2].

[1] Motzkin, T. S., Raiffa, H., Thompson, G. L. and Thrall, R. M.
The double description method
*Contribution to the Theory of Games*, *Princeton University Press*, **1953**

[2] Fukuda, K. and Prodon, A.
Double description method revisited
*Combinatorics and computer science*, *Springer*, **1996**, 91-111
"""
function doubledescription end

function print_v_summary(v::VRep)
    print("$(npoints(v)) points, $(nrays(v)) rays and $(nlines(v)) lines")
end

function intersect_and_remove_redundancy(v, hs, h; verbose=0)
    if eltype(hs) <: HalfSpace
        str = "halfspace"
    else
        str = "hyperplane"
    end
    for (i, hel) in enumerate(hs)
        if verbose >= 1
            println("Intersecting $str $i/$(length(hs))")
        end
        v_int = v ∩ hel
        if verbose >= 3
            print("Removing duplicates: ")
            print_v_summary(v_int)
            println(".")
        end
        # removeduplicates is cheaper than removevredundancy since the latter
        # needs to go through all the hrep element
        v_dup = removeduplicates(v_int)
        if verbose >= 3
            print("Removing redundancy: ")
            print_v_summary(v_dup)
            println(".")
        end
        v = removevredundancy(v_dup, h)
        if verbose >= 2
            print("After intersection:  ")
            print_v_summary(v)
            println(".")
        end
    end
    return v
end

function doubledescription(h::HRepresentation; kws...)
    # The redundancy of V-elements are measured using
    # the number of hyperplane they are in. If there are
    # redundant hyperplanes, it does not matter since no point
    # will be inside them but if there are duplicates it is a problem

    # FIXME Note that removevredundancy, uses `h` which contains all hyperplanes and halfspaces
    # already taken into account but also all the other ones. We should check that this
    # is the right idea.
    h = removeduplicates(h)
    v = dualfullspace(h)
    checkvconsistency(v)
    v = intersect_and_remove_redundancy(v, hyperplanes(h), h; kws...)
    v = intersect_and_remove_redundancy(v, halfspaces(h), h; kws...)
    return v
end

function doubledescription(v::VRepresentation{T}; kws...) where {T}
    checkvconsistency(v)
    lv = convert(LiftedVRepresentation{T, Matrix{T}}, v)
    R = -lv.R
    vl = doubledescription(MixedMatHRep{T}(R, zeros(T, size(R, 1)), lv.linset); kws...)
    LiftedHRepresentation{T}(vl.R, vl.Rlinset)
end


text:
@forwardfun convert_to_dataset
@forwardfun dict_to_dataset

"""
    dataset_to_dict(ds::PyObject)

Convert an `xarray.Dataset` to a dictionary of `Array`s. The function is the
returns a `Dict{String,Array}` and a `NamedTuple` with keyword arguments to
`dict_to_dataset`.
"""
function dataset_to_dict(ds::PyObject)
    ds_dict = ds.to_dict()
    data_vars = ds_dict["data_vars"]
    attrs = ds_dict["attrs"]

    coords = ds_dict["coords"]
    delete!(coords, "chain")
    delete!(coords, "draw")
    coords = Dict(k => v["data"] for (k, v) in coords)

    data = Dict{String,Array}()
    dims = Dict{String,Vector{String}}()
    for (k, v) in data_vars
        data[k] = v["data"]
        dim = v["dims"][3:end]
        if !isempty(dim)
            dims[k] = [dim...]
        end
    end

    return data, (attrs = attrs, coords = coords, dims = dims)
end


text:
using CaimCore
using CaimCore.ImageStack
using Colors, FixedPointNumbers

const STACKS_DIR = joinpath(@__DIR__, "data", "stacks")
const VIDEOS_DIR = joinpath(@__DIR__, "data", "videos")

@testset "Series" begin
    let stack = load(Frames, joinpath(STACKS_DIR, "png"))
        box = Box([1,1], [1,1])
        ts = series(stack, box)
        @test ts ≈ [RGB{N0f8}(1.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 1.0),
                    RGB{N0f8}(0.0, 1.0, 0.0),
                    RGB{N0f8}(1.0, 0.0, 0.0)]

        box = Box([1,1], [10,10])
        ts = series(stack, box)
        @test ts ≈ [RGB{N0f8}(1.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 1.0),
                    RGB{N0f8}(0.0, 1.0, 0.0),
                    RGB{N0f8}(1.0, 0.0, 0.0)]

		box = Box([5,5], [15,15])
		ts = series(stack, box)
		@test ts ≈ [RGB{N0f8}(6/11, sqrt(30)/11, 5/11),
		            RGB{N0f8}(sqrt(30)/11, 5/11, sqrt(30)/11),
					RGB{N0f8}(5/11, sqrt(30)/11, 6/11),
					RGB{N0f8}(sqrt(30)/11, 6/11, sqrt(30)/11),
					RGB{N0f8}(6/11, sqrt(30)/11, 5/11)]

		ff = FreeForm([[10,1], [10,10], [1,10]])
		ts = series(stack, ff)
		@test ts ≈ [RGB{N0f8}(1.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 0.0),
                    RGB{N0f8}(0.0, 0.0, 1.0),
                    RGB{N0f8}(0.0, 1.0, 0.0),
                    RGB{N0f8}(1.0, 0.0, 0.0)]
    end
end


text:
include(joinpath(dirname(pathof(Scallops)), "plotting.jl"))
using DataDeps, FileIO, Colors, CoordinateTransformations, ImageMorphology, Images, Interpolations
register(DataDep("scalloPupil", "a gif of scallop pupils", "https://vision-group-temporary.s3.amazonaws.com//scallop_pupil_Pm_green1_alt.gif", "4c8a194dd8a5df37833016d48a2f0f17cca321a8a6618984ddbc0a716b692063"))

folder = @datadep_str "scalloPupil"
file = first(readdir(folder))
gif = Gray.(load(joinpath(folder, file)))
bw = gif .< 0.1
bw = opening(bw)
lb = label_components(bw)
function fun(bb)
    XY = [zeros(2) for _ in bb]
    for (i,b) in enumerate(bb[2:2])
        xy, wh = b
        x,y = xy
        w,h = wh
        XY[i][1] = 2x + w
        XY[i][2] = 2y + h
    end
XY
end
wgif = [warp(gif[:,:,1], Translation(0,0))]
i = 1
bb = fun(component_boxes(lb[:,:,i]))
maxΔ = 0.0
for i in 2:size(lb, 3)
    global bb, maxΔ
    _bb = fun(component_boxes(lb[:,:,i]))
    Δ = -mean(bb .- _bb)
    maxΔ = max(maxΔ, abs.(Δ)...)
    tfm = Translation(Δ)
    push!(wgif, warp(gif[:,:,i], tfm))
end
w = cat(paddedviews(0, wgif...)..., dims = 3)
h = ceil(Int, maxΔ)
w = w[h+1:end-h, h+1:end-h, :]

bw = gif .> 0.1
bw = opening(bw)
lb = label_components(bw)

r = Vector{Float64}(undef, size(lb, 3))
for i in 1:size(lb,3)
    a = component_lengths(lb[:,:,i])[3]
    r[i] = sqrt.(a/π)
end
r .-= minimum(r)
r ./= maximum(r)
i = sortperm(r)
r = r[i]
sortedw = w[:,:,i]
itp = interpolate((1:size(w,1), 1:size(w,2), r), sortedw, Gridded(Linear()))
rl = range(0, 1, length=2size(w,3))
w2 = Array{Gray{Normed{UInt8,8}},3}(undef, size(w,1), size(w, 2), length(rl))
for i in 1:size(w,1), j in 1:size(w,2), (k,r) in enumerate(rl)
    w2[i,j,k] = itp(i,j,r)
end

min_aperture = 200
max_aperture = 400
aperture_n = size(w2,3)
min_morph = 0.8
max_morph = 1.48
morph_n = 100
min_distance = 0.1
max_distance = 1000
distance_n = 100
min_θ = -38
max_θ = 38
θ_n = max_θ - min_θ + 1
limits = FRect(-350,-210,700,420)
nvertices = 21
photoreceptor_radius = 5
absorption_coefficient = (dis_ret = 0.00667, pro_ret = 0.00667)
ranges = getranges(min_aperture, max_aperture, aperture_n, min_morph, max_morph, morph_n, min_distance, max_distance, distance_n, min_θ, max_θ, θ_n)
sliders, observables = getobservables(ranges...)
all2 = lift(createobj, observables.morphz, photoreceptor_radius, absorption_coefficient)
sc = Scene(scale_plot=false, limits=limits)
vertices2 = lift(all2) do all2
    populateθ0!(all2.layers)
    todo = [k for k in keys(all2.layers) if k ∉ (:water, :mirror)]
    map(todo) do k
        getvertices(all2.layers, k, nvertices)
    end
end
colors2 = lift(all2) do all2
    todo = [k for k in keys(all2.layers) if k ∉ (:water, :mirror)]
    m, M = extrema(all2.layers[k].tissue.ri for k in todo)
    Δ = M - m
    m -= 0.1Δ
    M += 0.1Δ
    if Δ == 0
        m = copy(M)
        M += 1
    end
    map(todo) do k
        layers2color(all2.layers, k, m, M)
    end
end
poly!(sc, vertices2, color = colors2)
rays = lift(observables.distance, observables.aperture, all2, observables.θ) do distance, aperture, all2, θ
    l = Light(mm2μm(distance), aperture, all2.rc, deg2rad(θ), all2.nodalz)
    getrays(all2.ous, l)
end
segments = lift(rays) do rays
    [ray[i].coord => ray[i+1].coord for ray in rays for i in 1:length(ray) - 1]
end
colors2 = lift(rays) do rays
    [mean((ray[i].color, ray[i+1].color)) for ray in rays for i in 1:length(ray) - 1]
end
linesegments!(sc, segments, color=colors2)
axis = sc[Makie.Axis];
axis[:names, :axisnames] = ("X (μm)", "Y (μm)");
observables.morphz[] = 1.1
t = Node(1)
aspect_ratio = size(w2, 1)/size(w2, 2)
image!(sc, range(-400,400,length = size(w2,2)), range(250,250+800*aspect_ratio, length = size(w2,1)), lift(i -> w2[end:-1:1,:,i]', t), scale_plot=false, show_axis = false)
txtcolor = :black
rctcolor = :white
h = 15.0
y = 190.0
poly!(sc, FRect(-h,y-h,2h,2h), color = rctcolor)
text!(sc, "C", position = (0.0, y), align = (:center,  :center), textsize = 2h, color = txtcolor)
y = 60.0
poly!(sc, FRect(-h,y-h,2h,2h), color = rctcolor)
text!(sc, "L", position = (0.0, y), align = (:center,  :center), textsize = 2h, color = txtcolor)
y = -62.5
poly!(sc, FRect(-h,y-h,2h,2h), color = rctcolor)
text!(sc, "D", position = (0.0, y), align = (:center,  :center), textsize = 2h, color = txtcolor)
y = -150.0
poly!(sc, FRect(-h,y-h,2h,2h), color = rctcolor)
text!(sc, "P", position = (0.0, y), align = (:center,  :center), textsize = 2h, color = txtcolor)
y = -220.0
poly!(sc, FRect(-h,y-h,2h,2h), color = rctcolor)
text!(sc, "M", position = (0.0, y), align = (:center,  :center), textsize = 2h, color = txtcolor)
apertures = ranges.aperture#[end:-1:1]
record(sc, "pupil.mp4", 1:aperture_n) do i
    if i == 1
        observables.distance[] = maximum(ranges.distance)
    end
    observables.aperture[] = apertures[i]
    t[] = i
end


text:
__precompile__(true)

module Backtester

using Distributed, Statistics, Printf, Dates, DelimitedFiles, Random
using Glob, DataStructures, Parameters, ProgressMeter
using PyCall, StandardMarketData, MLSuiteBase
using Iconv, HDF5Utils, PandasLite, PyCallUtils
using StatsBase: corspearman

export Strategy, backtest, combine
export eglob, rglob
export @NT, @redirect, @indir, @gc, @roll

include("util.jl")
include("strategy.jl")
include("backtest.jl")
include("output.jl")

end


text:
module VTKDataIO

using VTKDataTypes
using Printf
using WriteVTK
using PyCall
using LightXML
using Base.Iterators

@pyimport vtk.util.numpy_support as vtkns
@pyimport vtk as vtk
@pyimport numpy as np

include("vtkreaders.jl")
include("vtkwriters.jl")
include("vtpwriter.jl")
include("stlIO.jl")
include("extra_readers.jl")
include("PyVTK.jl")
include("visualize.jl")

export read_vtk, write_vtk, read_stl, write_stl, valid_to_write, PyVTK,
    _VTKDataTypes, visualize, write_x3d, write_ply, read_obj, read_ply, visualize_3ds

end


text:
const AWS_LAMBDA_RUNTIME_API = ENV["AWS_LAMBDA_RUNTIME_API"]
const LAMBDA_TASK_ROOT = ENV["LAMBDA_TASK_ROOT"]
const HANDLER = ENV["_HANDLER"]  # from Dockerfile's CMD command
const RUNTIME_URL = "http://$AWS_LAMBDA_RUNTIME_API/2018-06-01"

# Using curl as HTTP client.
# (Don't want to introduce HTTP.jl dependency for faster initialization)

function http_get(url)
    cmd = `curl --silent --show-error --location --include --request GET $url`
    response = read(cmd, String)
    headers_str, body = split(response, "\r\n\r\n")
    headers = [s[1] => s[2] for s in split.(split(headers_str, "\r\n"), r": *") if length(s) == 2]
    return headers, body
end

function http_post(url, headers, body)
    header_lines = ["$k: $v" for (k,v) in headers]
    header_options = collect(Iterators.flatten(("-H", v) for v in header_lines))
    cmd = `curl --silent --show-error $header_options --request POST $url --data $body`
    return read(cmd, String)
end

# Convenient function to post errors to the lamdba execution environment
function post_error(path::AbstractString, ex::Exception)
    headers = [
        "Content-type" => "application/json",
        "Lambda-Runtime-Function-Error-Type" => "Unhandled",
    ]
    local url
    try
        url = "$RUNTIME_URL/$path"
        http_post(url, headers, """{"errorType": "$(typeof(ex))", "errorMessage": "$(ex)"}""")
        @info "Notified lambda runtime about the error" url ex
    catch failure
        @error "Unable to notify lambda runtime about the error" url ex failure
    end
end

# Initialize lambda function by loading user module
try
    @info "Loading user module"
    global mod, func = Symbol.(split(HANDLER, "."))

    @eval using $(mod)

    @info "Completed loading user module" mod func
catch ex
    @error "Initializtion error" ex
    post_error("runtime/init/error", ex)
    @info "Shutting down container"
    exit(1)
end

@info "Resolve handler function"
const my_module = Base.eval(Main, mod)
const my_handler = Base.eval(my_module, func)

# An infinite loop to process events
@info "Start processing requests"
while true
    local state, request_id
    try
        state = :start
        request_headers, request_body = http_get("$RUNTIME_URL/runtime/invocation/next")
        request_id_idx = findfirst(x -> x[1] == "Lambda-Runtime-Aws-Request-Id", request_headers)
        request_id = request_headers[request_id_idx][2]
        @info "Received event" request_id
        state = :received

        response  = my_handler(request_body, request_headers)
        @info "Got response from handler" response
        state = :handled

        http_post("$RUNTIME_URL/runtime/invocation/$request_id/response", [], response)
        @info "notified lambda runtime"
        state = :finished
    catch ex
        if state == :start
            @error "Unable to receive request" ex
        elseif state == :received
            @error "Unable to handle event" ex request_id
            post_error("runtime/invocation/$request_id/error", ex)
        elseif state == :handled
            @error "Unable to notify lambda runtime (invocation response)" ex request_id
        else
            @error "Unknown state" ex state
        end
    end
end

@info "All done"

max_stars_repo_name: ramboldio/aws-lambda-container-julia

max_issues_repo_name: ramboldio/aws-lambda-container-julia

max_forks_repo_name: ramboldio/aws-lambda-container-julia


text:
"""
  Package: FishEBM
  File: simulate.jl

  Description: Brings together all of the functions necessary for a life cycle
    simulation

  Returns: Vector(EnviroAgent)

  Contributors: Justin Angevaare, Devin Rose

  Last update: October 2016
"""

function simulate(carrying_capacity::Vector{Float64}, effort::Vector{Float64}, initStock::Vector{Int64},
  stock_age::Vector{Int64}, e_a::EnvironmentAssumptions, adult_a::AdultAssumptions,
  age_a::AgentAssumptions; progress::Bool=true, plotPopDensity::Bool=false,
  plotPopDistribution::Bool=false, limit::Int64=1000000, simDescription::String="")

  # preconditions
  @assert(all(carrying_capacity .> 0.), "There is at least one negative carrying capacity")
  @assert(plotPopDensity == false || plotPopDistribution == false, "Only one plot can be run during
    a simulation")
  @assert(length(initStock) == length(stock_age), "Unmatching vector length for initializing the
    population, stock length = $(length(initStock)) & $(length(stock_age))")
  @assert(length(effort) == length(carrying_capacity) || (length(effort) == 1 && effort[1] == 0.0),
    "Effort must be the same length as the carrying capacity, or a length of one with a value of zero (effort[1] = 0) for no harvest")

  # initialize agent db and hash the enviro
  years = length(carrying_capacity)
  a_db = AgentDB(e_a); hashEnvironment!(a_db, e_a);
  if plotPopDensity
    popDensity = initPopulationDensity(e_a)
  end

  #initialize the stock with incoming paramaters
  for i = 1:length(initStock)
    injectAgents!(a_db, e_a.spawningHash, initStock[length(initStock)+1-i], stock_age[length(initStock)+1-i])
  end

  #Memory allocation for population data storage
  stageDataFrame = DataFrame(Week = 0, Stage1 = 0, Stage2 = 0, Stage3 = 0, Stage4 = 0, Total = 0)
  adultDataFrame = DataFrame(Week = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)

  #Memory allocation for mortality data storage (includes harvest)
  harvestDataFrame = DataFrame(Week = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)
  harvestZoneData = DataFrame(Week = 0, z1 = 0, z2 = 0, z3 = 0, z4 = 0, z5 = 0, z6 = 0, z7 = 0, z8 = 0, z9 = 0, z10 = 0,
                              z11 = 0, z12 = 0, z13 = 0, z14 = 0, z15 = 0, z16 = 0, z17 = 0, z18 = 0, Total = 0)
  killedDataFrame = DataFrame(Week = 0, Natural = 0, Extra = 0, Compensatory = 0, Total = 0)

  #Memory allocation for spawning data storage
  spawnDataFrame = DataFrame(Week = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)

  harvest_effort = fill(0., years)
  harvest_effort[1:length(effort)] = effort

  totalPopulation = sum(initStock)

  #initialize the progress meter
  if progress
    progressBar = Progress(years*52, 30, " $totalPopulation total agents, Year 1 (of $years), week 1 of simulation \n", 30)
  end

  # Specify the final week before beginning harvest season
  spawnMin = 39

  for y = 1:years
    for w = 1:52
      # Find total number of weeks passed in the simulation
      totalWeek = ((y-1)*52)+w

      # age specific population
      ageSpecificPop = fill(0, 7)
      for i = 1:length(a_db)
        for age = 2:8
          ageSpecificPop[age - 1] += getAgeSpecificPop(age, totalWeek, a_db[i].alive, a_db[i].weekNum, age_a)
        end #for age
      end #for i

      totalAdults = sum(ageSpecificPop)

      push!(adultDataFrame, vcat(totalWeek, ageSpecificPop..., sum(ageSpecificPop)))

      if progress
        progressBar.desc = " $totalAdults adults, $(ageSpecificPop[7]) age 8+, $totalPopulation total, Year $y (of $years), week $w of simulation "
        next!(progressBar)
      end

      # Harvest is all year but has peak periods
      harvest!(harvest_effort[y], totalWeek, a_db, e_a, adult_a, age_a, harvestDataFrame, harvestZoneData)

      # Spawn can be set to any week(s)
      if w > spawnMin
        spawn!(a_db, adult_a, age_a, e_a, totalWeek, carrying_capacity[y], spawnDataFrame)
      else
        push!(spawnDataFrame, (totalWeek, 0, 0, 0, 0, 0, 0, 0, 0))
      end

      #Agents are killed and moved weekly
      push!(killedDataFrame, (totalWeek, 0, 0, 0, 0))
      killAgeSpecific!(a_db, adult_a, ageSpecificPop, carrying_capacity[y], totalWeek, killedDataFrame)

      kill!(a_db, e_a, age_a, totalWeek, killedDataFrame)

      # Update agent locations weekly
      move!(a_db, age_a, e_a, totalWeek)

      #update population information
      stagePopulation = [0,0,0,0]; totalPopulation = 0;
      for j = 1:4
        stagePopulation[j] = getStagePopulation(j, totalWeek, a_db, age_a)
      end
      totalPopulation = sum(stagePopulation)

      push!(stageDataFrame, vcat(totalWeek,stagePopulation..., sum(stagePopulation)))

      #show a real time plot (every 10 weeks) of agent movement
      if plotPopDensity
        if w == 1 || w%10 == 0
          updatePopulationDensity!(a_db, popDensity, totalWeek)
          popPlot = spy(popDensity, Guide.title("Year = $y, week = $w, totalPop=$totalPopulation, totalAdult=$totalAdults"))
          display(popPlot)
        end
      end

      #show a real time plot (every 10 weeks) of adult age distribution
      if plotPopDistribution
        if w ==1 || w%10 == 0
          adultAge = [2,3,4,5,6,7,8]
          ageDistPlot = Gadfly.plot(x=adultAge, y=ageSpecificPop, Geom.point,
            Guide.xlabel("Adult age in years"), Guide.ylabel("Adult population"),
            Guide.title("Age Distribution Plot of Adult fish by age,\n y = $y, w = $w, totalPopulation = $totalPopulation"))
          display(ageDistPlot)
        end
      end

      #simulation failure protocol
      if totalPopulation == 0 || totalAdults > limit
        #simply for finishing the progress meter loops
        if progress
          for year = y:years
            for week = 1:52
              progressBar.desc = " $totalAdults adults, $totalPopulation total, Year $y (of $years), week $w of simulation "
              next!(progressBar)
            end
          end
        end

        removeEmptyClass!(a_db)
        description = simDescription*"\n\nSimulation population failure in year $y, week $w (total population
          = $totalPopulation, total adults = $totalAdults, population limit = $limit).\n\n"

        simSummary(adult_a, age_a, a_db, effort, ((length(carrying_capacity))*52), initStock, carrying_capacity,
                  stageDataFrame, adultDataFrame, harvestDataFrame, harvestZoneData, spawnDataFrame, killedDataFrame, description)
        return a_db
      end #population regulation failure
    end #end for week
    #Remove empty cohorts annually
    removeEmptyClass!(a_db)
  end #end for year

  description = string(simDescription, "\n\nSimulation was successfully completed.\n\n")

  simSummary(adult_a, age_a, a_db, effort, ((length(carrying_capacity))*52), initStock,
    carrying_capacity, stageDataFrame, adultDataFrame, harvestDataFrame, harvestZoneData,
    spawnDataFrame, killedDataFrame, description)

  return a_db
end


text:
function draw_axis!(ax_pos, specification, axis_attr = (;))
	@unpack P, df, args, layout_vars, group_dict, style_dict, kws, group_pairs, style_pairs = specification
	draw_axis!(P, ax_pos, df, args, layout_vars, group_dict, style_dict, kws, group_pairs, style_pairs, axis_attr)
end

function draw_axis!(P, figpos, df, args, layout_vars, group_dict, style_dict, kws, groups, styles, axis_attr)
	@unpack grp_x, grp_y, grp_wrap = layout_vars
	@unpack linkxaxes, linkyaxes, linkzcolor = layout_vars

	if isnothing(grp_wrap)
		I = isnothing(grp_y) ? 1 : length(unique(get(df, grp_y)))
		J = isnothing(grp_x) ? 1 : length(unique(get(df, grp_x)))
		N = I * J
	else
		N = length(unique(get(df, grp_wrap)))
		J = ceil(Int, sqrt(N))
		I = ceil(Int, N / J)
	end

	axs = [Axis(figpos[i,j]; axis_attr...) for i in 1:I, j in 1:J]
	for i in 1:I, j in 1:J
		ax = axs[i,j]
		if (i-1) * J + j > N
			hidespines!(ax)
			hidedecorations!(ax)
		end
	end

	grp = filter(!isnothing, collect((; grp_x, grp_y, grp_wrap)))

	out = combine(groupby(df, var_key.(grp))) do groupdf
		# Compute group key and index for layouting variables
		if !isnothing(grp_y)
			ykey = get(groupdf, grp_y) |> unique |> only
			i = categorical_positions(get(groupdf, grp_y)) |> unique |> only |> Int
		else
			i = 1
		end
		if !isnothing(grp_x)
			xkey = get(groupdf, grp_x) |> unique |> only
			j = categorical_positions(get(groupdf, grp_x)) |> unique |> only |> Int
		else
			j = 1
		end
		if !isnothing(grp_wrap)
			wrapkey = get(groupdf, grp_wrap) |> unique |> only
			ind = categorical_positions(get(groupdf, grp_wrap)) |> unique |> only |> Int

			i, j = fldmod1(ind, J)
		end

		# Do the plot
		plt = grouped_plot(P, axs[i, j], groupdf, group_dict, args, kws, groups, styles)

		let
			padding = (3f0, 3f0, 3f0, 3f0)
		# Add labels for faceting
		if !isnothing(grp_wrap)
			Box(  figpos[i, j, Top()], color=:lightgray)
			Label(figpos[i, j, Top()], string(wrapkey); padding)
		end
		if !isnothing(grp_x) && i == 1
			Box(  figpos[1, j, Top()], color=:lightgray)
			Label(figpos[1, j, Top()], string(xkey); padding)
		end
		if !isnothing(grp_y) && j == 1
	    	Box(  figpos[i, end, Right()], color=:lightgray)
			Label(figpos[i, end, Right()], string(ykey); padding, rotation = -pi/2)
		end
	end
		(; plt)#, color=:red)
	end

	# spanned labels
	span_label(:x, var_lab(args[1]), axs, figpos)
	span_label(:y, var_lab(args[2]), axs, figpos)

	# Link axes
	linkyaxes && linkyaxes!(axs...)
	linkxaxes && linkxaxes!(axs...)

    linkyaxes && hideydecorations!.(axs[:, 2:end], grid = false)

	if linkxaxes
		for i in 1:I, j in 1:J
			needs_deco = needs_xdecorations(i, j, J, N)
			if (i < I-1) || ((i == I-1) && !needs_deco)
				hidexdecorations!(axs[i,j], grid = false)
			end
			if (i < I) && needs_deco
				axs[i,j].alignmode = Mixed(bottom = MakieLayout.GridLayoutBase.Protrusion(0))
			end
		end
	end

	if linkzcolor && haskey(style_dict, :color)
		for p in out.plt
			if p.plt isa Vector
				for p_ in p.plt
					p_.colorrange = style_dict[:color]
				end
			else
				p.plt.colorrange = style_dict[:color]
			end
		end
	end

	# x and y labels
	#fig[1,1,Left()]   = Label(fig, string(y_var), rotation = pi/2, padding = (0,80,0,0))
	#fig[1,1,Bottom()] = Label(fig, string(x_var), padding = (0, 0, 0, 0))

end

# Check if there is a non-empty axis in the position just below
function needs_xdecorations(i, j, J, N)
	i * J + j > N
end

function span_label(axis, spanned_label, axs, figpos)

	positive_rotation = axis == :x ? 0f0 : π/2f0
	N_y, N_x = size(axs)

	itr = axis == :x ? axs[end, :] : axs[:, 1]
    group_protrusion = lift(
        (xs...) -> maximum(x -> axis == :x ? x.bottom : x.left, xs),
        (MakieLayout.protrusionsobservable(ax) for ax in itr)...
    )

    padding = lift(group_protrusion) do val
        val += 10f0
        axis == :x ? (0f0, 0f0, 0f0, val) : (0f0, val, 0f0, 0f0)
    end

    label = Label(figpos.fig, spanned_label, padding = padding, rotation = positive_rotation)
    pos = axis == :x ? (N_y, :, Bottom()) : (:, 1, Left())

	figpos[pos...] = label

end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
"""
    struct FileID <: Any
Hold a pointer to a `clang::FileID` object.

Note that, this ID is managed by source manager and should not be manually created.
"""
struct FileID
    ptr::CXFileID
end

Base.unsafe_convert(::Type{CXFileID}, x::FileID) = x.ptr
Base.cconvert(::Type{CXFileID}, x::FileID) = x

"""
    struct SourceLocation <: Any
Represent a Clang source location.

Note that, the underlying pointer is NOT a *pointer* to a `clang::SourceLocation` object.
Instead, it's the opaque pointer representation of the `clang::SourceLocation` itself.
"""
struct SourceLocation
    ptr::CXSourceLocation_
end

Base.unsafe_convert(::Type{CXSourceLocation_}, x::SourceLocation) = x.ptr
Base.cconvert(::Type{CXSourceLocation_}, x::SourceLocation) = x

"""
    struct SourceRange <: Any
Hold two `SourceLocation`s.
"""
struct SourceRange
    begin_loc::SourceLocation
    end_loc::SourceLocation
end


text:
"""
Uses Xfoil.jl to create a ::AirfoilPrep.Polar object for the specified airfoil.

Arguments:

* Re : Reynolds number at which a ::AirfoilPrep.Polar object is desired
* xs : vector containing x coordinates of the airfoil contour
* ys : vector containing y coordinates of the airfoil contour

Optional arguments:

* alphas              : angles of attack (DEGREES) at which to evaluate XFOIL for Polar construction
* M                   : desired Mach number
* `iter::Int`         : XFOIL convergence iterations, default to 100.
* `verbose::Bool`     : Verbose.

**NOTE: Airfoil points must go from trailing edge around the top, then the
bottom and end back at the trailing edge.**
"""
function runXFOIL(xs, ys, Re; alphas=range(-10, stop=25, step=0.5), Mach=0.0,
                                iter::Int=100, verbose=false,
                                npanels=160, ncrit=9, optargs...)

    xseps = zeros(2, length(alphas))

    # run XFOIL
    cls, cds, cdps, cms, convs = Xfoil.alpha_sweep(xs, ys, alphas, Re; iter=iter,
                                                    npan=npanels, mach=Mach,
                                                    ncrit=ncrit,
                                                    # xsep=xseps,
                                                    printdata=verbose,
                                                    optargs...)

    # remove unconverged points
    (alphas_converged, cls_converged,
        cds_converged, cms_converged,
        xseps_converged) = removeSingularities(alphas, cls, cds, cdps, cms, xseps, convs)

    # use AirfoilPrep to wrap airfoilpreppy
    polar = Polar(Re, alphas_converged, cls_converged, cds_converged, cms_converged;
                        Ma=Float64(Mach), npanels=npanels, ncrit=Float64(ncrit),
                            x=xs, y=ys, xsepup=xseps[1, :], xseplo=xseps[2, :])

    return polar
end

"""
Returns:

* alphas_converged  : vector of converged coefficients
* cls_converged : vector of converged coefficients
* cds_converged : vector of converged coefficients
* cdps_converged: vector of converged coefficients
* cms_converged : vector of converged coefficients

"""
function removeSingularities(alphas, cls, cds, cdps, cms, xseps, convs)
    iconverged = Bool.(convs)
    alphas_converged = alphas[iconverged]
    cls_converged = cls[iconverged]
    cds_converged = cds[iconverged]
    cdps_converged= cdps[iconverged]
    cms_converged = cms[iconverged]
    xseps_converged = xseps[:, iconverged]

    return alphas_converged, cls_converged, cds_converged, cms_converged, xseps_converged
end


text:
using Test
using Dictionaries
using Indexing

include("Indices.jl")
include("Dictionary.jl")
include("HashIndices.jl")
include("HashDictionary.jl")
include("PairDictionary.jl")
include("indexing.jl")
include("foreach.jl")
include("map.jl")
include("broadcast.jl")
include("filter.jl")
include("find.jl")


text:
abstract type DiscreteDynamics <: AbstractModel end
@inline evaluate(model::DiscreteDynamics, x, u, p) =
    discrete_dynamics(model, x, u, p.t, p.dt)
@inline evaluate!(model::DiscreteDynamics, xn, x, u, p) =
    discrete_dynamics!(model, xn, x, u, p.t, p.dt)

discrete_dynamics(model::DiscreteDynamics, z::AbstractKnotPoint) =
    discrete_dynamics(model, state(z), control(z), time(z), timestep(z))
discrete_dynamics!(model::DiscreteDynamics, xn, z::AbstractKnotPoint) =
    discrete_dynamics!(model, xn, state(z), control(z), time(z), timestep(z))

discrete_dynamics(model::DiscreteDynamics, x, u, t, dt) =
    error("Discrete dynamics not defined yet.")
discrete_dynamics!(model::DiscreteDynamics, xn, x, u, t, dt) =
    error("In-place discrete dynamics not defined yet.")

jacobian!(model::DiscreteDynamics, J, y, x, u, p) = jacobian!(model, J, y, x, u, p.t, p.dt)
jacobian!(model::DiscreteDynamics, J, y, x, u, t, dt) =
    error("User-defined discrete dynamics Jacobian not defined.")

dynamics_error(model::DiscreteDynamics, z2::AbstractKnotPoint, z1::AbstractKnotPoint) =
    discrete_dynamics(model, z1) - state(z2)
function dynamics_error!(
    model::DiscreteDynamics,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    discrete_dynamics!(model, y2, z1)
    y2 .-= state(z2)
    return nothing
end

dynamics_error_jacobian!(
    ::FunctionSignature,
    ::UserDefined,
    model::DiscreteDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error_jacobian!(model, J2, J1, y2, y1, z2, z1)

dynamics_error_jacobian!(
    model::DiscreteDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = error("User-defined dynamics error Jacobian not defined for $(typeof(model)).")

function propagate_dynamics!(
    ::InPlace,
    model::DiscreteDynamics,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    discrete_dynamics!(model, state(z2), z1)
    return nothing
end

function propagate_dynamics!(
    ::StaticReturn,
    model::DiscreteDynamics,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    setstate!(z2, discrete_dynamics(model, z1))
    return nothing
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Test
using FiniteDifferences
using StaticArrays
using LinearAlgebra
using Suppressor
using Random
using Unitful
using Plots

using OpticSim
# Geometry Module
using OpticSim.Geometry
# curve/surface imports
using OpticSim: findspan, makemesh, knotstoinsert, coefficients, inside, quadraticroots, tobeziersegments, evalcsg, makiemesh
# interval imports
using OpticSim: α, halfspaceintersection, positivehalfspace, lower, upper, EmptyInterval, rayorigininterval, intervalcomplement, intervalintersection, intervalunion, RayOrigin, Infinity, Intersection
include("TestData/TestData.jl")
# bounding box imports
using OpticSim: doesintersect
# RBT imports
using OpticSim: rotmat, rotmatd
# lens imports
using OpticSim: reflectedray, snell, refractedray, trace, intersection, pathlength, power
#Bounding volume hierarchy imports
# using OpticSim: partition!

const COMP_TOLERANCE = 25 * eps(Float64)
const RTOLERANCE = 1e-10
const ATOLERANCE = 10 * eps(Float64)
const SEED = 12312487
const ALL_TESTS = isempty(ARGS) || "all" in ARGS || "All" in ARGS || "ALL" in ARGS
const TESTSET_DIR = "testsets"

"""Evalute all functions not requiring arguments in a given module and test they don't throw anything"""
macro test_all_no_arg_functions(m)
    quote
        for n in names($m, all = true)
            # get all the valid functions
            if occursin("#", string(n)) || string(n) == string(nameof($m))
                continue
            end
            f = Core.eval($m, n)
            # get the methods of this function
            for meth in methods(f)
                # if this method has no args then try and evaluate it
                if meth.nargs == 1
                    # suppress STDOUT to prevent loads of stuff spamming the console
                    @test (@suppress_out begin
                        f()
                    end;
                    true)
                end
            end
        end
    end
end

include("Benchmarks/Benchmarks.jl")

alltestsets = [
    "JuliaLang",
    # "BVH",
    # "Examples", # slow
    "General",
    "SurfaceDefs",
    "Intersection",
    "Lenses",
    "OpticalSystem",
    "Emitters", # TODO
    "Comparison",
    "Visualization",
    "Allocations",
    "GlassCat",
    "Transform"
]

runtestsets = ALL_TESTS ? alltestsets : intersect(alltestsets, ARGS)
include.([joinpath(TESTSET_DIR, "$(testset).jl") for testset in runtestsets])


text:
module FeatureSelectors

using DataFrames
using Random

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Utils

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!, transform!
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

export feature_test

# generic way to extract num/cat features by specifying their columns
"""
    FeatureSelector(
       Dict(
         :name => "featureselector",
         :columns => [col1, col2, ...]
       )
    )

Returns a dataframe of the selected columns.

Implements `fit!` and `transform!`.
"""
mutable struct FeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function FeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "featureselector",
			    :columns => Int[],
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

"""
    FeatureSelector(cols::Vector{Int})

Helper function for FeatureSelector.
"""
FeatureSelector(cols::Vector{Int}) = FeatureSelector(Dict(:columns => cols))

"""
    FeatureSelector(cols::Vararg{Int})

Helper function for FeatureSelector.
"""
FeatureSelector(cols::Vararg{Int}) = FeatureSelector([cols...])

function fit!(ft::FeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
	error("empty dataframe")
    end
    ft.model = ft.args
end

function transform!(ft::FeatureSelector, features::DataFrame)
    nfeatures = deepcopy(features)
    if nfeatures == DataFrame()
	error("empty dataframe")
    end
    cols = ft.model[:columns]
    if  cols != []
	return nfeatures[:,cols]
    else
	return DataFrame()
    end
end

# ----------
"""
    CatFeatureSelector(Dict(:name => "catfeatsel"))

Automatically extract categorical columns based on
inferred element types.

Implements `fit!` and `transform!`.
"""
mutable struct CatFeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function CatFeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "catfeatsel",
			    :nominal_columns => []
		       )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function fit!(ft::CatFeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    catcols,_ = find_catnum_columns(features)

    # create model
    ft.model = Dict(
		    :nominal_columns => catcols
	       )
end

function transform!(ft::CatFeatureSelector, features::DataFrame)
    nfeatures = deepcopy(features)
    catcols = ft.model[:nominal_columns]
    if catcols != []
	return nfeatures[:,catcols]
    else
	return DataFrame()
    end
end

"""
    NumFeatureSelector(Dict(:name=>"numfeatsel"))

Automatically extracts numeric features based on their inferred element types.

Implements `fit!` and `transform!`.
"""
mutable struct NumFeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function NumFeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "numfeatsel",
			    :real_columns => []
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function fit!(ft::NumFeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    _,realcols = find_catnum_columns(features)

    # create model
    ft.model = Dict(
		    :real_columns => realcols
		    )
end

function transform!(ft::NumFeatureSelector, features::DataFrame)
    nfeatures = deepcopy(features)
    realcols = ft.model[:real_columns]
    if realcols != []
	return nfeatures[:,realcols]
    else
	return DataFrame()
    end
end

"""
    CatNumDiscriminator(
       Dict(
          :name => "catnumdisc",
          :maxcategories => 24
       )
    )

Transform numeric columns to string (as categories)
if the count of their unique elements <= maxcategories.

Implements `fit!` and `transform!`.
"""
mutable struct CatNumDiscriminator <: Transformer
    name::String
    model::Dict
    args::Dict

    function CatNumDiscriminator(args::Dict = Dict())
	default_args = Dict(
			    :name => "catnumdisc",
			    # default max categories for numeric-encoded categories
			    :maxcategories => 24,
			    :nominal_columns => Int[],
			    :real_columns => Int[]
		       )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

"""
    CatNumDiscriminator(maxcat::Int)

Helper function for CatNumDiscriminator.
"""
function CatNumDiscriminator(maxcat::Int)
    CatNumDiscriminator(Dict(:maxcategories=>maxcat))
end

function fit!(ft::CatNumDiscriminator, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    catcols,realcols = find_catnum_columns(features,ft.args[:maxcategories])

    # create model
    ft.model = Dict(
	:real_columns => realcols,
	:nominal_columns => catcols
    )
end

function transform!(ft::CatNumDiscriminator, features::DataFrame)
    nfeatures = features |> deepcopy
    catcols = ft.model[:nominal_columns]
    if catcols != []
		nfeatures[!,catcols] = nfeatures[!,catcols] .|> string
    end
    return nfeatures
end

end


text:
@json struct Note
    id::Int
    body::String
    attachment::String
    author::User
    created_at::DateTime
    updated_at::DateTime
    system::Bool
    noteable_id::Int
    noteable_type::String
    noteable_iid::Int
    resolvable::Bool
    confidential::Bool
end

## Merge Request Notes

# List All Merge Request Notes
# https://docs.gitlab.com/ee/api/notes.html#list-all-merge-request-notes
function endpoint(
    ::GitLabAPI,
    ::typeof(list_pull_request_comments),
    project::Integer,
    pull_request_id::Integer,
)
    return Endpoint(:GET, "/projects/$project/merge_requests/$pull_request_id/notes")
end
into(::GitLabAPI, ::typeof(list_pull_request_comments)) = Vector{Note}

# Get Single Merge Request Note
# https://docs.gitlab.com/ee/api/notes.html#get-single-merge-request-note
function endpoint(
    ::GitLabAPI,
    ::typeof(get_pull_request_comment),
    project::Integer,
    pull_request_id::Integer,
    comment_id::Integer,
)
    return Endpoint(
        :GET, "/projects/$project/merge_requests/$pull_request_id/notes/$comment_id"
    )
end
into(::GitLabAPI, ::typeof(get_pull_request_comment)) = Note

# Create Merge Request Note
# https://docs.gitlab.com/ee/api/notes.html#create-new-merge-request-note
function endpoint(
    ::GitLabAPI,
    ::typeof(create_pull_request_comment),
    project::Integer,
    pull_request_id::Integer,
)
    return Endpoint(:POST, "/projects/$project/merge_requests/$pull_request_id/notes")
end
into(::GitLabAPI, ::typeof(create_pull_request_comment)) = Note

# Update Merge Request Note
# https://docs.gitlab.com/ee/api/notes.html#modify-existing-merge-request-note
function endpoint(
    ::GitLabAPI,
    ::typeof(update_pull_request_comment),
    project::Integer,
    pull_request_id::Integer,
    comment_id::Integer,
)
    return Endpoint(
        :PUT, "/projects/$project/merge_requests/$pull_request_id/notes/$comment_id"
    )
end
into(::GitLabAPI, ::typeof(update_pull_request_comment)) = Note

# Delete Merge Request Note
# https://docs.gitlab.com/ee/api/notes.html#delete-a-merge-request-note
function endpoint(
    ::GitLabAPI,
    ::typeof(delete_pull_request_comment),
    project::Integer,
    pull_request_id::Integer,
    comment_id::Integer,
)
    return Endpoint(
        :DELETE, "/projects/$project/merge_requests/$pull_request_id/notes/$comment_id"
    )
end
into(::GitLabAPI, ::typeof(delete_pull_request_comment)) = Note


text:
@with_kw struct TMazeState
    x::Int64 = 1 # position in corridor
    g::Symbol = :north# goal north or south
end

@with_kw struct TMaze <: POMDP{Union{TMazeState,TerminalState}, Int64, Int64}
    n::Int64 = 10 # corridor length
    discount::Float64 = 0.99 # discount factor
end

# state space is (length of corr)*(north, south) + terminal
#                   |G|
# | | |x| | | | | | | |
#                   | |
function states(maze::TMaze)
    space = statetype(maze)[]
    for x in 1:(maze.n + 1), g in [:north, :south]
        push!(space, TMazeState(x, g))
    end
    push!(space, terminalstate) # terminal
    return space
end
stateindex(m::TMaze, s::TMazeState) = 2*s.x - (s.g==:north)
stateindex(m::TMaze, s::TerminalState) = 2*(m.n+1) + 1

# 4 actions: go North, East, South, West (1, 2, 3, 4)
actions(maze::TMaze) = 1:4
actionindex(maze::TMaze, i::Int) = i

# 5 observations: 2 for goal (left or right) + 2 for in corridor or at intersection + 1 term
observations(maze::TMaze) = 1:5
obsindex(maze::TMaze, i::Int) = i

function initialstate(maze::TMaze)
    s = states(maze)
    ns = length(s)
    p = zeros(ns) .+ 1.0 / (ns-1)
    p[end] = 0.0
    return SparseCat(s, p)
end

function transition(m::TMaze, s::TMazeState, a::Int64)
    if a == 1 || a == 3
        if s.x == m.n + 1
            Deterministic(terminalstate)
        else
            Deterministic(s)
        end
    elseif a == 2
        xp = min(s.x + 1, m.n + 1)
        return Deterministic(TMazeState(xp, s.g))
    elseif a == 4
        xp = max(s.x - 1, 1)
        return Deterministic(TMazeState(xp, s.g))
    end
end

transition(m::TMaze, s::TerminalState, a::Int64) = Deterministic(s)

function reward(m::TMaze, s::TMazeState, a::Int64)
    if s.x == m.n + 1
        # if at junction check action
        if (s.g == :north && a == 1) || (s.g == :south && a == 3)
            return 4.0
        else
            return -0.1
        end
    elseif a == 1 || a == 3
        # bump against wall
        return -0.1
    else
        return 0.0
    end
end

# observation mapping
#    1      2        3         4         5
# goal N  goal S  corridor  junction  terminal
function observation(m::TMaze, sp::TMazeState)
    if sp.x <= 2
        if sp.g == :north
            return Deterministic(1)
        else
            return Deterministic(2)
        end
    elseif sp.x == m.n+1
        return Deterministic(4)
    else
        return Deterministic(3)
    end
end

observation(m::TMaze, sp::TerminalState) = Deterministic(5)

initialobs(m::TMaze, s) = observation(m, s)

discount(m::TMaze) = m.discount

function POMDPs.convert_s(::Type{A}, s::Union{TMazeState,TerminalState}, m::TMaze) where A <: AbstractArray
    return convert(A, [stateindex(m, s)])
end

# inverse of stateindex(m::TMaze, s::TMazeState) = 2*s.x - (s.g==:north)
function POMDPs.convert_s(::Type{S}, v::AbstractVector, m::TMaze) where S <: Union{TMazeState,TerminalState}
    i = first(v)
    if i == 2*(m.n + 1) + 1
        return terminalstate
    end

    if i%2 == 0
        g = :south
    else
        g = :north
    end
    x = div(i-1, 2) + 1
    @assert x <= m.n + 1
    return TMazeState(x, g)
end

struct MazeBelief
    last_obs::Int64
    mem::Symbol # memory
end
MazeBelief() = MazeBelief(1, :none)

struct MazeUpdater <: Updater end
POMDPs.initialize_belief(bu::MazeUpdater, d::Any) = d

function POMDPs.update(bu::MazeUpdater, b::MazeBelief, a, o)
    mem = b.mem
    if o == 1
        mem = :north
    end
    if o == 2
        mem = :south
    end
    return MazeBelief(o, mem)
end

mutable struct MazeOptimal <: Policy end
POMDPs.updater(p::MazeOptimal) = MazeUpdater()

# 4 actions: go North, East, South, West (1, 2, 3, 4)
# observation mapping
#    1      2        3         4         5
# goal N  goal S  corridor  junction  terminal
function POMDPs.action(p::MazeOptimal, b::MazeBelief)
    # if don't know the goal go back
    if b.mem == :none
        return 4
    end
    if b.mem == :north && b.last_obs == 4
        return 1
    end
    if b.mem == :south && b.last_obs == 4
        return 3
    end
    return 2
end

# deprecated in POMDPs v0.9
POMDPs.initialstate_distribution(m::TMaze) = initialstate(m)


text:
module SampleJuliaPackage

greet() = print("Hello World!")

export test0, test0i, test1, test1i

function test0(a, b)
  return a[1] + b[1]
end

@inline function test0i(a, b)
  return a[1] + b[1]
end

function test1(a, b)
  @boundscheck (@assert a==b)
  return a[1] + b[1]
end

@inline function test1i(a ::Vector{S}, b ::Vector{T}) where {S, T}
  @boundscheck (@assert a==b)
  return a[1] + b[1]
end

@inline function test2i(a ::Vector{Int}, b ::Vector{Int})
  @boundscheck (@assert a==b)
  return a[1] + b[1]
end

end # module


text:
rot(n::Int) = M  -> reduce((M, _) -> rotr90(M), 1:n; init = M)
flip(M::Matrix)::Matrix = M[end:-1:1, :]

const TRANSFORMATIONS = [[rot(i) for i in 0:3];
                         [rot(i) ∘ flip for i in 0:3]]

function parse_input(file::AbstractString)
    tiles = open(file, "r") do io
        tiles = []
        while !eof(io)
            push!(tiles, parse_tile(io))
        end
        return tiles
    end

    tilemap = Dict()
    sidemap = Dict()
    for t in tiles, (tr_id, tr) in enumerate(TRANSFORMATIONS)
        t_rot_id = (id = t.id, tr = tr_id)
        t_rot = tr(t.tile)
        sides = tile_sides(t_rot)
        tilemap[t_rot_id] = (; inside = t_rot[2:end - 1, 2:end - 1], sides...)
        for side in pairs(sides)
            push!(get!(sidemap, side, NamedTuple[]), t_rot_id)
        end
    end
    return (tilemap, sidemap)
end

function parse_tile(io::IO)
    id_line = readline(io)
    id_m = match(r"^Tile ([0-9]+):$", id_line)
    @assert id_m !== nothing id_line
    id = parse(Int, id_m.captures[1])
    tile = hcat([split(readline(io), "") for _ in 1:10]...) |> permutedims
    @assert eof(io) || isempty(readline(io))
    return (; id, tile)
end

function tile_sides(tile::Matrix)
    side(v) = Symbol(join(v, ""))

    return (top    = side(tile[1, 1:end]),
            bottom = side(tile[end, 1:end]),
            left   = side(tile[1:end, 1]),
            right  = side(tile[1:end, end]))
end

up(ci::CartesianIndex) = ci - CartesianIndex(1, 0)
left(ci::CartesianIndex) = ci - CartesianIndex(0, 1)

function arrange(tiles, sides)::Matrix
    w = tiles |> keys .|> first |> unique |> length |> sqrt |> Int
    tileset = (unmatched_tops  = unmatched(tiles, :top),
               matched_bottoms = matched(tiles, :bottom),
               unmatched_lefts = unmatched(tiles, :left),
               matched_rights  = matched(tiles, :right),
               all             = keys(tiles),
               dict            = tiles)
    arr = Matrix{Any}(nothing, w, w)
    ok = arrange!(arr, keys(arr), [], tileset, sides)
    @assert ok
    println("found")
    return arr
end

unmatched(tiles, side::Symbol) = setdiff(keys(tiles), matched(tiles, side))

function matched(tiles, side::Symbol)
    opp_side::Symbol = opposite(side)
    return filter(keys(tiles)) do i
        tiles[i][side] ∈ [t[opp_side] for (j, t) in tiles if j.id != i.id]
    end
end

opposite(side::Symbol)::Symbol = opposite(Val(side))
opposite(::Val{:top}) = :bottom
opposite(::Val{:bottom}) = :top
opposite(::Val{:left}) = :right
opposite(::Val{:right}) = :left

function arrange!(arr::Matrix, cis, used, tileset, sides)::Bool
    isempty(cis) && return true
    (ci, cis...) = cis

    suitable = filter(∉(used) ∘ first, tileset.all)
    if ci[1] == 1
        intersect!(suitable, tileset.unmatched_tops)
    else
        border = tileset.dict[arr[up(ci)]].bottom
        intersect!(suitable, sides[:top => border])
    end
    if ci[1] < size(arr, 1)
        intersect!(suitable, tileset.matched_bottoms)
    end
    if ci[2] == 1
        intersect!(suitable, tileset.unmatched_lefts)
    else
        border = tileset.dict[arr[left(ci)]].right
        intersect!(suitable, sides[:left => border])
    end
    if ci[2] < size(arr, 2)
        intersect!(suitable, tileset.matched_rights)
    end

    for t in suitable
        arr[ci] = t
        arrange!(arr, cis, [used; t.id], tileset, sides) && return true
    end
    print(".")                  # backtrack
    return false
end

q1(arr) = map(t -> t.id, corners(arr)) |> prod

corners(M::Matrix)::Matrix = M[[1, end], [1, end]]

function monster(w::Int)::String
    m = ["..................#."
         "#....##....##....###"
         ".#..#..#..#..#..#..."]
    return join(m, repeat('.', w - length(m[1])))
end

function q2(arr, tiles)::Int
    return count(==('#'), mask_monsters(arr, tiles))
end

function mask_monsters(arr, tiles)::String
    images = assemble(arr, tiles)
    w = size(arr, 2) * size(first(tiles).second.inside, 2)
    m = monster(w)
    mre = Regex(m)
    image = images[only(findall(img -> count(mre, img) > 0, images))]
    mask = findall(==('#'), m)
    vimage = split(image, "")
    for i in findall(mre, image; overlap = true) .|> first
        @. vimage[i + mask - 1] = "O"
    end
    return join(vimage, "")
end

function assemble(arr, tiles)::Vector{String}
    grid = size(first(tiles).second.inside)
    image = fill("", size(arr) .* grid)
    for (i, t_id) in pairs(arr)
        tl = (i - one(i)) * grid + one(i)
        br  = i * grid
        image[tl:br] .= tiles[t_id].inside
    end
    return map(tr -> join(tr(image), ""), TRANSFORMATIONS)
end

Base.:*(ci::CartesianIndex, t::Tuple) = CartesianIndex(ci.I .* t)

let (tiles, sides) = parse_input("day20_test.in"),
    arr = arrange(tiles, sides)

    @assert q1(arr) == 20899048083289
    @assert q2(arr, tiles) == 273
end

(tiles, sides) = parse_input("day20.in")
arr = arrange(tiles, sides)
println("Q1: ", q1(arr))
println("Q2: ", q2(arr, tiles))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using LinearAlgebra
@testset "biconnected" begin
    n = 10
    O = ones(Int64,n-1)
    Z = zeros(Int64,n)
    A = Tridiagonal(O,Z,O)
    B = MatrixNetwork(dropzeros!(sparse(A)))
    obj = biconnected_components(B)

    for i = 2:n-1
        @test obj.articulation_points[i]!=0
    end
    components = biconnected_components(B;art=false, components=false).number
    number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
    @test 2*components == 2*number_of_components == length(obj.map)

    @test obj.articulation_points[1]==0
    @test obj.articulation_points[n]==0

    B = dropzeros!(sparse(A))
    obj = biconnected_components(B)

    for i = 2:n-1
        @test obj.articulation_points[i]!=0
    end
    components = biconnected_components(B;art=false, components=false).number

    @test obj.articulation_points[1]==0
    @test obj.articulation_points[n]==0

    B = sparse_to_csr(dropzeros!(sparse(A)))
    obj = biconnected_components(B...)
    for i = 2:n-1
        @test obj.articulation_points[i]!=0
    end
    components = biconnected_components(B...;art=false, components=false).number

    @test obj.articulation_points[1]==0
    @test obj.articulation_points[n]==0

    B = findnz(dropzeros!(sparse(A)))
    obj = biconnected_components(B[1], B[2])

    for i = 2:n-1
        @test obj.articulation_points[i]!=0
    end
    components = biconnected_components(B[1], B[2];art=false, components=false).number

    @test obj.articulation_points[1]==0
    @test obj.articulation_points[n]==0

    A = load_matrix_network("minnesota")
    B = MatrixNetwork(A)
    components = biconnected_components(B;art=false, components=false).number
    number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
    @test components == number_of_components == 142

    A = load_matrix_network("biconnected_example")
    B = MatrixNetwork(A)
    obj = biconnected_components(B)
    components = biconnected_components(B;art=false, components=false).number
    number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
    @test components == number_of_components == 5
    @test obj.articulation_points[5] == 1

    A = load_matrix_network("cores_example")
    B = MatrixNetwork(A)
    obj = biconnected_components(B)
    components = biconnected_components(B;art=false, components=false).number
    number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
    @test components == number_of_components == 6
    @test obj.articulation_points[20] == 1
    @test obj.articulation_points[2] == 1
    @test obj.articulation_points[11] == 1
    @test obj.articulation_points[12] == 0

    A = empty_graph()
    obj = biconnected_components(A)
    components = biconnected_components(A;art=false, components=false).number
    number_of_components = biconnected_components!(A, zeros(Bool,0), zeros(Int64,0))
    @test components == number_of_components == 0

    A = load_matrix_network("clique-10")
    B = MatrixNetwork(A)
    obj = biconnected_components(B)
    components = biconnected_components(B;art=false, components=false).number
    number_of_components = biconnected_components!(B, zeros(Bool,0), zeros(Int64,0))
    @test components == number_of_components == 1
end


text:
module Models

using ..OceanTurbulenceParameterEstimation
using ..OceanTurbulenceParameterEstimation.Grids
using ..OceanTurbulenceParameterEstimation.Observations

using Oceananigans
using Oceananigans: AbstractModel
using Oceananigans.Fields
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractEddyViscosityClosure
using Oceananigans.Models

using Printf

import Oceananigans.Fields: interpolate
import Base: size
import OceanTurbulenceParameterEstimation.Observations: set!

export
       get_model_field,
       set!,

       # ensemble_model.jl
       OneDimensionalEnsembleModel,
       ensemble_size, batch_size,
       initialize_forward_run!,

       # free_parameters.jl
       DefaultFreeParameters,
       get_free_parameters,
       FreeParameters,
       @free_parameters

#
# AbstractModel extensions
#

function get_model_field(m::AbstractModel, p)
    p ∈ propertynames(m.tracers) && return m.tracers[p]
    p ∈ propertynames(m.velocities) && return m.velocities[p]
    @error "$p is not a valid field name"
end

"""
set!(model; kwargs...)

Set velocity and tracer fields of `model` given data defined on `grid`.
The keyword arguments `kwargs...` take the form `name=data`, where `name` refers
to one of the fields of `model.velocities` or `model.tracers`, and the `data`
may be an array, a function with arguments `(x, y, z)`, or any data type for which a
`set!(ϕ::AbstractField, data)` function exists.

Equivalent to `Oceananigans.Models.HydrostaticFreeSurfaceModels.set!` but where the
data in `kwargs...` exist on an arbitrary-resolution grid `grid`
and may need to be interpolated to the model grid `model.grid`.
"""
function set!(model; kwargs...)

    for (fldname, value) in kwargs

        # Get the field `ϕ` corresponding to `fldname`.
        if fldname ∈ propertynames(model.velocities)
            ϕ = getproperty(model.velocities, fldname)
        elseif fldname ∈ propertynames(model.tracers)
            ϕ = getproperty(model.tracers, fldname)
        elseif fldname ∈ propertynames(model.free_surface)
            ϕ = getproperty(model.free_surface, fldname)
        else
            throw(ArgumentError("name $fldname not found in model.velocities or model.tracers."))
        end

        set!(ϕ, value)
    end

    return nothing
end

include("ensemble_model.jl")
include("free_parameters.jl")

end #module


text:
using GSL
using HypothesisTests
using Distributions
#=
  This test below is a port from MethPipe methdiff package.

  Which is a one directional Fishers test
    https://groups.google.com/forum/#!topic/methpipe/dBWjZTZjlV0

  In https://github.com/smithlabcode/methpipe/blob/master/src/radmeth/methdiff.cpp a pseudocount
  is added before calculating test_greater_population_s. Likely to avoid zero counts. If there
  are zero counts, this returns 1 - rather than error.

  Currently, I don't get values of a similar magnitude between this implementation and the
  Julia HypothesisTests (and R) based right sided Fisher's exact test
=#

function log_sum_log(p, q)
   if p == 0
     return q
   elseif q == 0
     return p
   end
   larger = (p > q) ? p : q
   smaller = (p > q) ? q : p
   return larger + log(1.0 + exp(smaller - larger))
end

function log_hyper_g_greater(meth_a, unmeth_a, meth_b, unmeth_b, k)
   return  sf_lnchoose(meth_b + unmeth_b - 1, k) +
     sf_lnchoose(meth_a + unmeth_a - 1, meth_a + meth_b - 1 - k) -
     sf_lnchoose(meth_a + unmeth_a + meth_b + unmeth_b - 2, meth_a + meth_b - 1)
end

function test_greater_population_s( meth_a, unmeth_a, meth_b, unmeth_b)
   p = 0;
   k_begin = (meth_b > unmeth_a) ? meth_b - unmeth_a : 0
   for k=k_begin:(meth_b-1)
        p = log_sum_log(p, log_hyper_g_greater(meth_a, unmeth_a, meth_b, unmeth_b, k))
   end
   return exp(p)
end

#=

test_greater_population(meth_a, unmeth_a, meth_b, unmeth_b; ci=false )

Uses the right sided Fisher's Test using HypthesisTests' FisherExactTest

       | Cond A | Cond B
---------------------------
  meth |   a     |  b
unmeth |   c     |  d

Test whether a/c == b/d using FisherExactTest(a,b,c,d) - so we have
 FisherExactTest(A_meth,B_meth,A_unmeth,B_unmeth)

Use the one sided test since our alternative is that A is greater than B.

=#
function test_greater_population_fisher(meth_a, unmeth_a, meth_b, unmeth_b; ci=false )
    FE=nothing
    try
        FE=FisherExactTest(meth_a, meth_b, unmeth_a, unmeth_b)
    catch e
         error("Caught error: $e. counts: $meth_a, $unmeth_a, $meth_b, $unmeth_b")
    end
    return pvalue(FE, tail=:right)
end

#=
 The same as above, except, when  HypothesisTest's FisherTest is instantiated it computes
 the odss ratio that currently is compute intensive for counts > 100 or more. However,
 we only want the p value for the one sided test - so steal this from the HypothesisTests package
=#

function test_greater_population(meth_a, unmeth_a, meth_b, unmeth_b)
     a = meth_a
     b = meth_b
     c = unmeth_a
     d = unmeth_b
     p = pvalue(Hypergeometric(a + b, c + d, a + c), a, tail=:right)
     p = max(min(p, 1.0), 0.0)
     return p
end


text:
using SerialPorts

function _end_of_command(ser)
    for i in 0:2
        write(stdout, 0xff)
        write(ser, 0xff)
    end
end

function _execute_command(ser, cmd)
    println(cmd)
    write(ser, cmd)
    _end_of_command(ser)
end

function main()
    ser = SerialPort("COM4", 9600)
    _execute_command(ser, "page 0")
    sleep(2)
    _execute_command(ser, "page 1")
    sleep(2)
    _execute_command(ser, "t0.txt=\"Hello\"")
    close(ser)
end

main()