text:
"""
    covsel(X, Y; nlv = nothing, scaly = true)
Variable (feature) selection with the covsel method
* `X` : X-data.
* `Y` : Y-data.
* `nlv` : Nb. variables to select.
* `scaly` : If `true`, columns of `Y` are scaled.

The function alse returns the absolute partial covariances (`covs`)
for the selected variables.

## References

Höskuldsson, A., 1992. The H-principle in modelling with applications
to chemometrics. Chemometrics and Intelligent Laboratory Systems,
Proceedings of the 2nd Scandinavian Symposium on Chemometrics 14,
139–153. https://doi.org/10.1016/0169-7439(92)80099-P

Roger, J.M., Palagos, B., Bertrand, D., Fernandez-Ahumada, E., 2011.
covsel: Variable selection for highly multivariate and multi-response
calibration: Application to IR spectroscopy.
Chem. Lab. Int. Syst. 106, 216-223.
"""
function covsel(X, Y; nlv = nothing, scaly = true)
    covsel!(copy(X), copy(Y); nlv = nlv, scaly = scaly)
end

function covsel!(X, Y; nlv = nothing, scaly = true)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    n, p = size(X)
    isnothing(nlv) ? nlv = p : nothing
    xmeans = colmeans(X)
    ymeans = colmeans(Y)
    center!(X, xmeans)
    center!(Y, ymeans)
    if scaly
        yscales = sqrt.(sum(Y.^2; dims = 1))
        Y = scale(Y, yscales)
    end
    xsstot = sum(X.^2)
    ysstot = sum(Y.^2)
    xss = zeros(nlv)
    yss = zeros(nlv)
    selvar = Int64.(zeros(nlv))
    selcov = zeros(nlv)
    covs = zeros(p)
    H = similar(X, n, n)
    for i = 1:nlv
        z = vec(sum(abs.(X' * Y); dims = 2))
        zsel = argmax(z)
        selvar[i] = zsel
        selcov[i] = z[zsel]
        covs[zsel] = z[zsel]
        x = vcol(X, zsel)
        H .= x * x' / sum(x.^2)
        X .= X .- H * X
        Y .= Y .- H * Y
        xss[i] = sum(X.^2)
        yss[i] = sum(Y.^2)
    end
    cumpvarx = 1 .- xss / xsstot
    cumpvary = 1 .- yss / ysstot
    sel = DataFrame((sel = selvar, cov = selcov,
        cumpvarx = cumpvarx, cumpvary = cumpvary))
    (sel = sel, covs = covs)
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#gravity.jl

#2-dimensional gravity simulation to show the limits of floating point arithmetic.

type Body{T <: Real}
  pos::Vector{T}
  vel::Vector{T}
  mass::T
end

#just use units where the gravitational constant is 1.

sqr(x::Float32) = x^2
sqr(x::Float64) = x^2
sqr(x::BigFloat) = x^2

#subjects object a to the gravitational influence of object b.
#note that only a is changed.

#gravitational constant.
Gc_std = 6.674e-11 #N m^2 kg^-2 = m^3 s^-2 kg^-1
#rescale distance to 10^9 m (1 AU = ~150 units)
# = 6.674e-38 d^3 s^-2 kg^-1
#rescale time to hours
# = 8.650e-31 d^3 h^-2 kg^-1
#rescale mass to 10^24 kg (1Me = 6.0 Mu)
# = 8.650e-7 d^3 h^-2 Mu-1
Gc = 8.650e-7

function gravity!{T}(a::Body{T}, b::Body{T})
  #calculate the displacement vector:
  disp = b.pos - a.pos
  sdisp = map(sqrt, sum(map(sqr, disp)))
  a.vel += Gc * disp * (b.mass / (sdisp) ^ 3)
  nothing
end

#executes the passage of time on object a.
function time!{T}(a::Body{T})
  a.pos += a.vel
  nothing
end

function adjust_system!{T}(a::Vector{Body{T}})
  total_mass = mapreduce((b) -> b.mass, + , zero(T), a)
  #number one, set the center of mass to be at the origin.
  delta_p::Vector{T} = mapreduce((b) -> b.pos * b.mass, +, zero(T), a) / total_mass
  map!((b)-> Body(b.pos - delta_p, b.vel, b.mass), a)
  #number two, set the frame of reference so that net momentum is zero.
  delta_v::Vector{T} = mapreduce((b) -> b.vel * b.mass, +, zero(T), a) / total_mass
  map!((b)-> Body(b.pos, b.vel - delta_v, b.mass), a)
  nothing
end

#applies gravity to a system s.
function apply_gravity!{T}(s::Array{Body{T}})
  l = length(s)
  for idx = 1:l
    for jdx = 1:l
      (idx != jdx) && gravity!(s[idx], s[jdx])
    end
  end
  nothing
end

function apply_time!{T}(s::Array{Body{T}})
  for idx = 1:length(s)
    time!(s[idx])
  end
end

#speed of earth:  ~ 150 d * 3.14 * 2 / 360 / 24

star = Body{Float64}([0.0,0.0], [0.0, 0.0], 1998000)
planet = Body{Float64}([0.0, -150.0], [0.11, 0.0], 6)

system = Body{Float64}[star, planet]
adjust_system!(system)

for idx=1:1000
  for idx = 1:100_000
    apply_gravity!(system)
    apply_time!(system)
  end
  println(system[2].pos[1], " " , system[2].pos[2], " r: ", sqrt(sum(map(sqr,system[2].pos))))
end

ext: jl

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Various physical constants
const G_FERMI = 0.000011663787;
const HIGGS_VEV = 246.21965;
const ALPHA_EM = 0.0078125; # at p^2 = 0
const SIN_THETA_WEAK = 0.480853;
const SIN_THETA_WEAK_SQRD = 0.23122;
const COS_THETA_WEAK = 0.876801;
const M_PLANK = 1.220890e19;
const RHO_CRIT = 0.0000105375;
const S_TODAY = 2891.2;
const T_CMB = 2.56215e-10;
const T_BBN = 0.0001; # 0.1 MeV in GeV
const OMEGA_H2_CDM = 0.1198;

# Masses
const ELECTRON_MASS = 0.0005109989461;
const MUON_MASS = 0.1056583745;
const TAU_MASS = 1.77686;
const UP_QUARK_MASS = 0.00216;
const DOWN_QUARK_MASS = 0.00467;
const STRANGE_QUARK_MASS = 0.093;
const CHARM_QUARK_MASS = 1.27;
const BOTTOM_QUARK_MASS = 4.18;
const TOP_QUARK_MASS = 172.9;
const W_BOSON_MASS = 80.379;
const Z_BOSON_MASS = 91.1876;
const HIGGS_MASS = 125.10;
const NEUTRAL_PION_MASS = 134.9766e-3;
const CHARGED_PION_MASS = 139.57018e-3;
const NEUTRAL_KAON_MASS = 497.61e-3;
const LONG_KAON_MASS = 497.614e-3;
const CHARGED_KAON_MASS = 493.68e-3;

# Boson widths
const W_BOSON_WIDTH = 2.085;
const Z_BOSON_WIDTH = 2.4952;
const HIGGS_WIDTH = 0.00407;

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
struct FilterSimilarBatchQs <: BatchPQs
    model::SVDD.OCClassifier
    sequential_qs::SequentialPQs
    div_measure::Function
    k::Int

    function FilterSimilarBatchQs(model::SVDD.OCClassifier, sequential_qs::SequentialPQs; diversity::Symbol, k::Int)::TopKBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))
        diversity_measure = get_pairwise_distance(diversity)
        return new(model, sequential_qs, diversity_measure, k)
    end
end

"""
Iteratively find most similar pairs and prune the one with less utility calculated
with a sequential strategy. This process is repeated until batch_size observations
are left.
"""
function select_batch(qs::FilterSimilarBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    candidate_scores = qs_score(qs.sequential_qs, x, labels)[candidate_indices]
    distances = qs.div_measure(qs.model, x, candidate_indices)

    mask = trues(num_observations)
    while sum(mask) > qs.k
        a, b = Tuple(findmin(distances[mask, mask])[2])
        eliminate_idx = candidate_scores[a] > candidate_scores[b] ? b : a
        mask[findall(mask)[eliminate_idx]] = false
    end

    batch = candidate_indices[findall(mask)]
    return batch
end


text:
#
# x, y: Arrays of dimensions (nx,3), (ny,3)
# ifmol, ilmol: first and last indexes of array x to be considered (first and last atoms of molecule)
# jfmol, jlmol: first and last indexes of array y to be considered
# returns the the minimum distance between the points considered and
# the indexes of these points in x and y vectors
#

# If x is only a vector
"""

Returns the the minimum distance between the points considered.

"""
function minimumdistance(x::T, y::Vector{T}) where T
  jatom = 0
  dmin = +Inf
  for (j,yj) in pairs(y)
    d = distance(x,yj)
    if d < dmin
      jatom = j
      dmin = d
    end
  end
  return dmin, 1, jatom
end

# If both are vectors of vectors

function minimumdistance(x::Vector{T}, y::Vector{T}) where T
  iatom = 0
  jatom = 0
  dmin = +Inf
  for (i,xi) in pairs(x)
     for (j,yj) in pairs(y)
       d = distance(xi,yj)
       if d < dmin
         iatom = i
         jatom = j
         dmin = d
       end
     end
  end
  return dmin, iatom, jatom
end

# Function that returns the distance of a reference atom as well, to be used for
# computation of the volume shell by Monte-Caro integration

function minimumdistance(x::Vector{T}, y::Vector{T},
                         jrefatom::Int) where T
  iatom = 0
  jatom = 0
  drefatom = +Inf
  dmin = +Inf
  for (i,xi) in pairs(x)
     for (j,yj) in pairs(y)
       d = distance(xi,yj)
       # Minimum distance of any solvent atom to the solute
       if d < dmin
         iatom = i
         jatom = j
         dmin = d
       end
       # Minimum distance of the reference atom to the solute
       if j == jrefatom
         if d < drefatom
           drefatom = d
         end
       end
     end
  end
  return dmin, iatom, jatom, drefatom
end

#
# With periodic boundary conditions
#

function minimumdistance(x::Vector{T}, y::Vector{T},
                         sides::T) where T
  iatom = 0
  jatom = 0
  dmin = +Inf
  for (i,xi) in pairs(x)
     for (j,yj) in pairs(y)
       d = distance(xi,yj,sides)
       if d < dmin
         iatom = i
         jatom = j
         dmin = d
       end
     end
  end
  return dmin, iatom, jatom
end

# Function that returns the distance of a reference atom as well, to be used for
# computation of the volume shell by Monte-Caro integration

function minimumdistance(x::Vector{T}, y::Vector{T},
                         jrefatom::Int, sides::T) where T
  iatom = 0
  jatom = 0
  drefatom = +Inf
  dmin = +Inf
  for (i,xi) in pairs(x)
     for (j,yj) in pairs(y)
       d = distance(xi,yj,sides)
       # Minimum distance of any solvent atom to the solute
       if d < dmin
         iatom = i
         jatom = j
         dmin = d
       end
       # Minimum distance of the reference atom to the solute
       if j == jrefatom
         if d < drefatom
           drefatom = d
         end
       end
     end
  end
  return dmin, iatom, jatom, drefatom
end


text:
include("header.jl")
Random.seed!(42)
struct M370; layer; end;

@testset "conv" begin

    conv41(a;o...)=conv4(a[1],a[2];o...)
    deconv41(a;o...)=deconv4(a[1],a[2];o...)
    rand41(d...)=reshape(0.01*collect(Float64,1:prod(d)),d)

    TOL=0.2
    ax = rand41(5,4,3,2)
    aw = rand41(3,3,3,4)
    ad = permutedims(aw, (1,2,4,3))
    ax32 = convert(Array{Float32}, ax)
    aw32 = convert(Array{Float32}, aw)
    ad32 = convert(Array{Float32}, ad)
    ax5 = rand41(6,5,4,3,2)
    aw5 = rand41(3,3,3,3,3)
    if gpu() >= 0
        kx = KnetArray(ax)
        kw = KnetArray(aw)
        kd = KnetArray(ad)
        kx32 = KnetArray(ax32)
        kw32 = KnetArray(aw32)
        kd32 = KnetArray(ad32)
        kx5 = KnetArray(ax5)
        kw5 = KnetArray(aw5)
    end

    @warn "cpuconv tests commented out until we figure out why they are failing on gitlab-ci"

#=
    @testset "cpuconv" begin
        ### Default
        @test gradcheck(pool, ax)
        @test gradcheck(unpool, ax)
        @test isapprox(pool(unpool(ax)),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL)
        @test gradcheck(deconv41, (ad,ax); rtol=TOL)

        ### Float32
        @test gradcheck(pool, ax32; rtol=TOL) # TODO: sensitive to seed
        @test gradcheck(unpool, ax32; rtol=TOL) # TODO: sensitive to seed
        @test isapprox(pool(unpool(ax32)),ax32)
        @test gradcheck(conv41, (aw32,ax32); rtol=0.5) # TODO: sensitive to seed
        @test gradcheck(deconv41, (ad32,ax32); rtol=TOL) # TODO: sensitive to seed

        ### 5D
        @test gradcheck(pool, ax5)
        @test gradcheck(unpool, ax5)
        @test isapprox(pool(unpool(ax5)),ax5)
        @test gradcheck(conv41, (aw5,ax5); rtol=TOL)
        @test gradcheck(deconv41, (aw5,ax5); rtol=TOL)

        ### window=3 (default=2) only for pool
        @test gradcheck(pool, ax; kw=[(:window,3)])
        @test gradcheck(unpool, ax; kw=[(:window,3)])
        @test isapprox(pool(unpool(ax;window=3);window=3),ax)
        @test gradcheck(pool, ax; kw=[(:window,(3,3))])
        @test gradcheck(unpool, ax; kw=[(:window,(3,3))])
        @test isapprox(pool(unpool(ax;window=(3,3));window=(3,3)),ax)

        ### padding=1 (default=0)
        @test gradcheck(pool, ax; kw=[(:padding,1)])
        @test gradcheck(unpool, ax; kw=[(:padding,1)])
        @test isapprox(pool(unpool(ax;padding=1);padding=1),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:padding,1)])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:padding,1)])
        @test gradcheck(pool, ax; kw=[(:padding,(1,1))])
        @test gradcheck(unpool, ax; kw=[(:padding,(1,1))])
        @test isapprox(pool(unpool(ax;padding=(1,1));padding=(1,1)),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:padding,(1,1))])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:padding,(1,1))])

        ### stride=3 (default=1 for conv, window=2 for pool)
        @test gradcheck(pool, ax; kw=[(:stride,3)])
        @test gradcheck(unpool, ax; kw=[(:stride,3)])
        @test isapprox(pool(unpool(ax;stride=3);stride=3),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:stride,3)])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:stride,3)])
        @test gradcheck(pool, ax; kw=[(:stride,(3,3))])
        @test gradcheck(unpool, ax; kw=[(:stride,(3,3))])
        @test isapprox(pool(unpool(ax;stride=(3,3));stride=(3,3)),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:stride,(3,3))])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:stride,(3,3))])

        ### mode=1 (default=0)
        @test gradcheck(pool, ax; kw=[(:mode,1),(:padding,1)])
        @test gradcheck(unpool, ax; kw=[(:mode,1),(:padding,1)])
        @test isapprox(pool(unpool(ax;mode=1);mode=1),ax)
        @test_broken isapprox(pool(unpool(ax;mode=1,padding=1);mode=1,padding=1),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:mode,1),(:padding,1)])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:mode,1),(:padding,1)])

        ### mode=2 (only for pool) -- is not supported in NNlib #218
        # @test gradcheck(pool, ax; kw=[(:mode,2),(:padding,1)])
        # @test gradcheck(unpool, ax; kw=[(:mode,2),(:padding,1)])
        # @test isapprox(pool(unpool(ax;mode=2);mode=2),ax)
        # @test isapprox(pool(unpool(ax;mode=2,padding=1);mode=2,padding=1),ax)

        ### alpha=2 (default=1)
        @test gradcheck(pool, ax; kw=[(:alpha,2)])
        @test gradcheck(unpool, ax; kw=[(:alpha,2)])
        @test isapprox(pool(unpool(ax;alpha=2);alpha=2),ax)
        @test gradcheck(pool, ax; kw=[(:alpha,2),(:mode,1),(:padding,1)])
        @test gradcheck(unpool, ax; kw=[(:alpha,2),(:mode,1),(:padding,1)])
        @test isapprox(pool(unpool(ax;alpha=2,mode=1);alpha=2,mode=1),ax)
        @test_broken isapprox(pool(unpool(ax;alpha=2,mode=1,padding=1);alpha=2,mode=1,padding=1),ax)
        @test gradcheck(conv41, (aw,ax); rtol=TOL, kw=[(:alpha,2)])
        @test gradcheck(deconv41, (ad,ax); rtol=TOL, kw=[(:alpha,2)])
    end
=#

    if gpu() >= 0; @testset "gpuconv" begin
        ### Default
        @test isapprox(pool(kx), pool(ax))
        @test gradcheck(pool, kx)
        @test isapprox(unpool(kx), unpool(ax))
        @test gradcheck(unpool, kx)
        @test isapprox(conv4(kw,kx), conv4(aw,ax))
        @test gradcheck(conv41, (kw,kx); rtol=TOL)
        @test isapprox(deconv4(kd,kx), deconv4(ad,ax))
        @test gradcheck(deconv41, (kd,kx); rtol=TOL)

        ### Float32
        @test isapprox(pool(kx32), pool(ax32))
        @test gradcheck(pool, kx32)
        @test isapprox(unpool(kx32), unpool(ax32))
        @test gradcheck(unpool, kx32)  # TODO: sensitive to seed
        @test isapprox(conv4(kw32,kx32), conv4(aw32,ax32))
        @test gradcheck(conv41, (kw32,kx32); rtol=TOL)
        @test isapprox(deconv4(kd32,kx32), deconv4(ad32,ax32))
        @test gradcheck(deconv41, (kd32,kx32); rtol=TOL)

        ### 5D
        @test isapprox(pool(kx5), pool(ax5))
        @test gradcheck(pool, kx5)
        @test isapprox(unpool(kx5), unpool(ax5))
        @test gradcheck(unpool, kx5)
        @test isapprox(conv4(kw5,kx5), conv4(aw5,ax5))
        @test gradcheck(conv41, (kw5,kx5); rtol=TOL)
        @test isapprox(deconv4(kw5,kx5), deconv4(aw5,ax5))
        @test gradcheck(deconv41, (kw5,kx5); rtol=TOL)

        ### window=3 (default=2) only for pool
        @test isapprox(pool(kx;window=3), pool(ax;window=3))
        @test gradcheck(pool, kx; kw=[(:window,3)])
        @test isapprox(unpool(kx;window=3), unpool(ax;window=3))
        @test gradcheck(unpool, kx; kw=[(:window,3)])
        @test isapprox(pool(kx;window=(3,3)), pool(ax;window=(3,3)))
        @test gradcheck(pool, kx; kw=[(:window,(3,3))])
        @test isapprox(unpool(kx;window=(3,3)), unpool(ax;window=(3,3)))
        @test gradcheck(unpool, kx; kw=[(:window,(3,3))])

        ### padding=1 (default=0)
        @test isapprox(pool(kx;padding=1), pool(ax;padding=1))
        @test gradcheck(pool, kx; kw=[(:padding,1)])
        @test isapprox(unpool(kx;padding=1), unpool(ax;padding=1))
        @test gradcheck(unpool, kx; kw=[(:padding,1)])
        @test isapprox(conv4(kw,kx;padding=1), conv4(aw,ax;padding=1))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:padding,1)])
        @test isapprox(deconv4(kd,kx;padding=1), deconv4(ad,ax;padding=1))
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:padding,1)])

        @test isapprox(pool(kx;padding=(1,1)), pool(ax;padding=(1,1)))
        @test gradcheck(pool, kx; kw=[(:padding,(1,1))])
        @test isapprox(unpool(kx;padding=(1,1)), unpool(ax;padding=(1,1)))
        @test gradcheck(unpool, kx; kw=[(:padding,(1,1))])
        @test isapprox(conv4(kw,kx;padding=(1,1)), conv4(aw,ax;padding=(1,1)))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:padding,(1,1))])
        @test isapprox(deconv4(kd,kx;padding=(1,1)), deconv4(ad,ax;padding=(1,1)))
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:padding,(1,1))])

        ### stride=3 (default=1 for conv, window=2 for pool)
        @test isapprox(pool(kx;stride=3), pool(ax;stride=3))
        @test gradcheck(pool, kx; kw=[(:stride,3)])
        @test isapprox(unpool(kx;stride=3), unpool(ax;stride=3))
        @test gradcheck(unpool, kx; kw=[(:stride,3)])
        @test isapprox(conv4(kw,kx;stride=3), conv4(aw,ax;stride=3))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:stride,3)])
        @test isapprox(deconv4(kd,kx;stride=3), deconv4(ad,ax;stride=3); rtol=1e-6)
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:stride,3)])

        @test isapprox(pool(kx;stride=(3,3)), pool(ax;stride=(3,3)))
        @test gradcheck(pool, kx; kw=[(:stride,(3,3))])
        @test isapprox(unpool(kx;stride=(3,3)), unpool(ax;stride=(3,3)))
        @test gradcheck(unpool, kx; kw=[(:stride,(3,3))])
        @test isapprox(conv4(kw,kx;stride=(3,3)), conv4(aw,ax;stride=(3,3)))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:stride,(3,3))])
        @test isapprox(deconv4(kd,kx;stride=(3,3)), deconv4(ad,ax;stride=(3,3)); rtol=1e-6)
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:stride,(3,3))])

        ### mode=1 (default=0)
        @test isapprox(pool(kx;mode=1,padding=1), pool(ax;mode=1,padding=1))
        @test gradcheck(pool, kx; kw=[(:mode,1),(:padding,1)])
        @test isapprox(unpool(kx;mode=1,padding=1), unpool(ax;mode=1,padding=1))
        @test gradcheck(unpool, kx; kw=[(:mode,1),(:padding,1)])
        @test isapprox(conv4(kw,kx;mode=1,padding=1), conv4(aw,ax;mode=1,padding=1))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:mode,1),(:padding,1)])
        @test isapprox(deconv4(kd,kx;mode=1,padding=1), deconv4(ad,ax;mode=1,padding=1))
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:mode,1),(:padding,1)])

        ### mode=2 (only for pool)
        # @test isapprox(pool(kx;mode=2,padding=1), pool(ax;mode=2,padding=1)) ## mode=2 is not supported in NNlib #218.
        @test gradcheck(pool, kx; kw=[(:mode,2),(:padding,1)])
        # @test isapprox(unpool(kx;mode=2,padding=1), unpool(ax;mode=2,padding=1))  ## mode=2 is not supported in NNlib #218.
        @test gradcheck(unpool, kx; kw=[(:mode,2),(:padding,1)])

        ### alpha=2 (default=1)
        @test isapprox(pool(kx;alpha=2), pool(ax;alpha=2))
        @test gradcheck(pool, kx; kw=[(:alpha,2)])
        @test isapprox(unpool(kx;alpha=2), unpool(ax;alpha=2))
        @test gradcheck(unpool, kx; kw=[(:alpha,2)])
        @test isapprox(pool(kx;alpha=2,mode=1,padding=1), pool(ax;alpha=2,mode=1,padding=1))
        @test gradcheck(pool, kx; kw=[(:alpha,2),(:mode,1),(:padding,1)])
        @test isapprox(unpool(kx;alpha=2,mode=1,padding=1), unpool(ax;alpha=2,mode=1,padding=1))
        @test gradcheck(unpool, kx; kw=[(:alpha,2),(:mode,1),(:padding,1)])

        # @test isapprox(pool(kx;alpha=2,mode=2,padding=1), pool(ax;alpha=2,mode=2,padding=1)) ## mode=2 is not supported in NNlib #218.
        @test gradcheck(pool, kx; kw=[(:alpha,2),(:mode,2),(:padding,1)])
        # @test isapprox(unpool(kx;alpha=2,mode=2,padding=1), unpool(ax;alpha=2,mode=2,padding=1)) ## mode=2 is not supported in NNlib #218.
        @test gradcheck(unpool, kx; kw=[(:alpha,2),(:mode,2),(:padding,1)])

        @test isapprox(conv4(kw,kx;alpha=2), conv4(aw,ax;alpha=2))
        @test gradcheck(conv41, (kw,kx); rtol=TOL, kw=[(:alpha,2)])
        @test isapprox(deconv4(kd,kx;alpha=2), deconv4(ad,ax;alpha=2))
        @test gradcheck(deconv41, (kd,kx); rtol=TOL, kw=[(:alpha,2)])

        # 370-3
        m = M370(param(5,5,1,1))
        path = tempname()*".jld2"
        Knet.save(path,"m",m)
        gpusave = gpu()
        gpu(-1)
        mcpu = Knet.load(path,"m")
        @test conv4(mcpu.layer,randn(Float32,20,20,1,1)) isa Array
        gpu(gpusave)
        rm(path)
    end
    end
end

nothing


text:
# 2D FWI on Overthrust model using minConf library
# Authors:
# Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
# Mathias Louboutin, mlouboutin3@gatech.edu
# Date: January 2022

using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO, PyPlot
using SetIntersectionProjection

# Load starting model
n,d,o,m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ m0)

# Load data
block = segy_read("$(JUDI.JUDI_DATA)/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0)  # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
niterations = parse(Int, get(ENV, "NITER", "10"))
batchsize = 10
fhistory_SGD = zeros(Float32,niterations)

########## Setup constraints
# with constraints:
options=PARSDMM_options()
options.FL=Float32
options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma = true
options.adjust_rho = true
options.adjust_feasibility_rho = true
options.Blas_active = true
options.maxit = 1000
options.feas_tol = 0.001
options.obj_tol = 0.001
options.evol_rel_tol = 0.00001

options.rho_ini=[1.0f0]

set_zero_subnormals(true)
BLAS.set_num_threads(2)
options.parallel=false
options.feasibility_only = false
options.zero_ini_guess=true

constraint = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
vmin = ones(Float32,model0.n) .* 1.3f0
vmax = ones(Float32,model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]

# Slowness squared [s^2/km^2]
m_min = vec((1f0 ./ vmax).^2)
m_max = vec((1f0 ./ vmin).^2)
set_type = "bounds"
TD_OP = "identity"
app_mode = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#TV
(TV,dummy1,dummy2,dummy3) = get_TD_operator(model0,"TV",options.FL)
m_min = 0.0
m_max = norm(TV*vec(v0),1) *2.0f0
set_type = "l1"
TD_OP = "TV"
app_mode = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#set up constraints, precompute some things and define projector
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,model0,options.FL)
(TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,set_Prop,model0,options)
options.rho_ini = ones(length(TD_OP))*10.0

proj_intersection = x-> PARSDMM(x, AtA, TD_OP, set_Prop, P_sub, model0, options)

function proj(input)
    input = Float32.(input)
    (x,dummy1,dummy2,dymmy3) = proj_intersection(vec(input.data))
    return reshape(x, model0.n)
end

########## Run
F0 = judiModeling(deepcopy(model0), q.geometry, d_obs.geometry)
ls = BackTracking(order=3, iterations=10)

# Main loop
for j=1:niterations

    # get fwi objective function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0,q[i],d_obs[i])
    p = -gradient/norm(gradient, Inf)
    println("FWI iteration no: ",j,"; function value: ",fval)
    fhistory_SGD[j] = fval

    # linesearch
    function ϕ(α)
        F0.model.m .= proj(model0.m .+ α * p)
        misfit = .5*norm(F0[i]*q[i] - d_obs[i])^2
        return misfit
    end
    step, fval = ls(ϕ, 1f0, fval, dot(gradient, p))

    # Update model and bound projection
    model0.m .= proj(model0.m .+ step .* p)
end

figure(); imshow(sqrt.(1f0./adjoint(model0.m))); title("FWI with SPG")


text:
using ZChop
function squeezing_param(N, ΔJ1, J2m, J3m)
    """
        ξ2 = squeezing_param(N, ΔJ1, J2m, J3m)

    Returns the squeezing parameter defined, e.g., as
    the inverse of Eq. (1) in Phys. Rev. A 65, 061801 (2002).

    The state is squeezed if greater than one.
    """
    return (J2m.^2 + J3m.^2) ./ (N * ΔJ1)
end

function squeezing_param(rho::QContinuousMeasurement.FixedjState)
    """
        ξ2 = squeezing_param(rho)

    Returns the squeezing parameter defined, e.g., as
    the inverse of Eq. (1) in Phys. Rev. A 65, 061801 (2002).

    The state is squeezed if greater than one.
    """
    Nj = nspins(rho)
    firstblock = block_sizes(Nj)[1]
    (Jx, Jy, Jz) =  map(x -> sparse(x[1:firstblock, 1:firstblock]), jspin(Nj))

    Jy2 = Jy^2
    Jzm = expectation_value!(rho, Jz)
    Jxm =  expectation_value!(rho, Jx)
    Jym =  expectation_value!(rho, Jy)
    Jy2m = expectation_value!(rho, Jy2)
    return (Jzm.^2 + Jxm.^2) ./ (Nj * (Jy2m - Jym.^2))
end

function squeezing_param(rho::QContinuousMeasurement.State)
    """
        ξ2 = squeezing_param(rho)

    Returns the squeezing parameter defined, e.g., as
    the inverse of Eq. (1) in Phys. Rev. A 65, 061801 (2002).

    The state is squeezed if greater than one.
    """
    Nj = nspins(rho)
    (Jx, Jy, Jz) = jspin(Nj)
    Jy2 = Jy^2
    Jzm = expectation_value!(rho, Jz)
    Jxm =  expectation_value!(rho, Jx)
    Jym =  expectation_value!(rho, Jy)
    Jy2m = expectation_value!(rho, Jy2)
    return (Jzm.^2 + Jxm.^2) ./ (Nj * (Jy2m - Jym.^2))
end

"""
    density(M)

Returns the density of a sparse matrix
"""
function density(s::SparseMatrixCSC)
    return length(s.nzval) / (s.n * s.m)
end

function Unconditional_QFI_Dicke(Nj::Int64, Tfinal::Real, dt::Real;
    κ::Real=1.,                    # Independent noise strength
    κcoll::Real=1.,                # Collective noise strength
    ω::Real=0.0                   # Frequency of the Hamiltonian
    )
    return Eff_QFI_HD_Dicke(Nj, 1, Tfinal, dt; κ=κ, κcoll=κcoll, ω=ω, η=0.0)
end

# Specialize zchop! to BlockDiagonal
function ZChop.zchop!(A::BlockDiagonal)
    for b in blocks(A)
        zchop!(b)
    end
end


text:
#RSDeltaSigma: mapCtoD (State space utility)
#-------------------------------------------------------------------------------

#==Helper functions
===============================================================================#
function B2formula(Ac, t1, t2, B2)
	local term=0
	local term1=0
	if t1==0 && t2==0
		term = B2
		return term
	end
	n = size(Ac, 1)
	tmp = eye(n) - exp(-Ac)
	if cond(tmp) < 1e6
		term = ( (exp(-Ac*t1) - exp(-Ac*t2)) * inv(tmp)) * B2
		return term
	end
	#Numerical trouble. Perturb slightly and check the result
	ntry = 0
	k = sqrt(eps(1.0))
	Ac0 = Ac
	#Make this operation repeatable
	#rng_state = rng; rng('default') #Original code
	Random.seed!(0)
	while ntry <2
		Ac = Ac0 + k*rand(n,n)
		tmp = eye(n) - exp(-Ac)
		if cond(tmp) < 1/sqrt(eps(1.0))
			ntry = ntry+1
			if ntry == 1
				term = ( (exp(-Ac*t1) - exp(-Ac*t2)) * inv(tmp)) * B2
			else
				term1 = ( (exp(-Ac*t1) - exp(-Ac*t2)) * inv(tmp)) * B2
			end
		end
		k = k*sqrt(2)
	end
	Random.seed!() #original code restores state: rng(rng_state)
	if norm(term1-term) > 1e-3
		@warn("Inaccurate calculation in mapCtoD.")
	end
	return term
end

"""`(sys, Gp) = mapCtoD(sys_c, t=[0 1], f0=0, calcGp=false)`

Map a MIMO continuous-time system to a SIMO discrete-time equivalent.
The criterion for equivalence is that the sampled pulse response
of the CT system must be identical to the impulse response of the DT system.
i.e. If yc is the output of the CT system with an input vc taken
from a set of DACs fed with a single DT input v, then y, the output
of the equivalent DT system with input v satisfies
y(n) = yc(n-) for integer n. The DACs are characterized by
rectangular impulse responses with edge times specified in the t matrix.

# Input
 - `sys_c`: The LTI description of the CT system.
 - `t`: The edge times of the DAC pulse used to make CT waveforms
   from DT inputs. Each row corresponds to one of the system
   inputs; `[-1 -1]` denotes a CT input. The default is `[0 1]`,
   for all inputs except the first.
 - `f0`: The (normalized) frequency at which the Gp filters' gains are
   to be set to unity. Default 0 (DC).

# Output
 - `sys`: The LTI description for the DT equivalent
 - `Gp`: The mixed CT/DT prefilters which form the samples
   fed to each state for the CT inputs.

# Reference
 - R. Schreier and B. Zhang, "Delta-sigma modulators employing continuous-time
   circuitry," IEEE Transactions on Circuits and Systems I, vol. 43, no. 4,
   pp. 324-332, April 1996.
"""
function mapCtoD(sys_c; t=[0 1], f0::Float64=0.0, calcGp::Bool = false)
	ni = size(sys_c.B, 2)

	if size(t) == (1,2) && (ni > 1) #VERIFYME
		t = [-1 -1; ones(ni-1,1)*t]
	end
	if size(t) != (ni, 2) #VERIFYME
		error("The t argument has the wrong dimensions.")
	end

	di = ones(ni)
	for i in 1:ni
		if t[[i],:]==[-1 -1]
			di[i] = 0
		end
	end
	di = (di .!= 0) #Convert to logical array.
	i!di = findall(.!di); idi = findall(di)

	#c2d assumes t1=0, t2=1.
	#Also c2d often complains about poor scaling and can even produce
	#incorrect results.
	NEEDS_DELAY = false #VERSIONSWITCH
	if NEEDS_DELAY
		#set(sys_c,'InputDelay',0)
	end
	sys = _c2d(sys_c, 1, :zoh)
	Ac = sys_c.A
	Bc1 = sys_c.B[:, i!di]
	A = sys.A; B = sys.B; C = sys.C; D = sys.D
	B1 = B[:, i!di]; D1 = D[:, i!di]

	#Examine the discrete-time inputs to see how big the
	#augmented matrices need to be.
	n = size(A,1)
	t2 = ceil.(Int, t[idi,2])
	esn = (t2 .== t[idi,2]) .& (D[[1],idi] .!= 0)' #extra states needed?
	np = n + maximum(t2 .- 1 + esn)

	#Augment A to np x np, B to np x 1, C to 1 x np.
	Ap = padb( padr(A,np), np )
	for i in n+2:np
		Ap[i, i-1] = 1
	end
	Bp = zeros(np, 1)
	if np>n
		Bp[n+1] = 1
	end
	Cp = padr(C, np)
	Dp = 0

	#Add in the contributions from each DAC
	for i in findall(di)
		t1 = t[i,1]
		t2 = t[i,2]
		B2 = B[:,i]
		D2 = D[:,i]
		if t1==0 && t2==1 && all(D2 .== 0) #No fancy stuff necessary
			Bp = Bp + padb(B2,np)
		else
			n1 = floor(Int, t1); n2 = ceil(Int, t2)-n1-1
			t1 = t1-n1;	t2 = t2-n2-n1
			if t2==1 && all(D2 .!= 0)
				n2 = n2+1
				extraStateNeeded = true
			else
				extraStateNeeded = false
			end
			nt = n+n1+n2
			if n2>0 && t2!=1
				Ap[1:n,nt] = Ap[1:n,nt] + B2formula(Ac, 0, t2, B2)
			end
			local Btmp
			if n2>0 #pulse extends to the next period
				Btmp = B2formula(Ac,t1,1,B2)
			else    #pulse ends in this period
				Btmp = B2formula(Ac,t1,t2,B2)
			end
			if n1>0
				Ap[1:n,n+n1] = Ap[1:n,n+n1] + Btmp
			else
				Bp = Bp + padb(Btmp,np)
			end
			if n2>0
				Cp = Cp + padr([zeros(size(D2,1),n+n1) D2*ones(1,n2)], np)
			end
		end
	end

	sys = _ss(Ap, Bp, Cp, Dp, 1)
	Gp = nothing
	if any(di .== 0)
		if calcGp
			#Compute the prefilters and add in the CT feed-ins.
			#Gp = inv(sI-Ac) * (zI-A)/z *Bc1
			(n, m) = size(Bc1)
			Gp = fill(CTDTPrefilters[], n,m)
			ztf = Array{Any}(nothing, size(Bc1)) #!!Make this like stf: an array of zpk objects
			#Compute the z-domain portions of the filters
			ABc1 = A*Bc1
			for h in 1:m
				for i in 1:n
					if Bc1[i,h]==0
						ztf[i,h] = _zpk([],0,-ABc1[i,h],1)
					else
						ztf[i,h] = _zpk(ABc1[i,h]/Bc1[i,h],0,Bc1[i,h],1)
					end
				end
			end
			#Compute the s-domain portions of each of the filters
			stf = _zpk(_ss(Ac,eye(n),eye(n),zeros(n,n))) #Doesn't do pole-zero cancelation
			for h in 1:m
				for i in 1:n
					for j=1:n
						if stf.k[i,j]!=0 && ztf[j].k!=0 #A non-zero term
							push!(Gp[i,h], CTDTPrefilters(stf(i,j), ztf[j,h]))
						end
					end
				end
			end
			if f0 != 0 #Need to correct the gain terms calculated by c2d
				#B1 = gains of Gp @f0;
				for h in 1:m
					for i in 1:n
						B1[i,h] = evalMixedTF(Gp[i,h],f0)
						#abs() used because ss() whines if B has complex entries...
						#This is clearly incorrect.
						#I've fudged the complex stuff by including a sign....
						B1[i,h] = abs(B1[i,h]) * sign(real(B1[i,h]))
						if abs(B1[i,h]) < 1e-9
							B1[i,h] = 1e-9 #This prevents NaN in line 174
						end
					end
				end
			end
			#Adjust the gains of the pre-filters
			for h in 1:m
				for i in 1:n
					for j in 1:length(Gp[i,h])
						Gp[i,h][j].Hs.k = Gp[i,h][j].Hs.k/B1[i,h] #This is line 174
					end
				end
			end
			newB = [padb(B1,np) sys.B]; newD = [D1 sys.D]
			sys = _ss(sys.A, newB, sys.C, newD)
		else #Cheat and just dublicate the B1 terms
			newB = [padb(Bc1,np) sys.B]; newD = [D1 sys.D]
			sys = _ss(sys.A, newB, sys.C, newD)
		end
	end

	return (sys, Gp)
end

#Last line


text:
module Qiskit

using PyCall

const qiskit = PyNULL()
const IBMQ = PyNULL()
const Aer = PyNULL()

function save_account(token)
    IBMQ.save_account(token)
end

function __init__()
    copy!(qiskit, pyimport("qiskit"))
    copy!(Aer, qiskit.Aer)
    copy!(IBMQ, qiskit.IBMQ)
end

include("register.jl")

end # module


text:
#= Test Constructors =#

for k in kernel_functions
    @testset "Testing $k" begin
        def_args, alt_args = get(kernel_functions_arguments, k, ((), ()))

        # Test Constructors
        for T in FloatingPointTypes
            for i = 1:length(alt_args)
                K = (k)([typeof(θ) <: AbstractFloat ? convert(T,θ) : θ for θ in alt_args[1:i]]...)
                for j = 1:length(alt_args)
                    @test MLK.getvalue(getfield(K,j)) == (j <= i ? alt_args[j] : def_args[j])
                end
            end
        end

        # Test Conversions
        K = (k)()
        for psi in (Kernel, MercerKernel, NegativeDefiniteKernel)
            if k <: psi
                for T1 in FloatingPointTypes
                    T2 = T1 == Float64 ? Float32 : T1
                    @test eltype(convert(psi{T2},K)) == T2
                end
            end
        end

        # Test Properties
        @test MLK.ismercer(K) == (typeof(K) <: MLK.MercerKernel ? true : false)
        @test MLK.isnegdef(K) == (typeof(K) <: MLK.NegativeDefiniteKernel ? true : false)
        @test MLK.isstationary(K) == MLK.isstationary(MLK.basefunction(K))
        @test MLK.isisotropic(K) == MLK.isisotropic(MLK.basefunction(K))

        # Test Display
        @test eval(Meta.parse(string(K))) == K
        @test show(devnull, K) == nothing
    end
end

@testset "Testing MLK.kappa" begin
    for k in kernel_functions
        k_tmp = get(kernel_functions_kappa, k, x->error(""))
        for T in FloatingPointTypes
            K = convert(Kernel{T}, (k)())
            args = T[MLK.getvalue(getfield(K,theta)) for theta in fieldnames(typeof(K))]

            for z in (zero(T), one(T), convert(T,2))
                v = MLK.kappa(K, z)
                v_tmp = (k_tmp)(z, args...)
                @test isapprox(v, v_tmp)
            end
        end
    end
end

@testset "Testing MLK.gettheta" begin
    for k in kernel_functions
        K = (k)()
        K_theta = MLK.gettheta(K)
        @test K_theta == [MLK.gettheta(getfield(K,field)) for field in MLK.thetafieldnames(K)]
        @test MLK.checktheta(K, K_theta) == true
    end
end

@testset "Testing MLK.settheta!" begin
    for T in FloatingPointTypes
        alpha1 = one(T)
        alpha2 = convert(T,0.6)
        gamma1 = one(T)
        gamma2 = convert(T,0.4)
        K = MLK.GammaExponentialKernel(alpha1,gamma1)
        MLK.settheta!(K, [log(alpha2); log(gamma2)])
        @test getvalue(K.alpha) == alpha2
        @test getvalue(K.gamma) == gamma2

        @test_throws DimensionMismatch MLK.settheta!(K, [one(T)])
        @test_throws DimensionMismatch MLK.settheta!(K, [one(T); one(T); one(T)])
    end
end

@testset "Testing MLK.checktheta" begin
    for T in FloatingPointTypes
        K = MLK.GammaExponentialKernel(one(T),one(T))

        @test MLK.checktheta(K, [log(one(T)); log(one(T))]) == true
        @test MLK.checktheta(K, [log(one(T)); log(convert(T,2))]) == false

        @test_throws DimensionMismatch MLK.checktheta(K, [one(T)])
        @test_throws DimensionMismatch MLK.checktheta(K, [one(T); one(T); one(T)])
    end
end


text:
module HackerNewsSlackApp

export handler, publish, monitor_hackernews

using HTTP
using Sockets
using JSON3
using Dates
using Memento

const ROUTER = HTTP.Router()
# const LOGGER = Memento.config!("debug"; fmt="[{level} | {name}]: {msg}")
info(logger, msg) = begin sleep(0.5); println(now(), " ", msg); end
error(logger, msg) = println(now(), " ", msg)
warn(logger, msg) = println(now(), " ", msg)
debug(logger, msg) = println(now(), " ", msg)
const LOGGER = nothing

# installed web hook
const WEB_HOOK_URL = ENV["HACKER_NEWS_SLACK_WEB_HOOK"]

# what to monitor
struct Monitor
    channel::Channel
    keywords::Any   # change to vector{string}?
    sent::Dict{String,DateTime}
end

# Slash command

function start_slash_command_server(port)
    HTTP.@register(ROUTER, "POST", "/", handler)
    HTTP.serve(ROUTER, IPv4(0), port, verbose = true)
end

function handler(req::HTTP.Request)
    println(req)
    try
        headers = ["Content-type" => "application/json"]
        body = JSON3.write(msg)
        return HTTP.Response(200, headers, body = body)
    catch ex
        error(LOGGER, "exception ex=$ex")
        return HTTP.Response(500, "Sorry, got an exception.")
    end
end

# Web hook

function publish(msg::AbstractString, monitor::Monitor)
    if already_published(msg, monitor)
        debug(LOGGER, "Already published message: $msg")
    else
        info(LOGGER, "Publishing message: $msg")
        headers = ["Content-type" => "application/json"]
        data = JSON3.write(Dict("text" => msg))
        HTTP.post(WEB_HOOK_URL, headers, data)
        monitor.sent[msg] = now()  # remember that we have already sent
    end
    return nothing
end

already_published(msg::AbstractString, monitor::Monitor) = haskey(monitor.sent, msg)

# Monitoring

function create_monitor(f::Function, keywords; interval = 60)
    info(LOGGER, "creating monitor")
    monitor = Monitor(Channel(32), keywords, Dict{String,DateTime}())
    task = @async monitor_loop(f, monitor, interval)
    return (monitor = monitor, task = task)
end

function monitor_loop(f::Function, monitor::Monitor, interval::Int)
    info(LOGGER, "started monitoring")
    while true
        f(monitor.keywords, monitor)  # execute client function
        if isready(monitor.channel)
            instruction = take!(monitor.channel)
            instruction == "stop" && break
        end
        info(LOGGER, "sleeping $interval seconds")
        sleep(interval)
    end
    info(LOGGER, "stopped monitoring")
end

function stop_monitor(monitor::Monitor)
    put!(monitor.channel, "stop")
    return nothing
end

# Hacker News API

# Container for a story. Underlying data is a Dict from parsing JSON.
struct Story
    by::String
    descendants::Union{Nothing,Int}
    score::Int
    time::Int
    id::Int
    title::String
    kids::Union{Nothing,Vector{Int}}
    url::Union{Nothing,String}
end

# Construct a Story from a Dict (or Dict-compatible) object
function Story(obj)
    value = (x) -> get(obj, x, nothing)
    return Story(
        obj[:by],
        value(:descendants),
        obj[:score],
        obj[:time],
        obj[:id],
        obj[:title],
        value(:kids),
        value(:url))
end

# Find top 500 stories from HackerNews.  Returns an array of story id's.
function fetch_top_stories()
    url = "https://hacker-news.firebaseio.com/v0/topstories.json"
    response = HTTP.request("GET", url)
    return JSON3.read(response.body)
end

# Get a specific story item
function fetch_story(id)
    url = "https://hacker-news.firebaseio.com/v0/item/$(id).json"
    response = HTTP.request("GET", url)
    return Story(JSON3.read(response.body))
end

# Get title for top N stories
top_stories(n::Int) = let stories = fetch_top_stories()
    fetch_story.(stories[1:min(n,end)])
end

# Top story titles
top_titles(n::Int) = getfield.(top_stories(n), :title)

# Interesting title?
is_interesting(title::AbstractString, keywords::AbstractVector{T}) where T <: AbstractString =
    match(Regex("\\b(" * join(keywords, "|") * ")\\b"), lowercase(title)) !== nothing

# Check if any of the top stories contains certain keywords
function monitor_hackernews(keywords::AbstractVector{T}; interval = 300) where T <: AbstractString
    return create_monitor(keywords; interval = interval) do kw, monitor
        info(LOGGER, "finding top 10")
        tt = top_titles(10)
        ii = is_interesting.(tt, Ref(kw))
        if any(ii)
            publish.(tt[ii], Ref(monitor))
        end
        info(LOGGER, "done")
    end
end

end # module


text:
using TimedSound
using Lazy: @>
using Base.Test

x = leftright(ramp(tone(1kHz,1s)),ramp(tone(1kHz,1s)))
rng() = MersenneTwister(1983)

show_str = "1.0 s 16 bit PCM stereo sound
Sampled at 44100 Hz
▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆
▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆"

@testset "Sound Indexing" begin
  @test isapprox(duration(x[0s .. 0.5s,:]),0.5s; atol = 1/samplerate(x))
  @test isapprox(duration(x[0.5s .. ends,:]),0.5s; atol = 1/samplerate(x))

  @test x[:,:left][0s .. 0.5s] == x[0s .. 0.5s,:left]
  @test x[:,:left] == x[:,:right]
  @test x[0.5s .. 0.75s] == x[0.5s .. 0.75s,:]
  @test x[0.5s .. ends] == x[0.5s .. ends,:]
  mylen = nsamples(x[0.1s .. 0.6s])
  newx = copy(x)
  @test (newx[0.1s .. 0.6s] = x[50:(50+mylen-1),:]) == x[50:(50+mylen-1),:]
  @test (newx[0.1s .. 0.6s,:] = x[50:(50+mylen-1),:]) == x[50:(50+mylen-1),:]
  myend = (0.1s + duration(x[0.5s .. ends]))
  @test (newx[0.5s .. ends] = x[0.1s .. myend]) == x[0.1s .. myend]
  @test (newx[0.5s .. ends,:] = x[0.1s .. myend,:]) == x[0.1s .. myend,:]
  @test (newx[0.1s .. 0.6s,:] = x[0.2s .. 0.700005s,:]) == x[0.2s .. 0.700005s,:]

  @test x[0s .. 22050samples,:] == x[0s .. 0.5s,:]
  @test x[22050samples .. 44100samples,:] == x[0.5s .. 1s,:]
  @test x[22050samples .. 1s,:] == x[0.5s .. 1s,:]
  @test x[22050samples .. ends,:] == x[0.5s .. ends,:]
  @test x[0.5s .. 1s,:] == (x[0s .. 22051samples,:] = x[0.5s .. 1s,:])
  @test x[0.5s .. ends,:] == (x[22050samples .. ends,:] = x[0.5s .. ends,:])

  @test x[22050samples .. 44100samples] == x[0.5s .. 1s]
  @test x[22050samples .. 1s] == x[0.5s .. 1s]
  @test x[22050samples .. ends] == x[0.5s .. ends]
  @test x[0.5s .. 1s] == (x[0s .. 22051samples] = x[0.5s .. 1s])
  @test x[0.5s .. ends] == (x[22050samples .. ends] = x[0.5s .. ends])

  @test_throws BoundsError x[0.5s .. 2s,:]
  @test_throws BoundsError x[-0.5s .. 0.5s,:]
  @test_throws BoundsError x[-0.5s .. ends,:]
  @test_throws BoundsError x[:,:doodle]
  @test_throws BoundsError x[0.5s .. 2s]
  @test_throws BoundsError x[-0.5s .. 0.5s]
  @test_throws BoundsError x[-0.5s .. ends]
  @test_throws BoundsError x[0.5s .. 2s,:] = 1:10
  @test_throws BoundsError x[-0.5s .. 0.5s,:] = 1:10
  @test_throws BoundsError x[-0.5s .. ends,:] = 1:10
  @test_throws BoundsError x[:,:doodle] = 1:10
  @test_throws BoundsError x[0.5s .. 2s] = 1:10
  @test_throws BoundsError x[-0.5s .. 0.5s] = 1:10
  @test_throws BoundsError x[-0.5s .. ends] = 1:10
  @test_throws DimensionMismatch x[0.5s .. 0.6s,:] = 1:10
  @test_throws DimensionMismatch x[0.2s .. 0.5s,:] = 1:10
  @test_throws DimensionMismatch x[0.2s .. ends,:] = 1:10
  @test_throws DimensionMismatch x[0.5s .. 0.6s] = 1:10
  @test_throws DimensionMismatch x[0.5s .. 0.6s] = 1:10
  @test_throws DimensionMismatch x[0.5s .. ends] = 1:10

  @test nsamples(x[0s .. 0.5s,:]) + nsamples(x[0.5s .. ends,:]) == nsamples(x)
  strbuff = IOBuffer()
  show(strbuff,playable(x))
  @test show_str == String(strbuff)
end

@inline same(x,y) = isapprox(x,y,rtol=1e-6)

@testset "Sound Construction" begin
  @test same(sound("sounds/tone.wav"),x)
  @test same(sound("sounds/two_tone.wav"),
             [tone(1kHz,100ms);silence(800ms);tone(1kHz,100ms)])
  a,b = tone(1kHz,200ms),tone(2kHz,200ms)
  @test [sound(zeros(10)); sound(zeros(Float32,10))[:,:]] ==
           [sound(zeros(10,1)); sound(zeros(10,1))]
  @test_throws ErrorException [sound(zeros(10)); sound(zeros(10);sample_rate=22050Hz)]
  @test leftright(a,b)[:,:left] == a
  @test leftright(a,b)[:,:right] == b
  @test size([tone(1kHz,0.2s); tone(2kHz,0.2s)],2) == 1
  @test [tone(1kHz,0.2s); leftright(tone(1.5kHz,0.2s),tone(0.5kHz,0.2s))] ==
    [leftright(tone(1kHz,0.2s),   tone(1kHz,0.2s));
     leftright(tone(1.5kHz,0.2s), tone(0.5kHz,0.2s))]
  @test_throws ErrorException ramp(tone(1kHz,50ms),100ms)
  @test same(sound("sounds/rampon.wav"),rampon(tone(1kHz,1s)))
  @test same(sound("sounds/rampoff.wav"),rampoff(tone(1kHz,1s)))
  @test same(sound("sounds/fadeto.wav"),
             fadeto(tone(1kHz,0.5s),tone(2kHz,0.5s)))
  @test fadeto(leftright(tone(1kHz,100ms),tone(1.5kHz,100ms)),
               tone(2kHz,100ms)) != 0
  @test same(sound("sounds/noise.wav"),noise(1s,rng=rng()))
  @test same(sound("sounds/bandpass.wav"),
             @> noise(1s,rng=rng()) bandpass(400Hz,800Hz))
  @test same(sound("sounds/complex.wav"),
             @> harmonic_complex(200Hz,0:5,ones(6),1s) attenuate(20))
end


text:
d = Deterministic(1)

@test rand(d) == 1
@test rand(MersenneTwister(4), d) == 1
@test collect(support(d)) == [1]
@test Random.gentype(d) == typeof(1)
@test Random.gentype(typeof(d)) == typeof(1)
@test pdf(d, 0) == 0.0
@test pdf(d, 1) == 1.0
@test mode(d) == 1
@test mean(d) == 1
@test typeof(mean(d)) == typeof(mean([1]))

d2 = Deterministic(:symbol)
@test rand(d2) == :symbol
@test rand(MersenneTwister(4), d2) == :symbol
@test collect(support(d2)) == [:symbol]
@test Random.gentype(d2) == typeof(:symbol)
@test Random.gentype(typeof(d2)) == typeof(:symbol)
@test pdf(d2, :another) == 0.0
@test pdf(d2, :symbol) == 1.0
@test mode(d2) == :symbol
@test mean(d2) == :symbol
@test typeof(mean(d2)) == Symbol


text:
# ## Font

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Measures, LaTeXStrings

plot(sin, 0, 2π,
    title=L"y = \sin(x)",
    titlefont=font(40), ## title

    xlabel=L"x",
    ylabel="y",
    xguidefont=font(30), ## x-label
    yguidefont=font(20), ## x-label
    ## guidefont=font(20), # x,y-label

    xtick=(0:0.5:2π, ["\$ $(i) \$" for i in 0:0.5:2π]),
    ytick=-1:0.5:1,
    xtickfont=font(15),
    ytickfont=font(20),
    ## tickfont=font(10), # x,y-ticks font

    label="Sin function",
    legendfont=font(12),

    xlims=(0,2π),
    ylims=(-1,1),
    bottom_margin=5mm,
    left_margin=10mm,
    top_margin=15mm
)


text:
function calc_seat_id(pass)
    digits = map(c -> c == 'F' || c == 'L' ? '0' : '1', pass)
    row = parse(Int64, digits[1:7]; base = 2)
    column = parse(Int64, digits[8:end]; base = 2)
    row * 8  + column
end

function seat_ids_on_list(boarding_passes)
    ids = Set{Int64}()
    for pass in boarding_passes
        push!(ids, calc_seat_id(pass))
    end
    ids
end

function find_missing_id(boarding_passes)
    existing_ids = seat_ids_on_list(boarding_passes)
    for row in 1:126, column in 0:7
        seat_id = row * 8 + column
        if seat_id ∉ existing_ids && (seat_id - 1) ∈ existing_ids && (seat_id + 1) ∈ existing_ids
            return seat_id
        end
    end
end

function main()
    lines = readlines("../resources/input_05.txt")
    result = find_missing_id(lines)
    println(result)
end

main()


text:
module ConstructionBaseExtras

using ConstructionBase, Requires

function __init__()
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
        include("staticarrays.jl")
    end
end

end


text:
## TOP LEVEL CALCULATION METHODS - KEEPING TRACK OF CALCULATION STATE ##
function update_pricing_engine(inst::Instrument, pe::PricingEngine)
  T = get_pricing_engine_type(inst)
  if typeof(pe) <: T
    inst.pricingEngine = pe
    inst.lazyMixin.calculated = false
  else
    # we have to clone
    inst = clone(inst, pe)
  end

  return inst
end

function npv(inst::Instrument)
  calculate!(inst)

  return inst.results.value
end

function perform_calculations!(inst::Instrument)
  reset!(inst.results)
  _calculate!(inst.pricingEngine, inst)

  return inst
end

## MISC TYPES ##
struct LongPosition <: PositionType end
struct ShortPosition <: PositionType end


text: @Base.deprecate predict(mm::LinearModel, newx::AbstractMatrix, interval::Symbol, level::Real = 0.95) predict(mm, newx; interval=interval, level=level)


text:
module FourierTransforms

using Compat: cispi
using Primes

################################################################################

function direct_ft!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    N = length(x)
    @assert length(X) == N
    RT = typeof(real(zero(T)))
    for k in 1:N
        s = zero(T)
        for n in 1:N
            s += cispi(-2 * ((n - 1) * (k - 1) % N) / RT(N)) * x[n]
        end
        X[k] = s
    end
    return X
end

direct_ft(x::AbstractVector) = direct_ft!(similar(x), x)

################################################################################

function direct_ft_1!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    @assert length(X) == length(x) == 1
    X[1] = x[1]
    return X
end

direct_ft_1(x::AbstractVector) = direct_ft_1!(similar(x), x)

################################################################################

function direct_ft_2!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    @assert length(X) == length(x) == 2
    X[1] = x[1] + x[2]
    X[2] = x[1] - x[2]
    return X
end

direct_ft_2(x::AbstractVector) = direct_ft_2!(similar(x), x)

################################################################################

function direct_ft_4!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    @assert length(X) == length(x) == 4
    X[1] = (x[1] + x[3]) + (x[2] + x[4])
    X[2] = (x[1] - x[3]) - im * (x[2] - x[4])
    X[3] = (x[1] + x[3]) - (x[2] + x[4])
    X[4] = (x[1] - x[3]) + im * (x[2] - x[4])
    return X
end

direct_ft_4(x::AbstractVector) = direct_ft_4!(similar(x), x)

################################################################################

# <https://en.wikipedia.org/w/index.php?title=Cooley%E2%80%93Tukey_FFTur_algorithm&oldid=1056450535>
function ditfft2!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    N = length(x)
    @assert length(X) == N
    RT = typeof(real(zero(T)))
    if N == 0
        # do nothing
    elseif N == 1
        X[1] = x[1]
    else
        # Ensure N is even
        @assert N % 2 == 0
        N2 = N ÷ 2
        ditfft2!((@view X[(0 * N2 + 1):(1 * N2)]), (@view x[1:2:end]))
        ditfft2!((@view X[(1 * N2 + 1):(2 * N2)]), (@view x[2:2:end]))
        for k in 1:N2
            ϕ = cispi(-2 * (k - 1) / RT(N))
            p = X[0 * N2 + k]
            q = ϕ * X[1 * N2 + k]
            X[0 * N2 + k] = p + q
            X[1 * N2 + k] = p - q
        end
    end
    return X
end

ditfft2(x::AbstractVector) = ditfft2!(similar(x), x)

################################################################################

function choose_radix_sqrt(n::Integer)
    @assert n ≥ 0
    n ≤ 1 && return n
    # Handle special cases efficiently
    if ispow2(n)
        # Power of 2
        return 2
    end
    # Power of 2 times a factor
    p2 = trailing_zeros(n)
    factor = n >> p2
    if factor ∈ (3, 5, 7, 11, 13, 17, 19)
        return factor
    end
    # Use a greedy algorithm to find the largest factor not larger than sqrt(n)
    prime_factor_counts = factor(n)
    prime_factors = sort!(collect(keys(prime_factor_counts)))
    sqrt_n = isqrt(n)
    radix = 1
    for prime_factor in Iterators.reverse(Iterators.filter(≤(sqrt_n), prime_factors))
        for count in 1:prime_factor_counts[prime_factor]
            if radix * prime_factor ≤ sqrt_n
                radix *= prime_factor
            else
                break
            end
        end
    end
    # Always return a radix larger than 1
    radix == 1 && return n
    return radix
end

# X: output
# Y: workspace
# x: input (will be destroyed)
function radix_fft!(X::AbstractVector{T}, Y::AbstractVector{T}, x::AbstractVector{T};
                    choose_radix=choose_radix_sqrt) where {T<:Complex}
    N = length(x)
    @assert length(X) == length(Y) == N
    RT = typeof(real(zero(T)))
    radix = choose_radix(N)
    if N == 0
        # do nothing
    elseif N ≤ radix
        direct_ft!(X, x)
    else
        N₁ = radix              # aka "decimation in time"
        @assert N % N₁ == 0
        N₂ = N ÷ N₁
        x2 = reshape(x, (N₁, N₂))
        X2 = reshape(X, (N₂, N₁))
        Y2 = reshape(Y, (N₂, N₁))
        Z2 = reshape(x, (N₂, N₁))
        for n₁ in 1:N₁
            radix_fft!((@view Y2[:, n₁]), (@view X2[:, n₁]), (@view x2[n₁, :]); choose_radix=choose_radix)
            for n₂ in 1:N₂
                Y2[n₂, n₁] *= cispi(-2 * (n₁ - 1) * (n₂ - 1) / RT(N))
            end
        end
        for n₂ in 1:N₂
            radix_fft!((@view X2[n₂, :]), (@view Z2[n₂, :]), (@view Y2[n₂, :]); choose_radix=choose_radix)
        end
    end
    return X
end

"x will contain output"
radix_fft!(x::AbstractVector; kws...) = radix_fft!(x, similar(x), copy(x); kws...)
"X will contain output, x will be preserved"
radix_fft!(X::AbstractVector, x::AbstractVector; kws...) = radix_fft!(X, similar(x), copy(x); kws...)
"x will be preserved"
radix_fft(x::AbstractVector; kws...) = radix_fft!(similar(x), x; kws...)

################################################################################

export fft!, fft, inv_fft!, inv_fft

"""
    fft!(X::AbstractVector, x::AbstractVector)

Calculate the Fourier transform of `x` and store it into `X`. `x` is
not modified. Both vectors must have the same length and must have
complex element types.

See also: [`fft`](@ref), [`inv_fft!`](@ref).
"""
fft!(X::AbstractVector, x::AbstractVector) = radix_fft!(X, x)

"""
    inv_fft!(x::AbstractVector, X::AbstractVector)

Calculate the inverse Fourier transform of `X` and store it into `x`.
`X` is not modified. Both vectors must have the same length and must
have complex element types.

See also: [`inv_fft`](@ref), [`fft!`](@ref).
"""
function inv_fft!(X::AbstractVector{T}, x::AbstractVector{T}) where {T<:Complex}
    # TODO: Don't modify `x`
    x .= conj(x)
    fft!(X, x)
    x .= conj(x)
    X .= conj(X) / length(X)
    return X
end

"""
    X = fft(x::AbstractVector)
    X::AbstractVector

Calculate the Fourier transform of `x` and return it in a newly
allocate vector. `x` is not modified. The element type of `x` must be
a complex number type.

See also: [`inv_fft`](@ref), [`fft!`](@ref).
"""
fft(x::AbstractVector) = fft!(similar(x), x)

"""
    x = inv_fft(X::AbstractVector)
    x::AbstractVector

Calculate the inverse Fourier transform of `X` and return it in a
newly allocate vector. `X` is not modified. The element type of `X`
must be a complex number type.

See also: [`fft`](@ref), [`inv_fft!`](@ref).
"""
inv_fft(x::AbstractVector) = inv_fft!(similar(x), x)

end


text:
"""
Wrap a single pnml net.

$(TYPEDEF)
$(TYPEDFIELDS)

# Details

"""
struct HLPetriNet{PNTD} <: PetriNet{PNTD}
    net::PnmlNet{PNTD}
end
"Construct from string of valid pnml XML using the first network"
HLPetriNet(str::AbstractString) = HLPetriNet(parse_str(str))
HLPetriNet(model::PnmlModel)  = HLPetriNet(first_net(model))

#-------------------------------------------------------------------------------
# Implement PNML Petri Net interface.
#-------------------------------------------------------------------------------

# Delegate to wrapped net.
pid(hlpn::HLPetriNet) = pid(hlpn.net)

# Flattened to page[1], so simple vectors.
places(hlpn::HLPetriNet)      = places(firstpage(hlpn.net))
transitions(hlpn::HLPetriNet) = transitions(firstpage(hlpn.net))
arcs(hlpn::HLPetriNet)        = arcs(firstpage(hlpn.net))
refplaces(hlpn::HLPetriNet)   = refPlaces(firstpage(hlpn.net))
reftransitions(hlpn::HLPetriNet) = refPlaces(firstpage(hlpn.net))


text:
module Experimental

using LightGraphs
using LightGraphs.SimpleGraphs

export description,
    #isomorphism
    VF2, vf2, IsomorphismProblem, SubgraphIsomorphismProblem, InducedSubgraphIsomorphismProblem,
    could_have_isomorph, has_isomorph, all_isomorph, count_isomorph,
    has_induced_subgraphisomorph, count_induced_subgraphisomorph, all_induced_subgraphisomorph,
    has_subgraphisomorph, count_subgraphisomorph, all_subgraphisomorph,

    ShortestPaths
description() = "This module contains experimental graph functions."

include("isomorphism.jl")
include("vf2.jl") # Julian implementation of VF2 algorithm
include("ShortestPaths/ShortestPaths.jl")

end


text:
# Here we use implicit resolution of symbols from FinEtools:
using FinEtools

# Alternatively one could explicitly import needed functionality:
# import FinEtools.MeshQuadrilateralModule: Q4quadrilateral
# import FinEtools.MeshModificationModule: mergemeshes
# import FinEtools.FESetModule: cat
# import FinEtools.FENodeSetModule: count

function onto_sphere!(fens, radius)
    # Project nodes onto a sphere of given radius.
    for j = 1:count(fens)
        fens.xyz[j,:] = fens.xyz[j,:]*radius/norm(fens.xyz[j,:]);
    end
    return fens;
end

"""
    Q4spheren(radius, nperradius)

Create a surface mesh of 1/8 of a sphere with a given number of elements per circumference.
"""
function Q4spheren(radius, nperradius)
    if (nperradius % 2 != 0)
        nperradius = nperradius+1;
    end
    nL=Int(nperradius/2); nW=Int(nperradius/2);
    a=sqrt(2.0)/2;
    b=1/sqrt(3.0);
    c=0.6*a;
    d=0.6*b;
    xyz=[1 0 0; 0 1 0; 0 0 1; a a 0; 0 a a; a 0 a; b b b];
    conn=[1 4 7 6; 4 2 5 7; 3 6 7 5];
    fens, fes = Q4quadrilateral(xyz[conn[1,:],:], nL, nW);
    fens1, fes1 = Q4quadrilateral(xyz[conn[2,:],:], nL, nW);
    fens, fes1, fes2 = mergemeshes(fens1, fes1, fens, fes, norm(xyz)/1000)
    fes = cat(fes1,fes2);
    fens1, fes1 = Q4quadrilateral(xyz[conn[3,:],:], nL, nW);
    fens, fes1, fes2 = mergemeshes(fens1, fes1, fens, fes, norm(xyz)/1000)
    fes = cat(fes1,fes2);
    fens= onto_sphere!(fens, radius)
    return fens, fes
end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_issues_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_forks_repo_name: PetrKryslUCSD/FinEtoolsUseCase


text:
@testset "LinearRegressor" begin
    # Define some linear, noise-free, synthetic data:
    # with intercept = 0.0
    n = 1000
    rng = StableRNG(1234)
    X, y = MLJBase.make_regression(n, 3, noise=0, intercept=false, rng=rng)
    # Train model with intercept on all data
    linear = LinearRegressor()
    yhat, fr = test_regression(linear, X, y)
    # Training error
    @test norm(yhat - y)/sqrt(n) < 1e-12
    # Get the true intercept?
    @test abs(fr.intercept) < 1e-10
    # Check metadata
end

@testset "MultitargetLinearRegressor" begin
    # Define some linear, noise-free, synthetic data:
    # with intercept = 0.0
    n = 1000
    rng = StableRNG(1234)
    X, Y = make_regression2(n, 3, noise=0, intercept=false, rng=rng)
    Y_mat = MLJBase.matrix(Y)
    # Train model on all data
    multi_linear = MultitargetLinearRegressor()
    Yhat, fr = test_regression(multi_linear, X, Y)
    Yhat_mat = MLJBase.matrix(Yhat)
    # Check if the column names is same after predict
    MLJBase.schema(Yhat).names == MLJBase.schema(Y).names
    # Training error
    @test norm(Yhat_mat - Y_mat)/sqrt(n) < 1e-12
    # Get the true intercept?
    @test norm(fr.intercept .- zeros(size(Y_mat, 2))) < 1e-10
end

@testset "RidgeRegressor" begin
    # Define some linear, noise-free, synthetic data:
    # with intercept = 0.0
    n = 1000
    rng = StableRNG(1234)
    X, y = MLJBase.make_regression(n, 3, noise=0, intercept=false, rng=rng)
    # Train model with intercept on all data with no regularization
    # and no standardization of target.
    ridge = RidgeRegressor(lambda=0.0)
    yhat, fr = test_regression(ridge, X, y)
    # Training error
    @test norm(yhat - y)/sqrt(n) < 1e-12
    # Get the true intercept?
    @test abs(fr.intercept) < 1e-10
    # Check metadata
end

@testset "MultitargetRidgeRegressor" begin
    # Define some linear, noise-free, synthetic data:
    # with intercept = 0.0
    n = 1000
    rng = StableRNG(1234)
    X, Y = make_regression2(n, 3, noise=0, intercept=false, rng=rng)
    Y_mat = MLJBase.matrix(Y)
    # Train model with intercept on all data with no regularization
    # and no standardization of target.
    multi_ridge = MultitargetRidgeRegressor(lambda=0.0)
    Yhat, fr = test_regression(multi_ridge, X, Y)
    Yhat_mat = MLJBase.matrix(Yhat)
    # Check if the column names is same after predict
    MLJBase.schema(Yhat).names == MLJBase.schema(Y).names
    # Training error
    @test norm(Yhat_mat - Y_mat)/sqrt(n) < 1e-12
    # Get the true intercept?
    @test norm(fr.intercept .- zeros(size(Y_mat, 2))) < 1e-10
end


text:
export identityActivation

"""
 identity activation A = identity(Y)

 Input:

   Y - array of features

 Optional Input:

   doDerivative - flag for computing derivative, set via varargin
                  Ex: identityActivation(Y,true);

 Output:

  A  - activation
  dA - derivatives
"""
function identityActivation(Y::Array{T},dA,doDerivative::Bool=false) where {T}

if doDerivative
    if isempty(dA)
        dA = ones(T,Y);
    else
        dA .= ones(T,Y);
    end
end

return A,dA #Depricated? A Isnt even declared lol
end


text:
include("../src/graph.jl")

@testset "Graph search algorithms" begin
    # vertice "7" and "8" is not connected
    examplegraph = (
        nv=8,
        edges=[
            (1, 6),
            (3, 5),
            (3, 4),
            (2, 3),
            (1, 2),
            (4, 5),
            (4, 6),
            (1, 3),
            (7, 8)
        ]
    )

    g = GraphAlgo.UndirectedGraph(examplegraph.nv, examplegraph.edges)

    @testset "Test depth first search" begin
        # source vertex is "1"
        result = GraphAlgo.search(g, 1, GraphAlgo.FindingPaths, GraphAlgo.DepthFirstSearch)
        @test length(result.marked) == g.numvertices
        @test sum(result.marked) == 6
        @test result.numvisited == 6
        @test result.edgeto == [nothing, 3, 4, 6, 3, 1, nothing, nothing]
        @test GraphAlgo.haspathto(result, 5) == true
        @test GraphAlgo.haspathto(result, 7) == false
        @test GraphAlgo.pathto(result, 3) == [1, 6, 4, 3]
    end

    @testset "Test breadth first search" begin
        result = GraphAlgo.search(g, 1, GraphAlgo.FindingPaths, GraphAlgo.BreadthFirstSearch)
        @test length(result.marked) == g.numvertices
        @test sum(result.marked) == 6
        @test result.numvisited == 6
        @test result.edgeto == [nothing, 1, 1, 6, 3, 1, nothing, nothing]
        @test GraphAlgo.haspathto(result, 2) == true
        @test GraphAlgo.haspathto(result, 7) == false
        @test GraphAlgo.pathto(result, 3) == [1, 3]
        @test GraphAlgo.pathto(result, 4) == [1, 6, 4]
        @test GraphAlgo.pathto(result, 5) == [1, 3, 5]
    end

    @testset "Connected components" begin
        cc = GraphAlgo.search(g, GraphAlgo.ConnectedComponents)
        @test sum(cc.marked) == g.numvertices
        @test cc.vids[6] == 1
        @test cc.vids[7] == 2
    end

    @testset "Cycle detection" begin
        @test GraphAlgo.search(g, GraphAlgo.Cycle).hascycle == true

        g = GraphAlgo.UndirectedGraph(3, [(1,2), (2,3)])
        @test GraphAlgo.search(g, GraphAlgo.Cycle).hascycle == false

        examplegraph = (
        nv=8,
        edges=[
            (1, 6),
            (3, 5),
            (1, 2),
            (4, 5),
            (4, 6),
            (7, 8)
            ]
        )

        g = GraphAlgo.UndirectedGraph(examplegraph.nv, examplegraph.edges)
        @test GraphAlgo.search(g, GraphAlgo.Cycle).hascycle == false
    end

end

@testset "Directed Graph" begin
    minidg = (
        nv = 6,
        ne = 4,
        edges = [
            (6, 5),
            (5, 4),
            (4, 6),
            (1, 6),
        ]
    )

    minidg2 = (
        nv = 6,
        ne = 4,
        edges = [
            (6, 5),
            (5, 4),
            (1, 6),
        ]
    )

    # tinyDG.txt from *Algorithm 4th* page 569
    tinydg = (
        nv = 13,
        ne = 22,
        edges = [
            ( 5,  3),
            ( 3,  4),
            ( 4,  3),
            ( 7,  1),
            ( 1,  2),
            ( 3,  1),
            (12, 13),
            (13, 10),
            (10, 11),
            (10, 12),
            ( 9, 10),
            (11, 13),
            (12,  5),
            ( 5,  4),
            ( 4,  6),
            ( 8,  9),
            ( 9,  8),
            ( 6,  5),
            ( 1,  6),
            ( 7,  5),
            ( 7, 10),
            ( 8,  7)
        ]
    )

    # DAG examle from p579 in *Algorithm 4th*.
    dag_p579 = (
        nv = 13,
        ne = 15,
        edges = [
            (1, 2),
            (1, 6),
            (1, 7),
            (3, 1),
            (3, 4),
            (4, 6),
            (6, 5),
            (7, 5),
            (7, 10),
            (8, 7),
            (9, 8),
            (10, 11),
            (10, 12),
            (10, 13),
            (12, 13)
        ]
    )

    g = GraphAlgo.DirectedGraph(tinydg.nv, tinydg.edges)

    @test g.numedges == tinydg.ne

    fp_result = GraphAlgo.search(g, 2, GraphAlgo.FindingPaths)
    @test length(fp_result.marked) == g.numvertices
    @test fp_result.numvisited == 1

    fp_result = GraphAlgo.search(g, 3, GraphAlgo.FindingPaths)
    @test fp_result.marked[1:6] == fill(true, 6)
    @test fp_result.marked[7:end] == fill(false, 7)

    fp_result = GraphAlgo.search(g, 3, GraphAlgo.FindingPaths, GraphAlgo.BreadthFirstSearch)
    @test fp_result.marked[1:6] == fill(true, 6)
    @test fp_result.marked[7:end] == fill(false, 7)

    # test whether breadth first search find the shortest path
    fp_result = GraphAlgo.search(g, 7, GraphAlgo.FindingPaths, GraphAlgo.DepthFirstSearch)
    @test (GraphAlgo.pathto(fp_result, 3) == [7, 1, 6, 5, 3] ||
           GraphAlgo.pathto(fp_result, 3) == [7, 5, 3])
    fp_result = GraphAlgo.search(g, 7, GraphAlgo.FindingPaths, GraphAlgo.BreadthFirstSearch)
    @test GraphAlgo.pathto(fp_result, 3) == [7, 5, 3]

    marked = GraphAlgo.search(g, [2, 3, 7], GraphAlgo.MultiFindingPaths)
    @test marked == Bool[1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1]

    g = GraphAlgo.DirectedGraph(minidg.nv, minidg.edges)
    cycle = GraphAlgo.search(g, GraphAlgo.Cycle)
    @test cycle.trace == [6, 5, 4, 6]
    @test cycle.hascycle == true

    g = GraphAlgo.DirectedGraph(minidg2.nv, minidg2.edges)
    cycle = GraphAlgo.search(g, GraphAlgo.Cycle)
    @test isnothing(cycle.trace)
    @test cycle.hascycle == false

    g = GraphAlgo.DirectedGraph(tinydg.nv, tinydg.edges)
    cycle = GraphAlgo.search(g, GraphAlgo.Cycle)
    @test isnothing(cycle.trace) == false
    @test cycle.hascycle == true

    # Test topological sort
    g = GraphAlgo.DirectedGraph(tinydg)
    result = GraphAlgo.search(g, GraphAlgo.Topological)
    @test isnothing(result.vorder)

    g = GraphAlgo.DirectedGraph(dag_p579)
    result = GraphAlgo.search(g, GraphAlgo.Topological)
    @test result.vorder == [9, 8, 3, 4, 1, 7, 10, 12, 13, 11, 6, 5, 2]

    g = GraphAlgo.DirectedGraph(tinydg)
    result = GraphAlgo.search(g, GraphAlgo.ConnectedComponents)
    @test Set(result.vids) == Set([1, 2, 3, 4, 5])
    @test length(Set(result.vids[10:13])) == 1
    @test length(Set(result.vids[[1, 3, 4, 5, 6]])) == 1
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
using Wavelets, DelimitedFiles, LinearAlgebra, PyPlot, MRIReco

function analyzeImage(x::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64};t0::Int64=size(D,2),tol=1e-3) where T
  nx,ny = xsize
  px,py = psize
  x = reshape(x,nx,ny)
  x_pad = repeat(x,2,2)[1:nx+px-1,1:ny+py-1] # pad image using periodic boundary conditions
  α = zeros(T,size(D,2),nx,ny)
  patch = zeros(T,px*py)
  for j=1:ny
    for i=1:nx
      patch[:] .= vec(x_pad[i:i+px-1,j:j+py-1])
      norm(patch)==0 && continue
      # matchingpursuit is contained in Wavelets.jl
      α[:,i,j] .= matchingpursuit(patch, x->D*x, x->transpose(D)*x, tol)
    end
  end
  return vec(α)
end

function synthesizeImage(α::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64}) where T
  nx,ny = xsize
  px,py = psize
  x = zeros(T,nx+px-1,ny+py-1)
  α = reshape(α,:,nx,ny)
  for j=1:ny
    for i=1:nx
      x[i:i+px-1,j:j+py-1] .+= reshape(D*α[:,i,j],px,py)
    end
  end
  return vec(x[1:nx,1:ny])/(px*py)
end

function dictOp(D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64},tol::Float64=1.e-3) where T
  produ = x->analyzeImage(x,D,xsize,psize,tol=tol)
  ctprodu = x->synthesizeImage(x,D,xsize,psize)
  return LinearOperator(prod(xsize)*size(D,2),prod(xsize),false,false
          , produ
          , nothing
          , ctprodu )
end

## To test our method, let us load some simulated data and subsample it

# phantom
img = readdlm("data/mribrain100.tsv")

acqData = AcquisitionData(ISMRMRDFile("data/acqDataBrainSim100.h5"))
nx,ny = acqData.encodingSize

# undersample kspace data
acqData = sample_kspace(acqData, 2.0, "poisson", calsize=25,profiles=false);

## Load the Dictionary, build the sparsifying transform and perform the reconstruction

# load the dictionary
D = ComplexF64.(readdlm("data/brainDict98.tsv"))

# some parameters
px, py = (6,6)  # patch size
K = px*py       # number of atoms

# CS reconstruction using Wavelets
params = Dict{Symbol,Any}()
params[:reco] = "standard"
params[:reconSize] = (nx,ny)
params[:iterations] = 50
params[:λ] = 2.e-2
params[:regularization] = "L1"
params[:sparseTrafo] = dictOp(D,(nx,ny),(px,py),2.e-2)
params[:ρ] = 0.1
params[:solver] = "admm"
params[:absTol] = 1.e-4
params[:relTol] = 1.e-2

img_d = reconstruction(acqData,params)
@info "relative error: $(norm(img-img_d)/norm(img))"

# For comparison, let us perform the same reconstruction as above but with a Wavelet transform

delete!(params, :sparseTrafo)
params[:sparseTrafoName] = "Wavelet"

img_w = reconstruction(acqData,params)
@info "relative error: $(norm(img-img_w)/norm(img))"

# use PyPlot for interactive display
figure(34)
clf()
subplot(2,2,1)
title("Original")
imshow(abs.(img[:,:,1,1,1]))
subplot(2,2,3)
title("Wavelet")
imshow(abs.(img_w[:,:,1,1,1]))
subplot(2,2,4)
title("Dictionary")
imshow(abs.(img_d[:,:,1,1,1]))

# export images
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/brainOrig.png")
exportImage(filename, abs.(img[:,:,1,1,1]) )
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/brainWavelet.png")
exportImage(filename, abs.(img_w[:,:,1,1,1]) )
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/brainDict.png")
exportImage(filename, abs.(img_d[:,:,1,1,1]) )


text:
using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z, default_indices
using Oceananigans.BoundaryConditions: CBC, PBC

struct XPartition{N} <: AbstractPartition
    div :: N
    function XPartition(sizes)
        if length(sizes) > 1 && all(y -> y == sizes[1], sizes)
            sizes = length(sizes)
        end
        return new{typeof(sizes)}(sizes)
    end
end

const EqualXPartition = XPartition{<:Number}

Base.length(p::XPartition)       = length(p.div)
Base.length(p::EqualXPartition)  = p.div

Base.summary(p::EqualXPartition) = "Equal partitioning in X with ($(p.div) regions)"
Base.summary(p::XPartition)      = "XPartition with [$(["$(p.div[i]) " for i in 1:length(p)]...)]"

function partition_size(p::EqualXPartition, grid)
    Nx, Ny, Nz = size(grid)
    @assert mod(Nx, p.div) == 0
    return Tuple((Nx ÷ p.div, Ny, Nz) for i in 1:length(p))
end

function partition_size(p::XPartition, grid)
    Nx, Ny, Nz = size(grid)
    @assert sum(p.div) != Nx
    return Tuple((p.div[i], Ny, Nz) for i in 1:length(p))
end

function partition_extent(p::XPartition, grid)
    x = cpu_face_constructor_x(grid)
    y = cpu_face_constructor_y(grid)
    z = cpu_face_constructor_z(grid)

    x = divide_direction(x, p)
    return Tuple((x = x[i], y = y, z = z) for i in 1:length(p))
end

function partition_topology(p::XPartition, grid)
    TX, TY, TZ = topology(grid)

    return Tuple(((TX == Periodic ? FullyConnected : i == 1 ?
                                    RightConnected : i == length(p) ?
                                    LeftConnected :
                                    FullyConnected), TY, TZ) for i in 1:length(p))
end

divide_direction(x::Tuple, p::EqualXPartition) =
    Tuple((x[1]+(i-1)*(x[2] - x[1])/length(p), x[1]+i*(x[2] - x[1])/length(p)) for i in 1:length(p))

function divide_direction(x::AbstractArray, p::EqualXPartition)
    nelem = (length(x)-1)÷length(p)
    return Tuple(x[1+(i-1)*nelem:1+i*nelem] for i in 1:length(p))
end

partition_global_array(a::Function, args...)  = a
partition_global_array(a::Field, p::EqualXPartition, args...) = partition_global_array(a.data, p, args...)

function partition_global_array(a::AbstractArray, ::EqualXPartition, local_size, region, arch)
    idxs = default_indices(length(size(a)))
    return arch_array(arch, a[local_size[1]*(region-1)+1:local_size[1]*region, idxs[2:end]...])
end

function partition_global_array(a::OffsetArray, ::EqualXPartition, local_size, region, arch)
    idxs    = default_indices(length(size(a)))
    offsets = (a.offsets[1], Tuple(0 for i in 1:length(idxs)-1)...)
    return arch_array(arch, OffsetArray(a[local_size[1]*(region-1)+1+offsets[1]:local_size[1]*region-offsets[1], idxs[2:end]...], offsets...))
end

#####
##### Global reconstruction utils
#####

function reconstruct_size(mrg, p::XPartition)
    Ny = mrg.region_grids[1].Ny
    Nz = mrg.region_grids[1].Nz
    Nx = sum([grid.Nx for grid in mrg.region_grids.regions])
    return (Nx, Ny, Nz)
end

function reconstruct_extent(mrg, p::XPartition)
    switch_device!(mrg.devices[1])
    y = cpu_face_constructor_y(mrg.region_grids.regions[1])
    z = cpu_face_constructor_z(mrg.region_grids.regions[1])

    if cpu_face_constructor_x(mrg.region_grids.regions[1]) isa Tuple
        x = (cpu_face_constructor_x(mrg.region_grids.regions[1])[1],
             cpu_face_constructor_x(mrg.region_grids.regions[length(p)])[end])
    else
        x = [cpu_face_constructor_x(mrg.region_grids.regions[1])...]
        for (idx, grid) in enumerate(mrg.region_grids.regions[2:end])
            switch_device!(mrg.devices[idx])
            x = [x..., cpu_face_constructor_x(grid)[2:end]...]
        end
    end
    return (; x = x, y = y, z = z)
end

const FunctionMRO     = MultiRegionObject{<:Tuple{Vararg{<:Function}}}
const ArrayMRO{T, N}  = MultiRegionObject{<:Tuple{Vararg{<:AbstractArray{T, N}}}} where {T, N}

reconstruct_global_array(ma::FunctionMRO, args...) = ma.regions[1]

function reconstruct_global_array(ma::ArrayMRO{T, N}, p::EqualXPartition, arch) where {T, N}
    local_size = size(first(ma.regions))
    global_Nx  = local_size[1] * length(p)
    idxs = default_indices(length(local_size))
    arr_out = zeros(eltype(first(ma.regions)), global_Nx, local_size[2:end]...)
    n = local_size[1]
    for r = 1:length(p)
        init = Int(n * (r - 1) + 1)
        fin  = Int(n * r)
        arr_out[init:fin, idxs[2:end]...] .= arch_array(CPU(), ma[r])[1:fin-init+1, idxs[2:end]...]
    end

    return arch_array(arch, arr_out)
end

function compact_data!(global_field, global_grid, data::MultiRegionObject, p::EqualXPartition)
    Nx, Ny, Nz = size(global_grid)
    n = Nx / length(p)
    for r = 1:length(p)
        init = Int(n * (r - 1) + 1)
        fin  = Int(n * r)
        interior(global_field)[init:fin, :, :] .= data[r][1:fin-init+1, :, :]
    end
    fill_halo_regions!(global_field)
end

#####
##### Boundary specific Utils
#####

struct Connectivity
    rank :: Int
    from_rank :: Int
end

inject_south_boundary(region, p::XPartition, bc) = bc
inject_north_boundary(region, p::XPartition, bc) = bc

function inject_west_boundary(region, p::XPartition, global_bc)
    if region == 1
        typeof(global_bc) <: Union{CBC, PBC} ?
                bc = CommunicationBoundaryCondition(Connectivity(region, length(p))) :
                bc = global_bc
    else
        bc = CommunicationBoundaryCondition(Connectivity(region, region - 1))
    end
    return bc
end

function inject_east_boundary(region, p::XPartition, global_bc)
    if region == length(p)
        typeof(global_bc) <: Union{CBC, PBC} ?
                bc = CommunicationBoundaryCondition(Connectivity(region, 1)) :
                bc = global_bc
    else
        bc = CommunicationBoundaryCondition(Connectivity(region, region + 1))
    end
    return bc
end

####
#### Global index flattening
####

@inline function displaced_xy_index(i, j, grid, region, p::XPartition)
    i′ = i + grid.Nx * (region - 1)
    t  = i′ + (j - 1) * grid.Nx * length(p)
    return t
end


text:
include("./Algos.jl")
using .Algos
using CSV
using DataFrames
using Plots

MAP_NAME = "map2"
TRAINING_EPS = 300000
EVALUATING_EPS = 100000
SAMPLING_RATE = 0.001

df = DataFrame(CSV.File("maps/$(MAP_NAME).csv"));
tiles = Matrix(df);
agent = Sarsa(tiles, learningRate = 0.1);
y1 = train(agent, TRAINING_EPS, samplingRate = SAMPLING_RATE);
Algos.visualizeOptimalPolicy(agent.env, agent.qTable, filename = "assets/sarsa_$(MAP_NAME)_optimal.png");

yt1 = Algos.evaluate(agent, EVALUATING_EPS, samplingRate = SAMPLING_RATE);

tiles = Matrix(df);
agent = QLearning(tiles, learningRate = 0.1);
y2 = train(agent, TRAINING_EPS, samplingRate = SAMPLING_RATE);
Algos.visualizeOptimalPolicy(agent.env, agent.qTable, filename = "assets/QL_$(MAP_NAME)_optimal.png");

yt2 = Algos.evaluate(agent, EVALUATING_EPS, samplingRate = SAMPLING_RATE);

fig = plot(collect(1:length(y1)) .* trunc(Int, 1 / SAMPLING_RATE), hcat(y1, y2), xlabel = "Episodes", ylabel = "Average Reward per Episode",
            label = ["SARSA" "Q-Learning"], lw = 2, dpi = 300, xlabelfontsize = 8, ylabelfontsize = 8,
            legend=:right, title = "Training", titlefontsize = 10);

savefig(fig, "assets/evaluation.png");

fig = plot(collect(1:length(yt1)) .* trunc(Int, 1 / SAMPLING_RATE), hcat(yt1, yt2), xlabel = "Episodes", ylabel = "Average Reward per Episode",
            label = ["SARSA" "Q-Learning"], lw = 2, dpi = 300, xlabelfontsize = 10, ylabelfontsize = 10,
            legend=:right, title = "Testing", titlefontsize = 10);

savefig(fig, "assets/evaluation_trained.png");


text:
using StatisticalRethinking, CmdStan, MCMCChain
gr(size=(500,500));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df)
first(df2, 5)

heightsmodel = "
// Inferring a Rate
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 20);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

stanmodel = Stanmodel(name="heights", model=heightsmodel);

heightsdata = Dict("N" => length(df2[:height]), "h" => df2[:height]);

rc, a3d, cnames = stan(stanmodel, heightsdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME);

pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
p = filter(p -> !(p in  pi), cnames)

chn = Chains(a3d,
  Symbol.(cnames),
  Dict(
    :parameters => Symbol.(p),
    :internals => Symbol.(pi)
  )
)

describe(chn)
describe(chn, section=:internals)


text:
### Linear Interpolation
struct LinearInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  LinearInterpolation{FT}(u,t) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t)
end
LinearInterpolation(u,t) = LinearInterpolation{true}(u,t)

### Quadratic Interpolation
struct QuadraticInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  QuadraticInterpolation{FT}(u,t) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t)
end
QuadraticInterpolation(u,t) = QuadraticInterpolation{true}(u,t)

### Lagrange Interpolation
struct LagrangeInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  n::Int
  LagrangeInterpolation{FT}(u,t,n) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t,n)
end

function LagrangeInterpolation(u,t,n=nothing)
  if n == nothing
    n = length(t) - 1 # degree
  end
  if n != length(t) - 1
    error("Currently only n=length(t) - 1 is supported")
  end
  LagrangeInterpolation{true}(u,t,n)
end

### ZeroSpline Interpolation
struct ZeroSpline{uType,tType,dirType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  dir::Symbol
  ZeroSpline{FT}(u,t,dir) where FT = new{typeof(u),typeof(t),typeof(dir),FT,eltype(u)}(u,t,dir)
end
 ZeroSpline(u,t;dir=:left) = ZeroSpline{true}(u,t,dir)

### QuadraticSpline Interpolation
struct QuadraticSpline{uType,tType,tAType,dType,zType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  tA::tAType
  d::dType
  z::zType
  QuadraticSpline{FT}(u,t,tA,d,z) where FT = new{typeof(u),typeof(t),typeof(tA),
                                                  typeof(d),typeof(z),FT,eltype(u)}(u,t,tA,d,z)
end

function QuadraticSpline(u,t)
  s = length(t)
  dl = ones(eltype(t),s-1)
  d = ones(eltype(t),s)
  du = zeros(eltype(t),s-1)
  tA = Tridiagonal(dl,d,du)
  d = zero(t)
  for i = 2:length(d)
    d[i] = 2//1 * (u[i] - u[i-1])/(t[i] - t[i-1])
  end
  z = tA\d
  QuadraticSpline{true}(u,t,tA,d,z)
end

# Cubic Spline Interpolation
struct CubicSpline{uType,tType,hType,zType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  h::hType
  z::zType
  CubicSpline{FT}(u,t,h,z) where FT = new{typeof(u),typeof(t),typeof(h),typeof(z),FT,eltype(u)}(u,t,h,z)
end

function CubicSpline(u,t)
  n = length(t) - 1
  h = vcat(0, diff(t), 0)
  dl = h[2:n+1]
  d = 2 .* (h[1:n+1] .+ h[2:n+2])
  du = h[2:n+1]
  tA = LinearAlgebra.Tridiagonal(dl,d,du)
  d = zero(t)
  for i = 2:n
    d[i] = 6(u[i+1] - u[i]) / h[i+1] - 6(u[i] - u[i-1]) / h[i]
  end
  z = tA\d
  CubicSpline{true}(u,t,h[1:n+1],z)
end

### BSpline Curve Interpolation
struct BSplineInterpolation{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  d::Int    # degree
  p::pType  # params vector
  k::kType  # knot vector
  c::cType  # control points
  pVecType::Symbol
  knotVecType::Symbol
  BSplineInterpolation{FT}(u,t,d,p,k,c,pVecType,knotVecType) where FT =  new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,p,k,c,pVecType,knotVecType)
end

function BSplineInterpolation(u,t,d,pVecType,knotVecType)
  n = length(t)
  s = zero(eltype(u))
  p = zero(t)
  k = zeros(eltype(t),n+d+1)
  l = zeros(eltype(u),n-1)
  p[1] = zero(eltype(t))
  p[end] = one(eltype(t))

  for i = 2:n
    s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2)
    l[i-1] = s
  end
  if pVecType == :Uniform
    for i = 2:(n-1)
      p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1)
    end
  elseif pVecType == :ArcLen
    for i = 2:(n-1)
      p[i] = p[1] + l[i-1]/s * (p[end]-p[1])
    end
  end

  lidx = 1
  ridx = length(k)
  while lidx <= (d+1) && ridx >= (length(k)-d)
    k[lidx] = p[1]
    k[ridx] = p[end]
    lidx += 1
    ridx -= 1
  end

  ps = zeros(eltype(t),n-2)
  s = zero(eltype(t))
  for i = 2:n-1
    s += p[i]
    ps[i-1] = s
  end

  if knotVecType == :Uniform
    # uniformly spaced knot vector
    # this method is not recommended because, if it is used with the chord length method for global interpolation,
    # the system of linear equations would be singular.
    for i = (d+2):n
      k[i] = k[1] + (i-d-1)//(n-d) * (k[end]-k[1])
    end
  elseif knotVecType == :Average
    # average spaced knot vector
    idx = 1
    if d+2 <= n
      k[d+2] = 1//d * ps[d]
    end
    for i = (d+3):n
      k[i] = 1//d * (ps[idx+d] - ps[idx])
      idx += 1
    end
  end
  # control points
  N = spline_coefficients(n,d,k,p)
  c = vec(N\u[:,:])
  BSplineInterpolation{true}(u,t,d,p,k,c,pVecType,knotVecType)
end

### BSpline Curve Approx
struct BSplineApprox{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  d::Int    # degree
  h::Int    # number of control points (n => h >= d >= 1)
  p::pType  # params vector
  k::kType  # knot vector
  c::cType  # control points
  pVecType::Symbol
  knotVecType::Symbol
  BSplineApprox{FT}(u,t,d,h,p,k,c,pVecType,knotVecType) where FT =  new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,h,p,k,c,pVecType,knotVecType)
end

function BSplineApprox(u,t,d,h,pVecType,knotVecType)
  n = length(t)
  s = zero(eltype(u))
  p = zero(t)
  k = zeros(eltype(t),h+d+1)
  l = zeros(eltype(u),n-1)
  p[1] = zero(eltype(t))
  p[end] = one(eltype(t))

  for i = 2:n
    s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2)
    l[i-1] = s
  end
  if pVecType == :Uniform
    for i = 2:(n-1)
      p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1)
    end
  elseif pVecType == :ArcLen
    for i = 2:(n-1)
      p[i] = p[1] + l[i-1]/s * (p[end]-p[1])
    end
  end

  lidx = 1
  ridx = length(k)
  while lidx <= (d+1) && ridx >= (length(k)-d)
    k[lidx] = p[1]
    k[ridx] = p[end]
    lidx += 1
    ridx -= 1
  end

  ps = zeros(eltype(t),n-2)
  s = zero(eltype(t))
  for i = 2:n-1
    s += p[i]
    ps[i-1] = s
  end

  if knotVecType == :Uniform
    # uniformly spaced knot vector
    # this method is not recommended because, if it is used with the chord length method for global interpolation,
    # the system of linear equations would be singular.
    for i = (d+2):h
      k[i] = k[1] + (i-d-1)//(h-d) * (k[end]-k[1])
    end
  elseif knotVecType == :Average
    # NOTE: verify that average method can be applied when size of k is less than size of p
    # average spaced knot vector
    idx = 1
    if d+2 <= h
      k[d+2] = 1//d * ps[d]
    end
    for i = (d+3):h
      k[i] = 1//d * (ps[idx+d] - ps[idx])
      idx += 1
    end
  end
  # control points
  c = zeros(eltype(u), h)
  c[1] = u[1]
  c[end] = u[end]
  q = zeros(eltype(u), n)
  N = zeros(eltype(t), n, h)
  for i = 1:n
    N[i, :] .= spline_coefficients(h,d,k,p[i])
  end
  for k = 2:n-1
    q[k] = u[k] - N[k, 1]*u[1] - N[k, h]*u[end]
  end
  Q = Matrix{eltype(u)}(undef, h-2, 1)
  for i = 2:h-1
    s = 0.0
    for k = 2:n-1
      s += N[k, i]*q[k]
    end
    Q[i-1] = s
  end
  N = N[2:end-1,2:h-1]
  M = transpose(N) * N
  P = M\Q
  c[2:end-1] .= vec(P)
  BSplineApprox{true}(u,t,d,h,p,k,c,pVecType,knotVecType)
end

### Loess
struct Loess{uType,tType,αType,xType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  d::Int
  α::αType
  q::Int
  x::xType
  Loess{FT}(u,t,d,α,q,x) where FT = new{typeof(u),typeof(t),typeof(α),typeof(x),FT,eltype(u)}(u,t,d,α,q,x)
end

function Loess(u,t,d,α)
  n = length(t)
  q = floor(Int,n*α)
  x = Matrix{eltype(t)}(undef,n,d+1)
  x[:,1] .= one(t[1])
  for i = 2:(d+1)
    x[:,i] = t .^ (i-1)
  end
  Loess{true}(u,t,d,α,q,x)
end

### GaussianProcesses
struct GPInterpolation{uType,tType,gpType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  gp::gpType
  GPInterpolation{FT}(u,t,gp) where FT = new{typeof(u),typeof(t),typeof(gp),FT,eltype(u)}(u,t,gp)
end

function GPInterpolation(u,t,m,k,n=-2.0)
  gp = GP(t,u,m,k,n)
  GPInterpolation{true}(u,t,gp)
end

### Curvefit
struct Curvefit{uType,tType,mType,cfType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  m::mType
  c_f::cfType
  Curvefit{FT}(u,t,m,c_f) where FT = new{typeof(u),typeof(t),typeof(m),typeof(c_f),FT,eltype(u)}(u,t,m,c_f)
end

function Curvefit(u,t,m,p)
  c_f = curve_fit(m,t,u,p)
  Curvefit{true}(u,t,m,c_f)
end


text:
function getProblemSpecificParameters(nx, ny, nquadpoints, ThisICType, pointsxyz)

    if ThisICType == "linesource"
        SigmaS = 1 * ones(ny + 4, nx + 4)
        SigmaA = 0 * ones(ny + 4, nx + 4)
        SigmaT = SigmaS + SigmaA
        x0 = -1.5
        x1 = 1.5
        y0 = -1.5
        y1 = 1.5
        dx = (x1 - x0) / (nx)
        dy = (y1 - y0) / (ny)
        tEnd = 1.0
        phi0 = setInitialConditionsGaussian(nx, ny, nquadpoints, dx, dy, x0, y0)
        Source = zeros(Float64, nquadpoints, ny + 4, nx + 4)
    elseif ThisICType == "beam"
        SigmaS = 0.0 * ones(ny + 4, nx + 4)
        SigmaA = 0.0 * ones(ny + 4, nx + 4)
        SigmaT = SigmaS + SigmaA
        x0 = -1.5
        x1 = 1.5
        y0 = -1.5
        y1 = 1.5
        dx = (x1 - x0) / (nx)
        dy = (y1 - y0) / (ny)
        tEnd = 10.0

        beamdirection = [1 / sqrt(2), 1 / sqrt(2), 0]
        beamdirection = [0, 1, 0]
        Source = setInitialConditionsBeam(
            nx,
            ny,
            nquadpoints,
            dx,
            dy,
            x0,
            y0,
            pointsxyz,
            beamdirection,
        )
        phi0 = zeros(Float64, nquadpoints, ny + 4, nx + 4)
        v = 10.0
        for j = 3:nx+2
            for i = 3:ny+2
                y = y0 + dy / 2 + (i - 3) * dy
                x = x0 + dx / 2 + (j - 3) * dx
                if x > 0.0 #&& x <0.1
                    SigmaS[i, j] = v
                end

                #				if x>0.4 && x<0.5
                #					SigmaT[i,j] = v
                #				end
                #
                #
                #				if x>0.7 && x<0.8
                #					SigmaT[i,j] = v
                #				end
            end
        end
        SigmaT = SigmaS + SigmaA

    elseif ThisICType == "checkerboard"
        x0 = 0.0
        x1 = 7.0
        y0 = 0.0
        y1 = 7.0
        dx = (x1 - x0) / (nx)
        dy = (y1 - y0) / (ny)
        SigmaS = zeros(ny + 4, nx + 4)
        SigmaT = zeros(ny + 4, nx + 4)
        phi0 = zeros(Float64, nquadpoints, ny + 4, nx + 4)
        tEnd = 3.2
        Source = zeros(Float64, nquadpoints, ny + 4, nx + 4)

        for j = 3:nx+2
            for i = 3:ny+2
                y = y0 + dy / 2 + (i - 3) * dy
                x = x0 + dx / 2 + (j - 3) * dx
                cx = ceil(x)
                cy = ceil(y)
                g =
                    float(ceil((cx + cy) / 2.0) * 2.0 == (cx + cy)) * (
                        (1.0 < cx) &&
                        (cx < 7.0) &&
                        (1.0 < cy) &&
                        ((cy - 2.0 * abs(cx - 4.0)) < 4.0)
                    )
                SigmaS[i, j] = (1 - g) * 1 + g * 0
                SigmaT[i, j] = (1 - g) * 0 + g * 10 + SigmaS[i, j]
                if ((3 < x) && (x < 4) && (3 < y) && (y < 4))
                    for q = 1:nquadpoints
                        Source[q, i, j] = 1.0
                    end
                    SigmaS[i, j] = 0.0
                    SigmaT[i, j] = 10.0
                end
            end
        end
    elseif ThisICType == "dim1"
        SigmaS = 1 * ones(ny + 4, nx + 4)
        SigmaA = 0 * ones(ny + 4, nx + 4)
        SigmaT = SigmaS + SigmaA
        x0 = -1.5
        x1 = 1.5
        y0 = -1.5
        y1 = 1.5
        dx = (x1 - x0) / (nx)
        dy = (y1 - y0) / (ny)
        tEnd = 10.0
        phi0 = setInitialConditionsGaussian1d(nx, ny, nquadpoints, dx, dy, x0, y0)
        Source = zeros(Float64, nquadpoints, ny + 4, nx + 4)
    else
        #warning('Not yet implemented')
    end
    return x0, x1, y0, y1, tEnd, SigmaS, SigmaT, Source, phi0
end

function setInitialConditionsGaussian(nx, ny, nq, dx, dy, domainx0, domainy0)
    # init to zero
    phi = zeros(Float64, nq, ny + 4, nx + 4)
    x0 = 0.0
    y0 = 0.0
    s2 = 0.03^2
    floor = 1e-4
    f(x, y) = max(
        floor,
        1.0 / (4.0 * pi * s2) *
        exp(-((x - x0) * (x - x0) + (y - y0) * (y - y0)) / 4.0 / s2),
    )
    for j = 3:nx+2
        for i = 3:ny+2
            y = domainy0 + dy / 2 + (i - 3) * dy
            x = domainx0 + dx / 2 + (j - 3) * dx
            for q = 1:nq
                phi[q, i, j] = f(x, y) / 4.0 / pi
            end
        end
    end
    return phi
end

function setInitialConditionsGaussian1d(nx, ny, nq, dx, dy, domainx0, domainy0)
    # init to zero
    phi = zeros(Float64, nq, ny + 4, nx + 4)
    x0 = 0.0
    y0 = 0.0
    s2 = 0.03^2
    floor = 1e-4
    f(x, y) = max(
        floor,
        1.0 / (4.0 * pi * s2) *
        exp(-((x - x0) * (x - x0) + (y - y0) * (y - y0)) / 4.0 / s2),
    )
    for j = 3:nx+2
        for i = 3:ny+2
            y = domainy0 + dy / 2 + (i - 3) * dy
            x = domainx0 + dx / 2 + (j - 3) * dx
            for q = 1:nq
                phi[q, i, j] = f(x, 0) / 4.0 / pi
            end
        end
    end
    return phi
end

function setInitialConditionsBeam(
    nx,
    ny,
    nq,
    dx,
    dy,
    domainx0,
    domainy0,
    quadpoints,
    beamdirection,
)
    # init to zero
    phi = 1e-14 .* ones(Float64, nq, ny + 4, nx + 4)
    # beam postion:
    x0 = -0.50
    y0 = 0.00
    # beam width, height
    height = 0.05
    width = 0.05

    for j = 3:nx+2
        for i = 3:ny+2
            y = domainy0 + dy / 2 + (i - 3) * dy
            x = domainx0 + dx / 2 + (j - 3) * dx
            if x > x0 - width && x < x0 + width && y > y0 - height && y < y0 + height # beam should be here!
                # beam is a gaussian with pole pointing into beamdirection
                for q = 1:nq
                    Omega = quadpoints[q, :]
                    eps = 1^2
                    x = dot(Omega, beamdirection) - 1
                    val = 1.0 / sqrt(2.0 * pi * eps) * exp(-x^2 / (2 * eps))
                    phi[q, i, j] = val
                end
            end
        end
    end
    return phi
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
# This file is a part of SimilaritySearch.jl

export LocalSearchAlgorithm, SearchGraph, SearchGraphPools, SearchGraphCallbacks, VisitedVertices, index!, push_item!
export Neighborhood, IdentityNeighborhood, DistalSatNeighborhood, SatNeighborhood
export find_neighborhood, push_neighborhood!
export callbacks
export BeamSearch, BeamSearchSpace, Callback
export KDisjointHints, DisjointHints, RandomHints

include("graph.jl")
include("callbacks.jl")
include("rebuild.jl")
include("insertions.jl")

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch