text:
#=

Burning phase for tribe.

Ignacio Quintero Mächler

t(-_-t)

April 27 2017

=#

"""
    burn_compete(...)

Burning & adaptive phase for MCMC.
"""
function burn_tribe(total_llf     ::Function,
                    ω1upd_llr     ::Function,
                    ω0upd_llr     ::Function,
                    λ1upd_llr     ::Function,
                    λ0upd_llr     ::Function,
                    Xupd_llr      ::Function,
                    Rupd_llr      ::Function,
                    ωxupd_llr     ::Function,
                    σ²upd_llr     ::Function,
                    bgiid         ::Function,
                    bgiid_br      ::Function,
                    mhr_upd_Xbr   ::Function,
                    mhr_upd_Xtrio ::Function,
                    mhr_upd_Ybr   ::Function,
                    mhr_upd_Ytrio ::Function,
                    mhr_upd_Ystem ::Function,
                    mhr_upd_XYbr  ::Function,
                    mhr_upd_XYtrio::Function,
                    Xc      ::Array{Float64,2},
                    Yc      ::Array{Int64,3},
                    δXc     ::Array{Float64,3},
                    δYc     ::Array{Float64,3},
                    LApc    ::Array{Float64,2},
                    LAnc    ::Array{Float64,2},
                    LDc     ::Array{Float64,3},
                    σ²c     ::Float64,
                    ωxc     ::Float64,
                    ω1c     ::Float64,
                    ω0c     ::Float64,
                    λ1c     ::Float64,
                    λ0c     ::Float64,
                    Xnc1    ::Array{Int64,1},
                    Xnc2    ::Array{Int64,1},
                    brl     ::Array{Float64,1},
                    wcol    ::Array{Array{Int64,1},1},
                    bridx_a ::Array{Array{UnitRange{Int64},1},1},
                    brs     ::Array{Int64,3},
                    stemevc ::Array{Array{Float64,1},1},
                    trios   ::Array{Array{Int64,1},1},
                    wXp     ::Array{Int64,1},
                    λprior  ::Float64,
                    ωxprior ::NTuple{2,Float64},
                    ω1prior ::NTuple{2,Float64},
                    ω0prior ::NTuple{2,Float64},
                    σ²prior ::Float64,
                    np      ::Int64,
                    parvec  ::Array{Int64,1},
                    nburn   ::Int64,
                    obj_ar  ::Float64 = 0.234,
                    tune_int::Int64   = 100)

  m, ntip, narea = size(Yc)

  nedge = size(brs, 1)

  # likelihood and prior
  if ωxc >= 0.0
    llc = total_llf(Xc, Yc, LApc, LDc, ωxc, ω1c, ω0c, λ1c, λ0c,
                    stemevc, brs, σ²c)
  else
    llc = total_llf(Xc, Yc, LAnc, LDc, ωxc, ω1c, ω0c, λ1c, λ0c,
                    stemevc, brs, σ²c)
  end

  prc = allλpr(   λ1c, λ0c, λprior)           +
        logdexp(  σ²c, σ²prior)               +
        logdnorm(ωxc, ωxprior[1], ωxprior[2]) +
        logdnorm(ω1c, ω1prior[1], ω1prior[2]) +
        logdnorm(ω0c, ω0prior[1], ω0prior[2])

  # make scaling function
  scalef = makescalef(obj_ar)

  # rest of tuning parameters
  ptn = fill(0.1, np)

  # initialize acceptance log
  ltn = zeros(Int64, np)
  lup = zeros(Int64, np)
  lac = zeros(Int64, np)

  # row i proposals for X
  xpi  = fill(NaN, ntip)             # proposal x slice
  δxi  = fill(NaN, ntip, ntip)       # lineage pairwise differences
  lapi = fill(NaN, ntip)             # lineage average for ωx > 0
  lani = fill(NaN, ntip)             # lineage average for ωx < 0
  ldi  = fill(NaN, ntip, narea)      # lineage difference

  rj = wcol[1][2]

  # progress bar
  p = Progress(nburn, dt=5, desc="burn...", barlen=20, color=:green)

  # print number of parameters
  printstyled(
    "\nωx updates per iter = ", lastindex(filter(x -> x == 2,parvec)),
    "\nω1 updates per iter = ", lastindex(filter(x -> x == 3,parvec)),
    "\nω0 updates per iter = ", lastindex(filter(x -> x == 4,parvec)),
    "\nσ² updates per iter = ", lastindex(filter(x -> x == 1,parvec)),
    "\nλ1 updates per iter = ", lastindex(filter(x -> x == 5,parvec)),
    "\nλ0 updates per iter = ", lastindex(filter(x -> x == 6,parvec)),
    "\nθ  updates per iter = ", length(parvec), "\n", color = :green)

  Xcidx = CartesianIndices(Xc)

  #start burnin
  for it = Base.OneTo(nburn)

    # Update vector
    shuffle!(parvec)

    for up = parvec

      # update X
      if up > 6

        # X updates
        @inbounds begin

          upx = wXp[up - 6]

          # if root
          if upx == 1

            # allocate
            @simd for i = Base.OneTo(ntip)
              xpi[i] =  Xc[1,i]
            end

            # update xi
            addupt!(xpi, ptn, 1, up)

            xpi[rj] = xpi[1]::Float64

            llr = Rupd_llr(xpi, Xc, σ²c)::Float64

            if -randexp() < llr
              llc     += llr::Float64
              Xc[1,:]  = xpi::Array{Float64,1}
              lac[up] += 1   # log acceptance
            end

          else

            xi, xj = Xcidx[upx].I

            # allocate
            for j = Base.OneTo(ntip)
              xpi[j]  = Xc[xi,j]
              lapi[j] = LApc[xi,j]
              lani[j] = LAnc[xi,j]
              @simd for k = Base.OneTo(narea)
                ldi[j,k] = LDc[xi,j,k]
              end
              @simd for i = Base.OneTo(ntip)
                δxi[i,j] = δXc[i,j,xi]
              end
            end

            # update xi
            addupt!(xpi, ptn, xj, up)

            if in(upx, Xnc1)        # if an internal node
              xpi[Xcidx[Xnc2[findfirst(isequal(upx),Xnc1)]][2]] = xpi[xj]
            end

            # calculate new averages
            Xupd_linavg!(δxi, lapi, lani, ldi, wcol,
                         xpi, xi, xj, Yc, δYc, narea)

            if ωxc >= 0.0
              llr = Xupd_llr(xi, xpi, Xc, lapi, ldi, LApc, LDc, Yc,
                             ωxc, ω1c, ω0c, λ1c, λ0c, σ²c)::Float64
            else
              llr = Xupd_llr(xi, xpi, Xc, lani, ldi, LAnc, LDc, Yc,
                             ωxc, ω1c, ω0c, λ1c, λ0c, σ²c)::Float64
            end

            if -randexp() < llr
              llc        += llr::Float64
              Xc[xi,:]    = xpi::Array{Float64,1}
              δXc[:,:,xi] = δxi::Array{Float64,2}
              LApc[xi,:]  = lapi::Array{Float64,1}
              LAnc[xi,:]  = lani::Array{Float64,1}
              LDc[xi,:,:] = ldi::Array{Float64,2}
              lac[up]    += 1   # log acceptance
            end
          end
        end

      # update λ1
      elseif up == 5

        # update λ
        λ1p = mulupt(λ1c, ptn[5])::Float64

        llr = λ1upd_llr(Yc, λ1c, λ1p, λ0c, ω1c, LDc, stemevc, brs)

        prr = llrdexp_x(λ1p, λ1c, λprior)

        if -randexp() < (llr + prr + log(λ1p/λ1c))
          llc    += llr::Float64
          prc    += prr::Float64
          λ1c     = λ1p::Float64
          lac[5] += 1
        end

        # update internal node
        if rand() < 0.4
          λϕ1, λϕ0 = λϕprop()
          llc = mhr_upd_Ytrio(rand(trios), Xc, Yc,
                      λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, λϕ1, λϕ0, llc, prc,
                      LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
        end

      # update λ0
      elseif up == 6

        # update λ
        λ0p = mulupt(λ0c, ptn[6])::Float64

        llr = λ0upd_llr(Yc, λ1c, λ0c, λ0p, ω0c, LDc, stemevc, brs)::Float64

        prr = llrdexp_x(λ0p, λ0c, λprior)

        if -randexp() < (llr + prr + log(λ0p/λ0c))
          llc     += llr::Float64
          prc     += prr::Float64
          λ0c      = λ0p::Float64
          lac[6] += 1
        end

        # update internal node
        if rand() < 0.4
          λϕ1, λϕ0 = λϕprop()
          llc = mhr_upd_Ytrio(rand(trios), Xc, Yc,
                    λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, λϕ1, λϕ0, llc, prc,
                    LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
        end

      # if σ² is updated
      elseif up == 1

        σ²p = mulupt(σ²c, ptn[1])::Float64

        #likelihood ratio
        if ωxc >= 0.0
          llr = σ²upd_llr(Xc, LApc, ωxc, σ²c, σ²p)::Float64
        else
          llr = σ²upd_llr(Xc, LAnc, ωxc, σ²c, σ²p)::Float64
        end

        # prior ratio
        prr = llrdexp_x(σ²p, σ²c, σ²prior)

        if -randexp() < (llr + prr + log(σ²p/σ²c))
          llc += llr::Float64
          prc += prr::Float64
          σ²c  = σ²p::Float64
          lac[1] += 1
        end

      #update ωx
      elseif up == 2

        ωxp = addupt(ωxc, ptn[2])::Float64

        #likelihood ratio
        llr = ωxupd_llr(Xc, LApc, LAnc, ωxc, ωxp, σ²c)

        # prior ratio
        prr = llrdnorm_x(ωxp, ωxc, ωxprior[1], ωxprior[2])

        if -randexp() < (llr + prr)
          llc += llr::Float64
          prc += prr::Float64
          ωxc  = ωxp::Float64
          lac[2] += 1
        end

      #update ω1
      elseif up == 3

        ω1p = addupt(ω1c, ptn[3])::Float64

        # proposal likelihood and prior
        llr = ω1upd_llr(Yc, λ1c, ω1c, ω1p, LDc)

        # prior ratio
        prr = llrdnorm_x(ω1p, ω1c, ω1prior[1], ω1prior[2])

        if -randexp() < (llr + prr)
          llc += llr::Float64
          prc += prr::Float64
          ω1c  = ω1p::Float64
          lac[3] += 1
        end

        # update internal node
        if rand() < 0.4
          λϕ1, λϕ0 = λϕprop()
          llc = mhr_upd_Ytrio(rand(trios), Xc, Yc,
                    λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, λϕ1, λϕ0, llc, prc,
                    LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
        end

      # update ω0
      else

        ω0p = addupt(ω0c, ptn[4])

        # proposal likelihood and prior
        llr = ω0upd_llr(Yc, λ0c, ω0c, ω0p, LDc)::Float64

        # prior ratio
        prr = llrdnorm_x(ω0p, ω0c, ω0prior[1], ω0prior[2])

        if -randexp() < (llr + prr)
          llc += llr
          prc += prr
          ω0c  = ω0p
          lac[4] += 1
        end

        # update internal node
        if rand() < 0.4
          λϕ1, λϕ0 = λϕprop()
          llc = mhr_upd_Ytrio(rand(trios), Xc, Yc,
                    λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, λϕ1, λϕ0, llc, prc,
                    LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
        end

      end

      ## make DA updates with some Pr
      # make Y branch update
      if rand() < 2e-3
        λϕ1, λϕ0 = λϕprop()
        llc = mhr_upd_Ybr(rand(Base.OneTo(nedge)),
                        Xc, Yc, λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, λϕ1, λϕ0,
                        llc, prc, LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
      end

      # make X branch update
      if rand() < 2e-3
        σ²ϕ = σ²ϕprop()
        llc = mhr_upd_Xbr(rand(Base.OneTo(nedge-1)),
                          Xc, Yc, λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, σ²ϕ, llc,
                          LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
      end

      # make X trio update
      if rand() < 2e-3
        σ²ϕ = σ²ϕprop()
        llc = mhr_upd_Xtrio(rand(trios),
                            Xc, Yc, λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, σ²ϕ, llc,
                            LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
      end

      # make joint X & Y branch update
      if rand() < 2e-3
        λϕ1, λϕ0 = λϕprop()
        σ²ϕ = σ²ϕprop()
        llc = mhr_upd_XYbr(rand(Base.OneTo(nedge-1)),
                           Xc, Yc, λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, σ²ϕ, λϕ1, λϕ0,
                           llc, prc, LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
      end

      # make joint X & Y trio update
      if rand() < 2e-3
        λϕ1, λϕ0 = λϕprop()
        σ²ϕ = σ²ϕprop()
        llc      = mhr_upd_XYtrio(rand(trios),
                           Xc, Yc, λ1c, λ0c, ωxc, ω1c, ω0c, σ²c, σ²ϕ, λϕ1, λϕ0,
                           llc, prc, LApc, LAnc, LDc, δXc, δYc, brs, stemevc)
      end

      # update stem branch
      if rand() < 2e-3
        λϕ1, λϕ0 = λϕprop()
        llc = mhr_upd_Ystem(λ1c, λ0c, λϕ1, λϕ0, llc, stemevc, brs)
      end

      # log number of updates
      ltn[up] += 1
      lup[up] += 1

      if (in(tune_int,ltn))
        wts = findall(isequal(tune_int),ltn)      # which to scale
        for j = wts
          ar     = lac[j]/lup[j]
          ptn[j] = scalef(ptn[j],ar)
          ltn[j] = 0
        end
      end
    end

    next!(p)
  end

  return llc, prc, Xc, Yc, LApc, LAnc, LDc, δXc, δYc, stemevc,
         brs, σ²c, ωxc, ω1c, ω0c, λ1c, λ0c, ptn
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
"""
Struct that holds the boundary arrays in grid-point space

    ϕ0::Array{NF,2}              # surface geopotential [m^2/s^2]
    ϕ0trunc::Array{NF,2}         # spectrally truncated surface geopotential [m^2/s^2]
    land_sea_mask::Array{NF,2}   # land-sea mask
    albedo::Array{NF,2}          # annual mean surface albedo
"""
struct Boundaries{NF<:AbstractFloat}    # number format NF
    ϕ0::Array{NF,2}              # surface geopotential (i.e. orography) [m^2/s^2]
    ϕ0trunc::Array{NF,2}         # spectrally truncated surface geopotential [m^2/s^2]
    land_sea_mask::Array{NF,2}   # land-sea mask
    albedo::Array{NF,2}          # annual mean surface albedo
end

""" Generator function for a Boundaries struct. Loads the boundary conditions,
orography, land-sea mask and albedo from an netCDF file and stores the in a
Boundaries-struct."""
function Boundaries{NF}(P::Params,
                        G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    @unpack boundary_path, boundary_file, g = P

    # LOAD NETCDF FILE
    if boundary_path == ""
        path = joinpath(@__DIR__,"../input_data",boundary_file)
    else
        path = joinpath(boundary_path,boundary_file)
    end
    ncfile = NetCDF.open(path)
    orog = ncfile.vars["orog"][:,end:-1:1]  # latitude is North to South in file
    lsm = ncfile.vars["lsm"][:,end:-1:1]
    alb = ncfile.vars["alb"][:,end:-1:1]

    # GEOPOTENTIAL, truncate in spectral space
    ϕ0 = g*orog
    ϕ0trunc = spectral_truncation(ϕ0,G)

    Boundaries{NF}(ϕ0,ϕ0trunc,lsm,alb)
end

ext: jl


text:
# # Trade-off curves
using Random
Random.seed!(1)
m = 25;
n = 10;
A = randn(m, n);
b = randn(m, 1);

#-

using Convex, SCS, LinearAlgebra

gammas = exp10.(range(-4, stop=2, length=100));

x_values = zeros(n, length(gammas));
x = Variable(n);
for i=1:length(gammas)
    cost = sumsquares(A*x - b) + gammas[i]*norm(x,1);
    problem = minimize(cost, [norm(x, Inf) <= 1]);
    solve!(problem, SCS.Optimizer(verbose=0));
    x_values[:,i] = evaluate(x);
end

#-

# Plot the regularization path.

using Plots
plot(title = "Entries of x vs lambda", xaxis=:log, xlabel="lambda", ylabel="x" )
for i = 1:n
    plot!(gammas, x_values[i,:], label="x$i")
end
plot!()


text:
function automatic_parallel_type()
    MPI.REFCOUNT[] == -1 && MPI.Init()

    if MPI.Comm_size(MPI.COMM_WORLD) == 1
        return NotParallel()
    else
        return ParallelMPI()
    end
end

function _iterator_blocks(interval, rank, n_par)
    rank         = rank -1
    n_min, extra = divrem(length(interval), n_par)
    iter_length  = n_min + (rank < extra ? 1 : 0)
    iter_start   = n_min * rank + min(rank, extra) + 1
    iter_end     = iter_start + iter_length - 1

    return iter_start:iter_end
end

function num_workers end

function worker_local_seed end

function worker_block end


text:
# ---
# title: 747. Largest Number At Least Twice of Others
# id: problem747
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/largest-number-at-least-twice-of-others/description/>
# hidden: true
# ---
#
# In a given integer array `nums`, there is always exactly one largest element.
#
# Find whether the largest element in the array is at least twice as much as
# every other number in the array.
#
# If it is, return the **index** of the largest element, otherwise return -1.
#
# **Example 1:**
#
#
#
#     Input: nums = [3, 6, 1, 0]
#     Output: 1
#     Explanation: 6 is the largest integer, and for every other number in the array x,
#     6 is more than twice as big as x.  The index of value 6 is 1, so we return 1.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1, 2, 3, 4]
#     Output: -1
#     Explanation: 4 isn't at least as big as twice the value of 3, so we return -1.
#
#
#
#
# **Note:**
#
#   1. `nums` will have a length in the range `[1, 50]`.
#   2. Every `nums[i]` will be an integer in the range `[0, 99]`.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
# Speed test for decomposing matrices in Julia.

function speed_test(func, orderN = 1000)
  sizes = Any[]
  for i in 0:int(log10(orderN))
    push!(sizes, 1 * 10^i)
    push!(sizes, 2 * 10^i)
    push!(sizes, 5 * 10^i)
  end

  y = zeros(length(sizes))
  for i in 1:length(sizes)
    n = sizes[i]
    A = randn(n, n)
    A = A'*A
    print("Testing with n = $(n)")
    tic()
    res = func(A)
    y[i] = toq()
    println(", time = $(y[i]), ")
  end

  return sizes, y
end

n_eig, t_eig = speed_test(eig, 1000)
n_chol, t_chol = speed_test(chol, 1000)

# Linear model
lin_eig = linreg(t_eig, n_eig)


text:
module SpatialData

using NamedDims, ImageCore, AxisIndices, NamedIndicesMeta

using NamedIndicesMeta.ColorData

using NamedIndicesMeta.TimeData

using NamedIndicesMeta.ObservationData

using Base: tail

export
    spatial_order,
    spatialdims,
    spatial_axes,
    spatial_units,
    spatial_offset,
    spatial_keys,
    spatial_indices,
    spatial_size,
    spatial_directions,
    spatial_units,
    pixel_spacing

"""
    spatial_order(x) -> Tuple{Vararg{Symbol}}
"""
spatial_order(x::X) where {X} = _spatial_order(Val(dimnames(X)))
@generated function _spatial_order(::Val{L}) where {L}
    # 0-Allocation see:
    # `@btime (()->remaining_dimnames_from_indexing((:a, :b, :c), (:,390,:)))()``
    keep_names = []
    for n in L
        if !(is_time(n) | is_color(n) | is_observation(n))
            push!(keep_names, n)
        end
    end
    out = (keep_names...,)
    quote
        return $out
    end
end

"""
    spatialdims(x) -> Tuple{Vararg{Int}}

Return a tuple listing the spatial dimensions of `img`.
Note that a better strategy may be to use ImagesAxes and take slices along the time axis.
"""
@inline spatialdims(x) = dim(dimnames(x), spatial_order(x))

"""
    spatial_axes(x) -> Tuple
"""
@inline spatial_axes(x) = _spatial_axes(named_axes(x), spatial_order(x))
function _spatial_axes(na::NamedTuple, spo::Tuple{Vararg{Symbol}})
    return map(spo_i -> getfield(na, spo_i), spo)
end

"""
    spatial_size(x) -> Tuple{Vararg{Int}}

Return a tuple listing the sizes of the spatial dimensions of the image.
"""
@inline spatial_size(x) = map(length, spatial_axes(x))

"""
    spatial_indices(x)

Return a tuple with the indices of the spatial dimensions of the
image. Defaults to the same as `indices`, but using ImagesAxes you can
mark some axes as being non-spatial.
"""
@inline spatial_indices(x) = map(values, spatial_axes(x))

"""
    spatial_keys(x)
"""
@inline spatial_keys(x) = map(keys, spatial_axes(x))

"""
    pixel_spacing(x)
"""
@inline pixel_spacing(x) = map(step, spatial_keys(x))
#@inline ImageCore.pixel_spacing(x::NamedDimsArray) = map(step, spatial_keys(x))

ImageCore.spacedirections(x::NamedDimsArray) = ImageCore._spacedirections(pixel_spacing(x))

spatial_eltype(x) = map(eltype, spatial_axes(x))

"""
    spatial_units(x)

Returns the units (i.e. Unitful.unit) that each spatial axis is measured in. If not
available `nothing` is returned for each spatial axis.
"""
spatial_units(x) = map(unit, spatial_eltype(x))

"""
    spatial_offset(x)

The offset of each dimension (i.e., where each spatial axis starts).
"""
spatial_offset(x) = map(first, spatial_keys(x))

"""
    spatial_directions(img) -> (axis1, axis2, ...)

Return a tuple-of-tuples, each `axis[i]` representing the displacement
vector between adjacent pixels along spatial axis `i` of the image
array, relative to some external coordinate system ("physical
coordinates").

By default this is computed from `pixel_spacing`, but you can set this
manually using ImagesMeta.
"""
@inline function spatial_directions(x::X) where {X}
    return _spatial_directions(HasProperties(X), x)
end

function _spatial_directions(::HasProperties{false}, x::AbstractArray{T,N}) where {T,N}
    return _default_spatial_directions(pixel_spacing(x))
end
function _spatial_directions(::HasProperties{true}, x::AbstractArray{T,N}) where {T,N}
    if hasproperty(x, :spatial_directions)
        return getproperty(x, :spatial_directions)
    else
        return _default_spatial_directions(pixel_spacing(x))
    end
end
function _default_spatial_directions(ps::NTuple{N,Any}) where N
    return ntuple(i->ntuple(d->d==i ? ps[d] : zero(ps[d]), Val(N)), Val(N))
end

end


text:
using Test
using JSON
using MianalyzerBackend
const GaussianMixturesController = MianalyzerBackend.GaussianMixturesController

include("../TestSupport.jl")

@testset "GaussianMixturesController" begin
  @testset ".call" begin
    request = Dict([
                    ("records", [[100.0], [200.0], [300.0], [200.0]]),
                    ("components", 1)
                   ])

    actual = JSON.parse(IOBuffer(GaussianMixturesController.call(request)))

    @test isapprox(actual["means"][1][1], 200.0; atol=eps(Float32), rtol=0)
    @test isapprox(actual["stds"][1][1], 70.71067812; atol=eps(Float32), rtol=0)
    @test isapprox(actual["weights"][1][1], 1.0; atol=eps(Float32), rtol=0)
  end
end


text:
export NSGA2Evolution,fastNonDominatedSort!,dominates,crowdingDistanceAssignement!,NSGA2Generation,NSGA2Populate,NSGA2Evaluate

import Cambrian.populate, Cambrian.evaluate,  Cambrian.selection, Cambrian.generation
import Cambrian.isless

mutable struct NSGA2Evolution{T<:Individual} <: Cambrian.AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{T}
    fitness::Function
    type::DataType
    rank::Dict{UInt64,Int64}
    distance::Dict{UInt64,Float64}
    gen::Int
    offsprings::Dict{UInt64,Bool}
end

function isless(i2::Individual, i1::Individual)
    k = length(i1.fitness)
    dom = false
    for i in 1:k
        if i1.fitness[i] > i2.fitness[i]
            return true
        end
    end
    return false
end

function tournament(pop::Array{T},t_size::Int) where {T <: Individual}
    inds = shuffle!(collect(1:length(pop)))
    sort(pop[inds[1:t_size]])[end]
end

populate(e::NSGA2Evolution) = NSGA2Populate(e)
evaluate(e::NSGA2Evolution) = NSGA2Evaluate(e)
generation(e::NSGA2Evolution) = NSGA2Generation(e)

function NSGA2Evolution{T}(cfg::NamedTuple, fitness::Function;
                      logfile=string("logs/", cfg.id, ".csv")) where {T <: Individual}
    logger = CambrianLogger(logfile)
    population = Cambrian.initialize(T, cfg)
    rank=Dict{UInt64,Int64}()
    distance=Dict{UInt64,Float64}()
    offsprings=Dict{UInt64,Bool}()
    for x in population
        offsprings[objectid(x)]=true
    end
    NSGA2Evolution(cfg, logger, population, fitness,T,rank,distance, 0, offsprings)
end

function NSGA2Evolution{T}(cfg::NamedTuple, fitness::Function, population::Array{T};
                      logfile=string("logs/", cfg.id, ".csv")) where {T <: Individual}
    logger = CambrianLogger(logfile)
    rank=Dict{UInt64,Int64}()
    distance=Dict{UInt64,Float64}()
    offsprings=Dict{UInt64,Bool}()
    for x in population
        offsprings[objectid(x)]=true
    end
    NSGA2Evolution(cfg, logger, population, fitness,T,rank,distance, 0, offsprings)
end

function NSGA2Evaluate(e::NSGA2Evolution)
    for i in eachindex(e.population)
        if e.offsprings[objectid(e.population[i])]
            e.population[i].fitness[:] = e.fitness(e.population[i])[:]
            e.offsprings[objectid(e.population[i])]=false
        end
    end
end

function NSGA2Populate(e::NSGA2Evolution)
    T=typeof(e.population[1])
    Qt=Array{T}(undef,0)
    for ind in e.population
        push!(Qt,ind)
    end
    i=0
    while i < e.config.n_offsprings
        parents =  [selection(e.population) for i in 1:2]
        if  rand() < e.config.p_crossover
            child1,child2 = crossover(parents...)
        else
            child1,child2= copy(parents[1]),copy(parents[2])
        end
        if rand() < e.config.p_mutation
            child1 = mutate(child1)
            child2 = mutate(child2)
        end
        push!(Qt,child1)
        e.offsprings[objectid(child1)]=true
        i+=1
        if i < e.config.n_offsprings
            push!(Qt,child2)
            e.offsprings[objectid(child2)]=true
            i+=1
        end
    end
    e.rank=Dict(objectid(x)=>0 for x in Qt)
    e.distance=Dict(objectid(x)=>0. for x in Qt)
    e.population=Qt
end

function dominates(e::NSGA2Evolution,ind1::T,ind2::T) where {T <: Individual}
    dom=false
    for i in 1:e.config.d_fitness
        if ind1.fitness[i]<ind2.fitness[i]
            return false
        elseif ind1.fitness[i]>ind2.fitness[i]
            dom=true
        end
    end
    return dom
end

function fastNonDominatedSort!(e::NSGA2Evolution)
    T=typeof(e.population[1])
    Fi=Array{T}(undef,0)

    n=Dict(objectid(x)=>0 for x in e.population)
    S=Dict(objectid(x)=>Array{T}(undef,0) for x in e.population)

    for ind1 in e.population
        for ind2 in e.population
            if dominates(e,ind1,ind2)
                push!(S[objectid(ind1)],ind2)
            elseif dominates(e,ind2,ind1)
                n[objectid(ind1)]+=1
            end
        end
        if n[objectid(ind1)]==0
            e.rank[objectid(ind1)]=1
            push!(Fi,ind1)
        end
    end

    i=1

    while isempty(Fi)==false
        Q=Array{T}(undef,0)
        for ind1 in Fi
            currentS=S[objectid(ind1)]
            for ind2 in currentS
                n[objectid(ind2)]-=1
                if n[objectid(ind2)]==0
                    e.rank[objectid(ind2)]=i+1
                    push!(Q,ind2)
                end
            end
        end
        i=i+1
        Fi=Q
    end
end

function crowdingDistanceAssignement!(e::NSGA2Evolution,I::Array{T}) where {T <: Individual}
    for x in I
        e.distance[objectid(x)]=0
    end
    l=length(I)
    for i in 1:e.config.d_fitness
        sort!(I,by=x->x.fitness[i])
        if I[1].fitness[i]!=I[end].fitness[i]
            e.distance[objectid(I[1])]=Inf
            e.distance[objectid(I[end])]=Inf
            quot=I[end].fitness[i]-I[1].fitness[i]
            for j in 2:l-1
                e.distance[objectid(I[j])]=e.distance[objectid(I[j])]+
                (I[j+1].fitness[i]-I[j-1].fitness[i])/quot
            end
        end
    end
end

function NSGA2Generation(e::NSGA2Evolution)
    if e.gen>1
        T=typeof(e.population[1])
        fastNonDominatedSort!(e)
        Pt1=Array{T}(undef,0)
        i=1
        sort!(e.population,by= x -> e.rank[objectid(x)])
        rank=1
        indIni=1
        indNext=findlast(x -> e.rank[objectid(x)] == rank , e.population)
        while indNext < e.config.n_population
            Pt1=[Pt1...,e.population[indIni:indNext]...]
            rank+=1
            indIni=indNext+1
            indNext=findlast(x -> e.rank[objectid(x)] == rank, e.population)
        end
        if isempty(Pt1)
            I=e.population[1:indNext]
            crowdingDistanceAssignement!(e,I)
            sort!(I, by= x->e.distance[objectid(x)],rev=true)
            Pt1=I[1:e.config.n_population]
        else
            I=e.population[indIni:indNext]
            crowdingDistanceAssignement!(e,I)
            sort!(I, by= x->e.distance[objectid(x)],rev=true)
            Pt1=[Pt1...,I[1:e.config.n_population-length(Pt1)]...]
        end
        e.population=Pt1
        e.offsprings=Dict(objectid(x)=>false for x in e.population)
    end
end


text:
using GZip

file_path = ARGS[1]

f = gzopen(joinpath(file_path))

for i=1:11
    print(readline(f))
end

close(f)

max_stars_repo_name: benhamner/avazu-ctr-prediction

max_issues_repo_name: benhamner/avazu-ctr-prediction

max_forks_repo_name: benhamner/avazu-ctr-prediction


text:
export ft_createdeviceinfolist

function ft_createdeviceinfolist()
  lpdwNumDevs = Ref{Cuint}()
  ft_status = ccall((:FT_CreateDeviceInfoList, d2xx),
                      Cuint,
                      (Ref{Cuint},),
                      lpdwNumDevs)
  checkstatus(ft_status)
  return lpdwNumDevs[]
end


text:
using Images, ImageView

function draw_seam(img, seam)
    img_w_seam = copy(img)
    for i = 1:size(img)[1]
        img_w_seam[i,seam[i]] = RGB(1,1,1)
    end

    return img_w_seam
end

function write_image(img, i; filebase = "out")
    save(filebase*lpad(string(i),5,string(0))*".png", img)
end

# function to return magnitude of image elements
function brightness(img_element)
    return img_element.r + img_element.g + img_element.b
end

# function to find energy of image
function find_energy(img)
    energy_x = imfilter(brightness.(img), Kernel.sobel()[2])
    energy_y = imfilter(brightness.(img), Kernel.sobel()[1])
    return sqrt.(energy_x.^2 + energy_y.^2)
    return abs2.(imfilter(brightness.(img), Kernel.sobel()))
end

# function to find direction to move in
function find_energy_map(energy)
    energy_map = zeros(size(energy))
    energy_map[end,:] .= energy[end,:]
    next_elements = zeros(Int, size(energy))

    for i = size(energy)[1]-1:-1:1
        for j = 1:size(energy)[2]
            left = max(1, j-1)
            right = min(size(energy)[2], j+1)
            if left < 1
                println("less than 1")
            end

            local_energy, next_element = findmin(energy_map[i+1, left:right])
            energy_map[i,j] += local_energy + energy[i,j]
            next_elements[i,j] = next_element -2

            # correct for only having 2 options on left edge.
            if left == 1 && size(energy)[2] > 3
                next_elements[i,j] += 1
            end
        end
    end

    return energy_map, next_elements
end

function generate_seam(energy, next_elements, element)
    seam = zeros(Int, size(next_elements)[1])
    seam[1] = element

    seam_energy = energy[element]

    for i = 2:length(seam)
        seam[i] = seam[i-1] + next_elements[i, seam[i-1]]
        seam_energy += energy[i, seam[i]]
    end

    return seam, seam_energy
end

# function to create seams and return seam of minimum energy
# n is number of seams
function find_seam(energy)

    energy_map, next_elements = find_energy_map(energy)

    kept_energy = Inf
    kept_seam = zeros(size(energy)[1])

    for i = 1:size(energy)[2]
        seam, seam_energy = generate_seam(energy, next_elements, i)
        if seam_energy < kept_energy
            kept_seam = seam
            kept_energy = seam_energy
        end
    end

    return kept_seam, kept_energy
end

# function to remove seams
function remove_seam(img, seam)
    img_res = (size(img)[1], size(img)[2]-1)

    # preallocate image
    new_img = Array{RGB}(undef, img_res)

    for i = 1:length(seam)
        if i > 1 && i < size(img)[2]
            new_img[i, :] .= vcat(img[i, 1:seam[i]-1],
                                  img[i, seam[i]+1:end])
        elseif i == 1
            new_img[i,:] .= img[i,2:end]
        elseif i == size(img)[2]
            new_img[i,:] .= img[i,1:end-1]
        end
    end

    return new_img
end

# Putting it all together
function seam_carving(img, res)

    if res < 0 || res > size(img)[2]
        error("resolution not acceptable")
    end

    for i = (1:size(img)[2] - res)
        energy = find_energy(img)
        seam, energy = find_seam(energy)
        img = remove_seam(img, seam)
        write_image(img, res, i)
    end
end


text:
struct SurrogateTransistor{T} <: TransistorModel
    W
    L
    T
    BaseModel::T
    Wrange
    Lrange
    Trange
    Voltagerange
    Idf
    gmf
    gdsf
end

function SurrogateTransistor(Wrange,Lrange,Trange,VoltageRange,model::TransistorModel)
    N_samples = 1000
    sampling_algorithm = LatinHypercubeSample()
    surrogate = LinearSurrogate
    opt = SRBF
    W = model.W
    L = model.L
    T = model.T0
    ## Create Ids Surrotage
    function fId(x)
        setW!(x[1],model)
        setL!(x[2],model)
        setT!(x[3],model)
        Id(x[4],x[5],x[6],model)
    end
    function fgm(x)
        setW!(x[1],model)
        setL!(x[2],model)
        setT!(x[3],model)
        gm(x[4],x[5],x[6],model)
    end
    function fgds(x)
        setW!(x[1],model)
        setL!(x[2],model)
        setT!(x[3],model)
        gds(x[4],x[5],x[6],model)
    end
    lb = [Wrange[1],Lrange[1],Trange[1],VoltageRange[1],VoltageRange[1],VoltageRange[1]]
    ub = [Wrange[2],Lrange[2],Trange[2],VoltageRange[2],VoltageRange[2],VoltageRange[2]]
    sx = sample(N_samples,lb,ub,sampling_algorithm)

    sId = fId.(sx)
    suId = surrogate(sx,sId,lb,ub)
    surrogate_optimize(fId,opt(),lb,ub,suId,sampling_algorithm;maxiters=N_samples)

    sgm = fgm.(sx)
    sugm = surrogate(sx,sgm,lb,ub)
    surrogate_optimize(fgm,opt(),lb,ub,sugm,sampling_algorithm;maxiters=N_samples)

    sgds = fId.(sx)
    sugds = surrogate(sx,sgds,lb,ub)
    surrogate_optimize(fgds,opt(),lb,ub,sugds,sampling_algorithm;maxiters=N_samples)
    SurrogateTransistor(W,L,model.T0,model,Wrange,Lrange,Trange,VoltageRange,suId,sugm,sugds)

end

Id(Vg,Vd,Vs,mod::SurrogateTransistor) = mod.Idf([mod.W,mod.L,mod.T,Vg,Vd,Vs])


text:
using SuiteSparseGraphBLAS
function pagerank(
    A,
    d = reduce(Monoids.PLUS_MONOID, A; dims=2),
    α = 0.85,
    maxiters = 100,
    ϵ = 1.0e-4
)
    n = size(A, 1)
    r = GBVector{Float32}(n)
    t = GBVector{Float32}(n)
    d[:, accum=BinaryOps.DIV] = α
    r[:] = 1.0 / n
    teleport = (1 - α) / n
    rdiff = 1.0
    i = 0
    for outer i ∈ 1:maxiters
        temp = t; t = r; r = temp
        w = t ./ d
        r[:] = teleport
        mul!(r, A', w, Semirings.PLUS_SECOND, accum=BinaryOps.PLUS)
        t .-= r
        map!(UnaryOps.ABS, t)
        rdiff = reduce(Monoids.PLUS_MONOID, t)
        if rdiff <= ϵ
            break
        end
    end
    return r, i
end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 2
coefs = [1;0.5;0.25]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 2, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def, m_step=2)  #don't combine time steps
st_single = NewStencilInstance(st_def, m_step=false)

## Input Data size Definition
radius = 2
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)

t_steps = 10
## GPU Implementation
gpu_out = ApplyStencil(st_inst, data, t_steps)
gpu_out = ApplyStencil(st_single, data, t_steps)

CUDA.cuProfilerStart()
gpu_out = ApplyStencil(st_inst, data, t_steps)
gpu_out = ApplyStencil(st_single, data, t_steps)

exit()
## CPU result
cpu_out = cpu_star_stencil(data, radius, coefs, t_steps)

## Check
mean_er = sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out)))
heatmap(gpu_out[:,:,Int(round(dz/2))+shift]-cpu_out[:,:,Int(round(dz/2))+shift])

@assert(mean_er/max(cpu_out...) < 0.000001)
## Visualize Error
using Plots

shift = 0
Int(round(dz/2))+shift
heatmap(gpu_out[:,:,Int(round(dz/2))+shift]-cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])

function init_coeff_arr(cf)
    ex = :(cfarr = SA_F32[])
    for i in cf
        push!(ex.args[2].args , i)
    end
    return ex
end

ex = :(cfarr = SA_F32[1])

using StaticArrays

a = SA_F32[1,2,3]

function ch_ker(arr)
    a = SA_F32[1,2,3,0.5]
    for i = 1:4
        arr[threadIdx().x] += a[i]
    end
    return nothing
end

ar = CUDA.ones(Float32, 32)
@cuda threads=32 ch_ker(ar)


text:
using Makie

function plot_curve()
    x = range(0, stop = 2π, length = 120)
    f(x) = sin.(x)
    y = f(x)
    scene = lines(x, y, color = :blue)
end

function plot_multiple()
    x = range(0, stop = 2π, length = 120)
    f1(x) = sin.(x)
    f2(x) = cos.(x)
    scene = lines(x, f2, color = :blue)
    lines!(scene, x, f1, color = :red, markersize = 0.1)
end

function plot_3dim()
    t = range(0,stop=2π,length=120)
    x = cos.(t)
    y = sin.(t)
    z = t
    lines(x,y,z)
end

function main()
    plot_3dim()
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module DataEnvelopmentAnalysis

    """
        DataEnvelopmentAnalysis
    A Julia package for efficiency and productivity measurement using Data Envelopment Analysis (DEA).
    [DataEnvelopmentAnalysis repository](https://github.com/javierbarbero/DataEnvelopmentAnalysis.jl).
    """

    using JuMP
    using GLPK
    using Ipopt
    using SparseArrays
    using LinearAlgebra

    using Printf: @sprintf
    using Statistics: std
    using StatsBase: CoefTable

    import StatsBase: nobs

    export
    # Types
    AbstractTechnicalDEAModel, AbstractRadialDEAModel,
    TechnicalDEAModel,
    RadialDEAModel, AdditiveDEAModel, DirectionalDEAModel, GeneralizedDFDEAModel,
    AbstractEconomicDEAModel,
    AbstractCostDEAModel, AbstractRevenueDEAModel, AbstractProfitDEAModel, AbstractProfitabilityDEAModel,
    CostDEAModel, RevenueDEAModel, ProfitDEAModel, ProfitabilityDEAModel,
    AbstractProductivityDEAModel,
    MalmquistDEAModel,
    # Technical models
    dea, deaadd, deaaddweights, deaddf, deagdf,
    efficiency, slacks,
    nobs, ninputs, noutputs, peers,
    # Economic models
    deacost, dearevenue, deaprofit, deaprofitability,
    # Productivity models
    malmquist,
    nperiods, prodchange

    # Include code of functions
    include("technical.jl")
    include("dea.jl")
    include("deaadd.jl")
    include("deaddf.jl")
    include("deagdf.jl")

    include("economic.jl")
    include("deacost.jl")
    include("dearevenue.jl")
    include("deaprofit.jl")
    include("deaprofitability.jl")

    include("productivity.jl")
    include("malmquist.jl")

    function __init__()

        nothing
    end

end # module


text:
using Test, LinearAlgebra

@testset "test/integration/facet_degeneracies.jl" begin

using BellScenario: Degeneracy

@testset "31->13 polytope" begin
    α = (3,1)
    β = (1,3)

    facets = [
        [0 -1 0 0 0 0 0], [0 0 -1 0 0 0 0], [0 0 0 -1 0 0 0],
        [0 0 0 0 -1 0 0], [0 0 0 0 0 -1 0], [0 0 0 0 0 0 -1],
        [-1 0 0 1 0 0 1], [-1 0 1 0 0 1 0], [-1 1 0 0 1 0 0],
        [-1 -1 0 1 -1 1 0], [-1 -1 1 0 -1 0 1], [-1 0 -1 1 1 -1 0],
        [-1 0 1 -1 1 0 -1], [-1 1 -1 0 0 -1 1], [-1 1 0 -1 0 1 -1]
    ]

    out_sym1 = facets[[1,4,9]]
    out_sym2 = facets[[2,5,8]]
    out_sym3 = facets[[3,6,7]]
    out_sym4 = facets[10:15]

    in_sym1 = facets[1:9]
    in_sym2 = facets[10:15]

    input_relabels = Degeneracy.bipartite_input_relabels(α,β,"normalized")
    output_relabels = Degeneracy.bipartite_output_relabels(α,β,"normalized")

    @testset "verify output relabeling symmetry groups" for sym_group in [out_sym1, out_sym2, out_sym3, out_sym4]
        for i in [1,2]
            for relabel in output_relabels[i]
                for facet in sym_group
                    @test in(facet*relabel,sym_group)
                end
            end
        end
    end

    @testset "verify symmetry groups 1,2,3 are closed under input relabelings" for sym_group in [in_sym1, in_sym2]
        for relabel in input_relabels
            for facet in sym_group
                @test in(facet*relabel,sym_group)
            end
        end
    end
end

end


text:
using ConvDiffMIPDECO
using Test
using jInv.Mesh

println("=== 2D test === ")

domain = [0. 3. 0 1. ]
n      = 2*[7 9 ]
M      = getRegularMesh(domain,n)
bc = (:dir,:neu,:neu,:neu)
gd = X->sin.(X[:,2]*pi)
pFor = getConvDiffFEMParam(M,bc=bc,gd=gd)

A,B,b = getConvDiffFEMConstraintsAMPL(M,gd=gd)

mc = rand(M.nc);
dobs = getData(mc,pFor)
utrue = pFor.Fields;
utest = A\(B*mc+b); # -->    A *u = B*m + b
@test norm(utrue-utest)/norm(utrue) < 1e-10

println("=== 3D test === ")

domain = [0. 3. 0 1. 0 2.]
n      = 2*[7 9 12]
M      = getRegularMesh(domain,n)
bc = (:dir,:neu,:neu,:neu,:neu,:neu)
gd = X->sin.(X[:,2]*pi)
pFor = getConvDiffFEMParam(M,bc=bc,gd=gd)

A,B,b = getConvDiffFEMConstraintsAMPL(M,gd=gd)

mc = rand(M.nc);
dobs = getData(mc,pFor)
utrue = pFor.Fields;
utest = A\(B*mc+b);
@test norm(utrue-utest)/norm(utrue) < 1e-10

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
function load(fname::AbstractString,model,channel)
    @info("fname $fname")

    df = DataFrame(CSV.File(fname, types = Dict(:channel => Int)))

    if isempty(df)
        return df
    else
        df.time = DateTime.(df.time,dateformat"yyyy-mm-dd HH:MM:SS")

        dfT = df[(df.model .== model) .& (df.channel .== channel),:]
        @show size(dfT)
        return dfT
    end
end

function load(fnames,model,channel)
    return reduce(vcat,load.(fnames,model,channel))
end

function loadavg(fnames::AbstractVector{<:AbstractString},parameter,taverage;
                 timerange = (typemin(DateTime),typemax(DateTime)))
    if parameter == "temperature"
        model = "Hideki-TS04"
        channel = 2
        param = "temperature_C"
    elseif parameter == "humidity"
        model = "Hideki-TS04"
        channel = 2
        param = "humidity"
    elseif parameter == "wind_speed"
        model = "Hideki-Wind"
        param = "wind_avg_mi_h"
        channel = 4
    elseif parameter == "wind_direction"
        model = "Hideki-Wind"
        param = "wind_dir_deg"
        channel = 4
    elseif parameter == "rain"
        model = "Hideki-Rain"
        param = "rain_mm"
        channel = 4
    end

    dfT = load(fnames,model,channel)
    dfT[!,:time] = Dates.epochms2datetime.(taverage * round.(Int64,Dates.datetime2epochms.(dfT[:,:time]) / taverage));
    df2 = by(dfT,:time; Dict(Symbol(parameter) => (Symbol(param) => mean))...)

    filter!(row -> timerange[1] <= row.time <= timerange[2],df2)

    return df2
end

function filetimerange(fname)
    time_start = DateTime(replace(replace(basename(fname),"out-" => ""),".csv"  => ""),dateformat"yyyy-mm-ddTHH:MM:SS")
    # file can stop earlier
    time_end = time_start + Dates.Day(1)
    return (time_start,time_end)
end

function listfile(basedir,timerange)
    fnames = String[]

    for fname in sort(glob("out-*csv",basedir))
        time_start,time_end = filetimerange(fname)
        if (time_start <= timerange[2]) & (timerange[1] <= time_end)
            @info "keep $fname"
            push!(fnames,fname)
        else
            @info "ignore $fname"
        end
    end

    return fnames
end

function loadavg(basedir::AbstractString,parameter,taverage;
                 timerange = (typemin(DateTime),typemax(DateTime)),
                 kwargs...)

    Δt = Dates.Millisecond(taverage)
    extendedtimerange = (timerange[1]-Δt,timerange[2]+Δt)

    fnames = listfile(basedir,extendedtimerange)
    return loadavg(fnames,parameter,taverage; timerange = extendedtimerange, kwargs...)
end

function combine(; basedir = get(ENV,"WEATHERSTATION_DIR","/var/lib/WeatherStation.jl/"))

taverage = 1000 * 60 * 60
taverage = 1000 * 60 * 10

channel = 3
channel = 2
model = "Hideki-Wind"
model = "Hideki-TS04"

filename = joinpath(basedir,"dataframe.csv")

dfT = load(fnames,model,channel)

end


text:
module TestDeprecated

using Test, DataFrames
using DataStructures: OrderedDict, SortedDict

const ≅ = isequal

@testset "by and aggregate" begin
    @test_throws ArgumentError by()
    @test_throws ArgumentError aggregate()
end

@testset "deprecated broadcasting assignment" begin
    df = DataFrame(a=1:4, b=1, c=2)
    df.a .= 'a':'d'
    @test df == DataFrame(a=97:100, b=1, c=2)
    dfv = view(df, 2:3, 2:3)
    dfv.b .= 0
    @test df.b == [1, 0, 0, 1]
end

@testset "All indexing" begin
    df = DataFrame(a=1, b=2, c=3)

    @test select(df, All(1, 2)) == df[:, 1:2]
    @test select(df, All(1, :b)) == df[:, 1:2]
    @test select(df, All(:a, 2)) == df[:, 1:2]
    @test select(df, All(:a, :b)) == df[:, 1:2]
    @test select(df, All(2, 1)) == df[:, [2, 1]]
    @test select(df, All(:b, 1)) == df[:, [2, 1]]
    @test select(df, All(2, :a)) == df[:, [2, 1]]
    @test select(df, All(:b, :a)) == df[:, [2, 1]]

    @test df[:, All(1, 2)] == df[:, 1:2]
    @test df[:, All(1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2)] == df[:, 1:2]
    @test df[:, All(:a, :b)] == df[:, 1:2]
    @test df[:, All(2, 1)] == df[:, [2, 1]]
    @test df[:, All(:b, 1)] == df[:, [2, 1]]
    @test df[:, All(2, :a)] == df[:, [2, 1]]
    @test df[:, All(:b, :a)] == df[:, [2, 1]]

    @test df[:, All(1, 1, 2)] == df[:, 1:2]
    @test df[:, All(:a, 1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2, :b)] == df[:, 1:2]
    @test df[:, All(:a, :b, 2)] == df[:, 1:2]
    @test df[:, All(2, 1, :a)] == df[:, [2, 1]]

    @test select(df, All(1, "b")) == df[:, 1:2]
    @test select(df, All("a", 2)) == df[:, 1:2]
    @test select(df, All("a", "b")) == df[:, 1:2]
    @test select(df, All("b", 1)) == df[:, [2, 1]]
    @test select(df, All(2, "a")) == df[:, [2, 1]]
    @test select(df, All("b", "a")) == df[:, [2, 1]]

    @test df[:, All(1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2)] == df[:, 1:2]
    @test df[:, All("a", "b")] == df[:, 1:2]
    @test df[:, All("b", 1)] == df[:, [2, 1]]
    @test df[:, All(2, "a")] == df[:, [2, 1]]
    @test df[:, All("b", "a")] == df[:, [2, 1]]

    @test df[:, All("a", 1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2, "b")] == df[:, 1:2]
    @test df[:, All("a", "b", 2)] == df[:, 1:2]
    @test df[:, All(2, 1, "a")] == df[:, [2, 1]]

    df = DataFrame(a1=1, a2=2, b1=3, b2=4)
    @test df[:, All(r"a", Not(r"1"))] == df[:, [1, 2, 4]]
    @test df[:, All(Not(r"1"), r"a")] == df[:, [2, 4, 1]]
end

@testset "indicator in joins" begin
    name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
    job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

    @test outerjoin(name, job, on = :ID, indicator=:source) ≅
          outerjoin(name, job, on = :ID, source=:source)
    @test leftjoin(name, job, on = :ID, indicator=:source) ≅
          leftjoin(name, job, on = :ID, source=:source)
    @test rightjoin(name, job, on = :ID, indicator=:source) ≅
          rightjoin(name, job, on = :ID, source=:source)

    @test_throws ArgumentError outerjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
    @test_throws ArgumentError leftjoin(name, job, on = :ID,
                                       indicator=:source, source=:source)
    @test_throws ArgumentError rightjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
end

@testset "map on GroupedDataFrame" begin
    df = DataFrame(a=1:3, b=4:6, c=7:9)
    dfv = @view df[1:3, 1:3]
    gdf = groupby(df, :a)
    gdfv = groupby(dfv, :a)

    for x in (gdf, gdfv)
        @test collect(x) == map(identity, x)
    end
end

@testset "new map behavior" begin
    df = DataFrame(g=[1, 2, 3])
    gdf = groupby(df, :g)
    @test map(nrow, gdf) == [1, 1, 1]
end

@testset "Conversion tests" begin
    df = DataFrame()
    df[!, :A] = 1:5
    df[!, :B] = [:A, :B, :C, :D, :E]
    @test isa(convert(Matrix, df), Matrix{Any})
    @test isa(convert(Matrix{Any}, df), Matrix{Any})
    @test isa(convert(Array, df), Matrix{Any})
    @test isa(convert(Array{Any}, df), Matrix{Any})

    df = DataFrame()
    df[!, :A] = 1:5
    df[!, :B] = 1.0:5.0
    @test isa(convert(Matrix, df), Matrix{Float64})
    @test isa(convert(Matrix{Any}, df), Matrix{Any})
    @test isa(convert(Matrix{Float64}, df), Matrix{Float64})
    @test isa(convert(Array, df), Matrix{Float64})
    @test isa(convert(Array{Any}, df), Matrix{Any})
    @test isa(convert(Array{Float64}, df), Matrix{Float64})

    df = DataFrame()
    df[!, :A] = Vector{Union{Float64, Missing}}(1.0:5.0)
    df[!, :B] = Vector{Union{Float64, Missing}}(1.0:5.0)
    a = convert(Matrix, df)
    aa = convert(Matrix{Any}, df)
    ai = convert(Matrix{Int}, df)
    @test isa(a, Matrix{Union{Float64, Missing}})
    @test a == convert(Array, convert(Matrix{Union{Float64, Missing}}, df))
    @test a == convert(Matrix, df)
    @test a == Matrix(df)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Matrix{Any}, df)
    @test aa == Matrix{Any}(df)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Matrix{Int}, df)
    @test ai == Matrix{Int}(df)

    df[1, 1] = missing
    @test_throws ArgumentError convert(Matrix{Float64}, df)
    na = convert(Matrix{Union{Float64, Missing}}, df)
    naa = convert(Matrix{Any}, df)
    nai = convert(Matrix{Union{Int, Missing}}, df)
    @test isa(na, Matrix{Union{Float64, Missing}})
    @test na ≅ convert(Matrix, df)
    @test na ≅ Matrix(df)
    @test isa(naa, Matrix{Union{Any, Missing}})
    @test naa ≅ convert(Matrix{Any}, df)
    @test naa ≅ Matrix{Any}(df)
    @test isa(nai, Matrix{Union{Int, Missing}})
    @test nai ≅ convert(Matrix{Union{Int, Missing}}, df)
    @test nai ≅ Matrix{Union{Int, Missing}}(df)

    df = DataFrame()
    df[!, :A] = Vector{Union{Float64, Missing}}(1.0:5.0)
    df[!, :B] = Vector{Union{Float64, Missing}}(1.0:5.0)
    a = convert(Array, df)
    aa = convert(Array{Any}, df)
    ai = convert(Array{Int}, df)
    @test isa(a, Matrix{Union{Float64, Missing}})
    @test a == convert(Array, convert(Array{Union{Float64, Missing}}, df))
    @test a == convert(Array, df)
    @test a == Array(df)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Array{Any}, df)
    @test aa == Array{Any}(df)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Array{Int}, df)
    @test ai == Array{Int}(df)

    df[1, 1] = missing
    @test_throws ArgumentError convert(Array{Float64}, df)
    na = convert(Array{Union{Float64, Missing}}, df)
    naa = convert(Array{Any}, df)
    nai = convert(Array{Union{Int, Missing}}, df)
    @test isa(na, Matrix{Union{Float64, Missing}})
    @test na ≅ convert(Array, df)
    @test na ≅ Array(df)
    @test isa(naa, Matrix{Union{Any, Missing}})
    @test naa ≅ convert(Array{Any}, df)
    @test naa ≅ Array{Any}(df)
    @test isa(nai, Matrix{Union{Int, Missing}})
    @test nai ≅ convert(Array{Union{Int, Missing}}, df)
    @test nai ≅ Array{Union{Int, Missing}}(df)

    a = Union{Float64, Missing}[1.0, 2.0]
    b = Union{Float64, Missing}[-0.1, 3]
    c = Union{Float64, Missing}[-3.1, 7]
    di = Dict("a"=>a, "b"=>b, "c"=>c)

    df = convert(DataFrame, di)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in sort(collect(keys(di)))]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame, od)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in keys(od)]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame, sd)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in keys(sd)]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    a = [1.0]
    di = Dict("a"=>a, "b"=>b, "c"=>c)
    @test_throws DimensionMismatch convert(DataFrame, di)
end

@testset "deprecated conversion for DataFrameRow and GroupKey" begin
    df = DataFrame(a=1)
    dfr = df[1, :]
    key = keys(groupby(df, :a))[1]

    @test convert(Vector, dfr) == Vector(dfr) == [1]
    @test convert(Vector{Any}, dfr) == Vector{Any}(dfr) == Any[1]
    @test convert(Array, dfr) == Vector(dfr) == [1]
    @test convert(Array{Any}, dfr) == Vector{Any}(dfr) == Any[1]
    @test convert(Tuple, dfr) == Tuple(dfr) == (1,)
    @test convert(Vector, key) == Vector(key) == [1]
    @test convert(Vector{Any}, key) == Vector{Any}(key) == Any[1]
    @test convert(Array, key) == Vector(key) == [1]
    @test convert(Array{Any}, key) == Vector{Any}(key) == Any[1]
    @test convert(Tuple, key) == Tuple(key) == (1,)

    @test convert(Vector, dfr) isa Vector{Int}
    @test convert(Vector{Any}, dfr) isa Vector{Any}
    @test convert(Array, dfr) isa Vector{Int}
    @test convert(Array{Any}, dfr) isa Vector{Any}
    @test convert(Tuple, dfr) isa Tuple{Int}
    @test convert(Vector, key) isa Vector{Int}
    @test convert(Vector{Any}, key) isa Vector{Any}
    @test convert(Array, key) isa Vector{Int}
    @test convert(Array{Any}, key) isa Vector{Any}
    @test convert(Tuple, key) isa Tuple{Int}
end

@testset "DataFrameRow convert" begin
    df = DataFrame(a=[1, missing, missing], b=[2.0, 3.0, 0.0])
    dfr = DataFrameRow(df, 1, :)
    @test convert(Vector, dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
    @test convert(Vector{Int}, dfr)::Vector{Int} == [1, 2]
    @test convert(Array, dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
    @test convert(Array{Int}, dfr)::Vector{Int} == [1, 2]

    dfr = DataFrame(a=1, b=2)[1, :]
    dfr2 = DataFrame(c=3, d=4)[1, :]
    @test convert(Tuple, dfr) == (1, 2)
end

end # module


text:
#=
    Access standard AbstractNode internals without using .std syntax
=#
@inline function Base.getproperty(c::AbstractNode, s::Symbol)
    if s in (
        :id,
        :lb,
        :ub,
    )
        Core.getproperty(Core.getproperty(c, :std), s)
    else
        getfield(c, s)
    end
end

@inline function Base.setproperty!(c::AbstractNode, s::Symbol, v)
    if s in (
        :id,
        :lb,
        :ub,
    )
        Core.setproperty!(c.std, s, v)
    else
        Core.setproperty!(c, s, v)
    end
end

"""
    is_approx_feasible(tree::BnBTree, value)

Return whether a given `value` is approximately feasible based on the tolerances defined in the tree options.
"""
function is_approx_feasible(tree::BnBTree, value::Number)
    return is_approx_feasible(value; atol=tree.options.atol, rtol=tree.options.rtol)
end

function is_approx_feasible(value::Number; atol=1e-6, rtol=1e-6)
    return isapprox(value, round(value); atol, rtol)
end

"""
    get_distance_to_feasible(tree::BnBTree, value)

Return the distance of feasibility for the given value.

- if `value::Number` this returns the distance to the nearest discrete value
"""
function get_distance_to_feasible(tree::BnBTree, value::Number)
    return abs(round(value)-value)
end


text:
Base.MainInclude.include( "../fileFind/filesfind.jl")
Base.MainInclude.include( "../RawFileofflineProc\\BaseLib.jl")
function getWindStatus!(WindStatusChangeLogfile,windinfoAofA1)
    # tm0 = Base.Filesystem.mtime(WindStatusChangeLogfile )+8*60*60
    # dt0 = Dates.unix2datetime( tm0)
    # jzdt = Dates.DateTime(baseDatetime ,Dates.DateFormat("yyyymmddHHMMSS" ))

    # relt= floor(UInt,(dt0-jzdt).value/1000 )
    # println( join( (dt0,jzdt,relt),"; ") )
    # println( dt0 ,WindStatusChangeLogfile)
    fiobj = open(WindStatusChangeLogfile,"r")
    skip(fiobj,0)
    relt=0
    filename = basename( WindStatusChangeLogfile)
    # println( filename )
    for line in readlines( fiobj )
        #line :67170634,ID=5991:6238,EventFirsttime:67170603,L0time：2019.03.30 10.30.34, Winding: 1
        # println( line)
        a = split( line,",")
        if length(a)  == 1
            time1 =replace(a[1],Pair("CurDatatime=","") )
            # println(time1)
            relt = parse(UInt,time1)
            dt0= UInt2DTStr( relt )
            append!(windinfoAofA1,[[ relt,0,0 ,false , relt ,filename,dt0]])
            continue
        end

        # for a0 in a
            ID =replace(a[2],Pair("ID=","") )

            segID1str ,segID2str= split(ID,":")

            segID1 =parse(UInt,segID1str )

            segID2 =parse(UInt,segID2str )
            # println(ID,segID2str)
            a3 = replace(a[3],"EventFirsttime:"=>"")
            a3 = replace(a3,"EventFirsttime："=>"")

            t0 = parse(UInt,a3)
            isWind =  occursin("1" , a[5]) #? true : false
            # println( join( (a[2],a[3],a[5] ,segID1,segID2 ,t0,isWind),";"))
            relt = parse( UInt,a[1])
            dt0= UInt2DTStr( relt )
            append!(windinfoAofA1,[[ t0,segID1,segID2 ,isWind , relt ,filename,dt0]])
        # end
    end
    # return a
end

function findwindofftime( seg1,seg2,tLog, windinfoAofA2)
    logt0 = tLog
    for idx =1 : length(windinfoAofA2 )
        a0= windinfoAofA2[idx]
        t1,segid1,segid2,iswind2,tLog2 = a0[1],a0[2],a0[3],a0[4],a0[5]
        logt0  = tLog2
        if iswind2 == false  && segid1 == seg1 && segid2 ==seg2
            return t1
        end
    end #for
    println("                NOT find  the offtime of $tLog, return the  last time: $logt0  in: $tLog  ")
    return logt0
end

function CalwindPeriod( windinfoAofA1 )
    a=[]
    for idx =1 : length( windinfoAofA1)
        a0= windinfoAofA1[idx]
        t0,seg1,seg2,iswind,tLog,filename = a0[1],a0[2],a0[3],a0[4],a0[5],a0[6]
        if iswind == true
            # if 67531208 == t0
            #    println( join( (t0,seg1,seg2,iswind,tLog ,filename),",") )
            # end
            @views piecedata = windinfoAofA1[idx+1:end]
            toff = findwindofftime( seg1,seg2,tLog,piecedata )
            append!(a, [[ t0,toff ,seg1,seg2,iswind,tLog, ]] )
        end
    end
    return a
end

# curtime,SegRange,Winding,StartTime,EndTime,LastingTime,CurDateTime
# 67089378,0199_0496,N,67085141,67089347,4206,2019.03.29 11.56.18
function SaveWindPeriods(windPeriod ,filename::String )
    fobj =open(filename,"w")
    for idx =1:length( windPeriod)
        a0= windPeriod[idx]
        # println( length( a0 ),typeof(a0))
        t0,toff ,seg1,seg2,iswind,tLog  = a0[1],a0[2],a0[3],a0[4],a0[5],a0[6]
        seg1Str = string(seg1,base=10,pad=4 )
        Segs= join( (string(seg1,base=10,pad=4 ),string(seg2,base=10,pad=4 ) ),"_" )

        if iswind == true && toff > t0
            str1= join( (t0,Segs,"Y",t0,toff ,string(toff -t0,base=10,pad=8 ) ,UInt2DTStr( t0,"yyyy.mm.dd HH.MM.SS" )),",")

        else
            str1 = join( (t0,Segs,"N",t0,toff ,string(toff -t0,base=10,pad=8 ) ,UInt2DTStr( t0,"yyyy.mm.dd HH.MM.SS" )),",")
        end

        write( fobj, str1 *"\r\n")
    end
    close( fobj)
end

function Cal_SaveWindperiod(rootdir , windLogfile)
    windinfoAofA=[]
    files = findAllFiles(rootdir,".log","","" )
    sort!( files)
    for file1 in files
        println(basename(file1) )
    end
    for file1 in files
        getWindStatus!( file1,windinfoAofA )
    end
# for a in windinfoAofA
#     println( a )
# end
    windPeriod1 = CalwindPeriod( windinfoAofA )
#     println("windows ")
# for a in windPeriod1
#     for b in a
#       print(b,", ")
#     end
#     println()
# end

    SaveWindPeriods(windPeriod1 ,windLogfile )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function InterpolationTableLoader()
    return InterpolationTableLoader(())
end

function get_abscissa_grid(obj::InterpolationTableLoader)
    return jcall(obj, "getAbscissaGrid", Vector{jdouble}, ())
end

function get_ordinate_grid(obj::InterpolationTableLoader)
    return jcall(obj, "getOrdinateGrid", Vector{jdouble}, ())
end

function get_values_samples(obj::InterpolationTableLoader)
    return jcall(obj, "getValuesSamples", Vector{Vector{jdouble}}, ())
end

function load_data(obj::InterpolationTableLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::InterpolationTableLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((2,),),
    ]
insts["fast"] = [
    ((5,),),
    ((10,),),
    ]
insts["various"] = [
    ((5,),),
    ((10,),),
    ((15,),),
    ]
return (NearestCorrelationJuMP, insts)


text:
module PlutoHooks

include("./notebook.jl")

export @use_state, @use_effect, @use_memo, @use_ref, @use_deps
export @use_is_pluto_cell, @only_as_script, @skip_as_script

end


text:
using Abaco
using Test

interval=10
ts = nowts()

sn1 = "tn01"
sn2 = "tn02"
sn3 = "tn03"

val1 = 100
val2 = 101
val3 = 102
x = 10

expected_triggers = 12
actual_triggers = 0

function onresult(ts, ne, name, value, inputs)
    global actual_triggers
    @debug "timestamp [$ts]: ne: [$ne] $name = $value"
    if name === "offset_footprint"
        @test value == x + val1 + val2 + val3
    end
    actual_triggers += 1
end

function onresult(timestamp, ne, name, value::Abaco.PValue, inputs)
    global actual_triggers
    @debug "progressive timestamp [$ts]: ne: [$ne] $name = $value"
    if name === "footprint_sum"
        if value.contribs == 1
            @test ne == "trento"
            @test timestamp == span(ts, interval)
            @test value.value == val1
        elseif value.contribs == 2
            @test value.value == val1 + val2
        elseif value.contribs == 3
            @test value.value == val1 + val2 + val3
        end
    end
    if name === "footprint_mean"
        if value.contribs == 1
            @test ne == "trento"
            @test timestamp == span(ts, interval)
            @test value.value == val1
        elseif value.contribs == 2
            @test value.value == (val1 + val2)/2
        elseif value.contribs == 3
            @test value.value == (val1 + val2 + val3)/3
        end
    end
    actual_triggers += 1
end

abaco = abaco_init(onresult, interval=interval)

@debug "starting demo with time span: $interval, values timestamp: $ts"

setup_settings(abaco, "hub", oncomplete=onresult)

formula(abaco, "footprint_sum", "sum(sensor.footprint)", "hub")
formula(abaco, "footprint_mean", "mean(sensor.footprint)", "hub")
formula(abaco, "footprint_std", "std(sensor.footprint)", "hub")
formula(abaco, "progressive_fake", "x", "hub")
formula(abaco, "foobar", "x+y", "hub")
formula(abaco, "offset_footprint", "x + footprint_sum", "hub")

city = node(abaco, "trento", "hub")

node(abaco, city, sn1, "sensor")
node(abaco, city, sn2, "sensor")
node(abaco, city, sn3, "sensor")

current_footprint = get_collected(abaco, city.ne, "sensor.footprint")
@debug "initial carbon footprints: $current_footprint"

ingest(abaco, ts, sn1, Dict("footprint" => val1))
ingest(abaco, ts, sn2, Dict("footprint" => val2))

current_footprint = get_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "current carbon footprint: $current_footprint"

footprint_sum = sum_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "footprint sum: $footprint_sum"
@test footprint_sum.contribs == 2
@test footprint_sum.expected == 3
@test footprint_sum.value == val1 + val2

ingest(abaco, ts, sn3, Dict("footprint" => val3))

footprint_sum = sum_collected(abaco, city.ne, "sensor.footprint", ts)
@debug "footprint sum: $footprint_sum"
@test footprint_sum.contribs == 3
@test footprint_sum.value == val1 + val2 + val3

ingest(abaco, ts, city.ne, Dict(
    "x" => x,
    "y" => 20
))

@test actual_triggers === expected_triggers


text:
#=
For details, please visit:
https://computationalmindset.com/en/neural-networks/experiments-with-neural-odes-in-julia.html#exp3
=#

using Plots
using Flux
using DiffEqFlux

tbegin = 0.0
tend = 4.0
datasize = 51
dataset_in = range(tbegin, tend, length=datasize)
dataset_out = log.(1 .+ dataset_in) .- 3 * sqrt.(dataset_in)

function neural_network(data_dim)
    fc = FastChain(FastDense(data_dim, 64, swish),
                  FastDense(64, 32, swish),
                  FastDense(32, data_dim))
end

nn = neural_network(1)
theta = initial_params(nn)

predict(t, p) = nn(t', p)'

loss(p) = begin
  yhat = predict(dataset_in, p)
  l = Flux.mse(yhat, dataset_out)
end

learning_rate=1e-2
opt = ADAMW(learning_rate)
epochs = 500

function cb_train(theta, loss)
    println("Loss: ", loss)
    false
end

res_train = DiffEqFlux.sciml_train(
    loss, theta, opt,
    maxiters = epochs,
    cb = cb_train)

y_pred = predict(dataset_in, res_train.minimizer)

pl = plot(
    dataset_in,
    dataset_out,
    linewidth=2, ls=:dash,
    title="MLP by FastChain without Neural ODEs",
    xaxis="t",
    label="original y(t)",
    legend=:topright)
display(pl)

pl = plot!(
    dataset_in,
    y_pred,
    linewidth=1,
    label="predicted y(t)")
display(pl)

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
# This file was generated, do not modify it. # hide
fit!(y1_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y1_oos(), ls="none", marker="o", label="linear oos")
legend()

savefig(joinpath(@OUTPUT, "s2.svg")) # hide


text:
module LassoPlot

export plot

using Lasso, Plots, StatsPlots, DataFrames, MLBase

include("plots.jl")

end # module


text:
function run()
    E = Ensamble()
    run(E)
end
function run(filename::String)
    E = Ensamble(filename)
    run(E)
end

function run(E::Ensamble)
    #glfw_window = to_native(display(E.scene))

    #btn = Button(E.scene, label="Selection")

    #on(btn.clicks) do
    #    E.mode = SelectionMode()
    #end

    # counts = Node([1, 4, 3, 7, 2])

    # fig[2, 1] = buttongrid = GridLayout(tellwidth = false)

    # buttonlabels = [@lift("Count: $($counts[i])") for i in 1:5]

    # buttons = buttongrid[1, 1:5] = [Button(fig, label = l) for l in buttonlabels]

    on(E.scene.events.keyboardbutton) do event
        respond_keyboard_event(E, E.mode, event)
    end

    on(events(E.scene).mousebutton, priority=Int8(20)) do event
        respond_mouse_event(E, E.mode, event)
    end

    E.scene
end


text:
using Continuables
import BenchmarkTools.@benchmark

@cont function crange(n::Int)
  for i in 1:n
    cont(i)
  end
end

function trange(n::Int)
  c = Channel{Int}(1)
  task = @async for i ∈ 1:n
    put!(c, i)
  end
  bind(c, task)
end

@Ref function sum_continuable(continuable)
  a = Ref(0)
  continuable() do i
    a += i
  end
  a
end

function sum_continuable_withoutref(continuable)
  a = 0
  continuable() do i
    a += i
  end
  a
end

function sum_iterable(it)
  a = 0
  for i in it
    a += i
  end
  a
end

function collect_continuable(continuable)
  a = []
  continuable() do i
    push!(a, i)
  end
  a
end

function collect_iterable(it)
  a = []
  for i in it
    push!(a, i)
  end
  a
end

@benchmark sum_continuable(crange(1000))
#=
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.185 ns (0.00% GC)
  median time:      1.580 ns (0.00% GC)
  mean time:        1.756 ns (0.00% GC)
  maximum time:     57.679 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
=#

@benchmark sum_continuable(@i2c 1:1000)
#=
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.185 ns (0.00% GC)
  median time:      1.580 ns (0.00% GC)
  mean time:        1.877 ns (0.00% GC)
  maximum time:     31.210 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
=#

@benchmark sum_continuable_withoutref(crange(1000))
#=
BenchmarkTools.Trial:
  memory estimate:  22.81 KiB
  allocs estimate:  1460
  --------------
  minimum time:     26.074 μs (0.00% GC)
  median time:      27.654 μs (0.00% GC)
  mean time:        39.603 μs (15.68% GC)
  maximum time:     42.758 ms (99.85% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#

@benchmark sum_continuable_withoutref(@i2c 1:1000)
#=
BenchmarkTools.Trial:
  memory estimate:  22.81 KiB
  allocs estimate:  1460
  --------------
  minimum time:     26.074 μs (0.00% GC)
  median time:      27.654 μs (0.00% GC)
  mean time:        39.910 μs (16.12% GC)
  maximum time:     46.847 ms (99.92% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#

@benchmark sum_iterable(trange(1000))
#=
BenchmarkTools.Trial:
  memory estimate:  33.05 KiB
  allocs estimate:  2019
  --------------
  minimum time:     4.968 ms (0.00% GC)
  median time:      5.839 ms (0.00% GC)
  mean time:        6.089 ms (0.34% GC)
  maximum time:     23.466 ms (73.27% GC)
  --------------
  samples:          821
  evals/sample:     1
=#

@benchmark sum_iterable(1:1000)
#=
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.185 ns (0.00% GC)
  median time:      1.580 ns (0.00% GC)
  mean time:        1.742 ns (0.00% GC)
  maximum time:     56.494 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
=#

@benchmark collect_continuable(crange(1000))
#=
BenchmarkTools.Trial:
  memory estimate:  24.03 KiB
  allocs estimate:  499
  --------------
  minimum time:     11.456 μs (0.00% GC)
  median time:      12.247 μs (0.00% GC)
  mean time:        20.871 μs (30.69% GC)
  maximum time:     43.985 ms (99.97% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#

@benchmark collect_continuable(@i2c 1:1000)
#=
BenchmarkTools.Trial:
  memory estimate:  24.03 KiB
  allocs estimate:  499
  --------------
  minimum time:     11.456 μs (0.00% GC)
  median time:      12.247 μs (0.00% GC)
  mean time:        25.120 μs (26.38% GC)
  maximum time:     44.245 ms (99.96% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#

@benchmark collect_iterable(trange(1000))
#=
BenchmarkTools.Trial:
  memory estimate:  57.08 KiB
  allocs estimate:  2518
  --------------
  minimum time:     4.968 ms (0.00% GC)
  median time:      5.715 ms (0.00% GC)
  mean time:        6.082 ms (1.31% GC)
  maximum time:     57.567 ms (89.97% GC)
  --------------
  samples:          822
  evals/sample:     1
=#

@benchmark collect_iterable(1:1000)
#=
BenchmarkTools.Trial:
  memory estimate:  24.03 KiB
  allocs estimate:  499
  --------------
  minimum time:     11.456 μs (0.00% GC)
  median time:      12.642 μs (0.00% GC)
  mean time:        23.706 μs (27.20% GC)
  maximum time:     45.308 ms (99.97% GC)
  --------------
  samples:          10000
  evals/sample:     1
=#


text:
function Base.broadcast!(::typeof(identity), ta::Tuple{A1,A2}, tb::Tuple{B1,B2}) where {A1,A2,B1,B2}
    ta1, ta2 = ta
    tb1, tb2 = tb
    @assert indices(ta1) == indices(ta2) == indices(tb1) == indices(tb2)
    for i in eachindex(ta1)
        @inbounds ta1[i], ta2[i] = tb1[i], tb2[i]
    end
    ta
end

function Base.broadcast!(ta::Tuple{A1,A2}, tb::Tuple{B1,B2}) where {A1,A2,B1,B2}
    ta1, ta2 = ta
    tb1, tb2 = tb
    @assert indices(ta1) == indices(ta2) == indices(tb1) == indices(tb2)
    for i in eachindex(ta1)
        @inbounds ta1[i], ta2[i] = tb1[i], tb2[i]
    end
    ta
end

function Base.broadcast!(ta::Tuple{A1,A2}, tb::Tuple{B1,B2}) where {A1,A2,B1,B2}
    ta1, ta2 = ta
    tb1, tb2 = tb

    @assert indices(ta1) == indices(ta2) == indices(tb1) == indices(tb2)
    for i in eachindex(ta1)
        @inbounds ta1[i], ta2[i] = tb1[i], tb2[i]
    end
    ta
end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
"""
`module  JAC.Einstein`
    ... a submodel of JAC that contains all methods for computing Einstein A and B coefficients;
        these computations are treated as level properties with regard to a single multiplet.
"""
module Einstein

    using Printf, ..AngularMomentum, ..Basics, ..Defaults, ..InteractionStrength, ..ManyElectron, ..PhotoEmission, ..Radial, ..TableStrings

    """
    `struct  Einstein.Settings  <:  AbstractPropertySettings`  ... defines a type for the details and parameters of computing Einstein lines and coefficients

        + multipoles              ::Array{EmMultipoles}          ... Specifies the multipoles of the radiation field that are to be included.
        + printBefore             ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection           ::LineSelection                ... Specifies the selected levels, if any.
        + photonEnergyShift       ::Float64                      ... An overall energy shift for all photon energies.
        + mimimumPhotonEnergy     ::Float64                      ... minimum transition energy for which (photon) transitions are included into the
                                                                     computation.
        + maximumPhotonEnergy     ::Float64                      ... maximum transition energy for which (photon) transitions are included.
    """
    struct Settings  <:  AbstractPropertySettings
        multipoles                ::Array{EmMultipole,1}
        printBefore               ::Bool
        lineSelection             ::LineSelection
        photonEnergyShift         ::Float64
        mimimumPhotonEnergy       ::Float64
        maximumPhotonEnergy       ::Float64
    end

    """
    `Einstein.Settings(settings::Einstein.Settings;`

                multipoles=..,            printBefore=..,      lineSelection=..,
                photonEnergyShift=..,     mimimumPhotonEnergy=..,         maximumPhotonEnergy=..)

        ... constructor for re-defining the settings::Einstein.Settings.
    """
    function Settings(settings::Einstein.Settings;                            multipoles::Union{Nothing,Array{EmMultipole,1}}=nothing,
        printBefore::Union{Nothing,Bool}=nothing,                             lineSelection::Union{Nothing,LineSelection}=nothing,
        photonEnergyShift::Union{Nothing,Float64}=nothing,
        mimimumPhotonEnergy::Union{Nothing,Float64}=nothing,                  maximumPhotonEnergy::Union{Nothing,Float64}=nothing)

        if  multipoles         == nothing   multipolesx          = settings.multipoles            else  multipolesx = multipoles                   end
        if  printBefore        == nothing   printBeforex         = settings.printBefore           else  printBeforex = printBefore                 end
        if  lineSelection      == nothing   lineSelectionx       = set.lineSelection              else  lineSelectionx = lineSelection             end
        if  photonEnergyShift  == nothing   photonEnergyShiftx   = settings.photonEnergyShift     else  photonEnergyShiftx = photonEnergyShift     end
        if  mimimumPhotonEnergy== nothing   mimimumPhotonEnergyx = settings.mimimumPhotonEnergy   else  mimimumPhotonEnergyx = mimimumPhotonEnergy end
        if  maximumPhotonEnergy== nothing   maximumPhotonEnergyx = settings.maximumPhotonEnergy   else  maximumPhotonEnergyx = maximumPhotonEnergy end

    	Settings(multipolesx, printBeforex, selectLinesx, selectedLinesx, photonEnergyShiftx, mimimumPhotonEnergyx, maximumPhotonEnergyx)
    end

    """
    `Einstein.Settings()`  ... constructor for the default values of Einstein line computations.
    """
    function Settings()
        Settings(EmMultipole[], false, LineSelection(), 0., 0., 0.)
    end

    #== """
    `Einstein.Settings(multipoles::Array{EmMultipole,1};` printBefore::Bool=false,
                       selectLines::Bool=false, selectedLines::Array{Tuple{Int64,Int64},1}=Tuple{Int64,Int64}[],
                       photonEnergyShift::Float64=0., mimimumPhotonEnergy::Float64=0., maximumPhotonEnergy::Float64=0.)
        ... keyword constructor to overwrite selected value of Einstein line computations; DON'T use it anymore (11/2019).
    """
    function Settings(multipoles::Array{EmMultipole,1}; printBefore::Bool=false,
                      selectLines::Bool=false, selectedLines::Array{Tuple{Int64,Int64},1}=Tuple{Int64,Int64}[],
                      photonEnergyShift::Float64=0., mimimumPhotonEnergy::Float64=0., maximumPhotonEnergy::Float64=0.)
        println("!!! UPDATE: Remove this call and use Einstein.Settings(settings::Einstein.Settings) instead.")
        Settings(multipoles, printBefore, selectLines, selectedLines, photonEnergyShift, mimimumPhotonEnergy, maximumPhotonEnergy)
    end  ==#

    # `Base.show(io::IO, settings::Einstein.Settings)`  ... prepares a proper printout of the variable settings::Einstein.Settings.
    function Base.show(io::IO, settings::Einstein.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "lineSelection:            $(settings.lineSelection)  ")
        println(io, "photonEnergyShift:        $(settings.photonEnergyShift)  ")
        println(io, "mimimumPhotonEnergy:      $(settings.mimimumPhotonEnergy)  ")
        println(io, "maximumPhotonEnergy:      $(settings.maximumPhotonEnergy)  ")
    end

    """
    `struct  Einstein.Channel`  ... defines a type for a single Einstein emission/absorption channel that specifies the multipole, gauge and amplitude.

        + multipole         ::EmMultipole        ... Multipole of the photon emission/absorption.
        + gauge             ::EmGauge            ... Gauge for dealing with the (coupled) radiation field.
        + amplitude         ::Complex{Float64}   ... Amplitude of this multiple channel.
    """
    struct Channel
        multipole           ::EmMultipole
        gauge               ::EmGauge
        amplitude           ::Complex{Float64}
    end

    # `Base.show(io::IO, channel::Einstein.Channel)`  ... prepares a proper printout of the variable channel::EinsteinChannel.
    function Base.show(io::IO, channel::Einstein.Channel)
        print(io, "$(channel.multipole) [$(channel.gauge); $(channel.amplitude)")
    end

    """
    `struct  Einstein.Line`
        ... defines a type for a Einstein line that may include the definition of sublines and their corresponding amplitudes.

        + initialLevel   ::Level            ... initial-(state) level
        + finalLevel     ::Level            ... final-(state) level
        + omega          ::Float64          ... photon energy = transition energy
        + photonRate     ::EmProperty       ... Total rate of this line.
        + angularBeta    ::EmProperty       ... Angular beta_2 coefficient.
        + hasSublines    ::Bool             ... Determines whether the individual sublines are defined in terms of
                                                their multipolarity, amplitude, or not.
        + channels       ::Array{Einstein.Channel,1}   ... List of radiative (photon) channels.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        omega            ::Float64
        photonRate       ::EmProperty
        angularBeta      ::EmProperty
        hasSublines      ::Bool
        channels         ::Array{Einstein.Channel,1}
   end

    """
    `Einstein.Line(initialLevel::Level, finalLevel::Level, photonRate::Float64)`
        ... constructor an Einstein line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, photonRate::EmProperty)
        Line(initialLevel, finalLevel, 0, 0, photonRate, EmProperty(0., 0.), false, Einstein.Channel[])
   end

    # `Base.show(io::IO, line::Einstein.Line)`  ... prepares a proper printout of the variable line::Einstein.Line.
    function Base.show(io::IO, line::Einstein.Line)
        println(io, "initialLevel:         $(line.initialLevel)  ")
        println(io, "finalLevel:           $(line.finalLevel)  ")
        println(io, "omega:                $(line.omega)  ")
        println(io, "photonRate:           $(line.photonRate)  ")
        println(io, "angularBeta:          $(line.angularBeta)  ")
        println(io, "hasSublines:          $(line.hasSublines)  ")
        println(io, "channels:             $(line.channels)  ")
    end

    """
    `Einstein.computeLines(multiplet::Multiplet, grid::Radial.Grid, settings::Einstein.Settings; output::Bool=true)`
        ... to compute the Einstein transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{Einstein.Lines} is returned.
    """
    function  computeLines(multiplet::Multiplet, grid::Radial.Grid, settings::Einstein.Settings; output::Bool=true)
        println("")
        printstyled("Einstein.computeLines(): The computation of the transition amplitudes and properties starts now ... \n", color=:light_green)
        printstyled("--------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        lines = Einstein.determineLines(multiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    Einstein.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = Einstein.Line[]
        for  line in lines
            newLine = Einstein.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        Einstein.displayRates(stdout, newLines)
        Einstein.displayLifetimes(stdout, newLines)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    Einstein.displayRates(iostream, newLines);    Einstein.displayLifetimes(iostream, newLines)   end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `Einstein.computeAmplitudesProperties(line::Einstein.Line, grid::Radial.Grid, settings::Einstein.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::Einstein.Line is returned for which the
            amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::Einstein.Line, grid::Radial.Grid, settings::Einstein.Settings)
        global JAC_counter
        newChannels = Einstein.Channel[];    rateC = 0.;    rateB = 0.
        for channel in line.channels
            #
            amplitude = PhotoEmission.amplitude("emission", channel.multipole, channel.gauge, line.omega, line.finalLevel, line.initialLevel, grid)
            #
            push!( newChannels, Einstein.Channel( channel.multipole, channel.gauge, amplitude) )
            if       channel.gauge == Basics.Coulomb     rateC = rateC + abs(amplitude)^2
            elseif   channel.gauge == Basics.Babushkin   rateB = rateB + abs(amplitude)^2
            elseif   channel.gauge == Basics.Magnetic    rateB = rateB + abs(amplitude)^2;   rateC = rateC + abs(amplitude)^2
            end
        end
        # Calculate the photonrate and angular beta if requested
        wa = 8.0pi * Defaults.getDefaults("alpha") * line.omega / (Basics.twice(line.initialLevel.J) + 1) *
                                                      (Basics.twice(line.finalLevel.J) + 1)
        photonrate  = EmProperty(wa * rateC, wa * rateB)
        angularBeta = EmProperty(-9., -9.)
        line = Einstein.Line( line.initialLevel, line.finalLevel, line.omega, photonrate, angularBeta, true, newChannels)
        return( line )
    end

    """
    `Einstein.determineChannels(finalLevel::Level, initialLevel::Level, settings::Einstein.Settings)`
        ... to determine a list of Einstein.Channel for a transitions from the initial to final level and by taking into
            account the particular settings of for this computation; an Array{Einstein.Channel,1} is returned.
    """
    function  determineChannels(finalLevel::Level, initialLevel::Level, settings::Einstein.Settings)
        channels = Einstein.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        for  mp in settings.multipoles
            if   AngularMomentum.isAllowedMultipole(symi, mp, symf)
                if     string(mp)[1] == 'E'   push!(channels, Einstein.Channel(mp, Basics.Coulomb,   0.) )
                                              push!(channels, Einstein.Channel(mp, Basics.Babushkin, 0.) )
                elseif string(mp)[1] == 'M'   push!(channels, Einstein.Channel(mp, Basics.Magnetic,  0.) )
                else   error("stop a")
                end
            end
        end
        return( channels )
    end

    """
    `Einstein.determineLines(multiplet::Multiplet, settings::Einstein.Settings)`
        ... to determine a list of EinsteinLine's for transitions between the lines from the given multiplet and by taking
            into account the particular selections and settings of for this computation; an Array{Einstein.Line,1} is
            returned. Apart from the level specification, all physical properties are set to zero during the initialization
            process.
    """
    function  determineLines(multiplet::Multiplet, settings::Einstein.Settings)
        lines = Einstein.Line[]
        for  iLevel  in  multiplet.levels
            for  fLevel  in  multiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    omega    = iLevel.energy - fLevel.energy + settings.photonEnergyShift
                    if  omega <= settings.mimimumPhotonEnergy  ||  omega > settings.maximumPhotonEnergy    continue   end
                    channels = determineChannels(fLevel, iLevel, settings)
                    push!( lines, Einstein.Line(iLevel, fLevel, omega, EmProperty(0., 0.), EmProperty(0., 0.), true, channels) )
                end
            end
        end
        return( lines )
    end

    #== """
    `Einstein.determineLines(multiplet::Multiplet, settings::Einstein.Settings)`
        ... to determine a list of EinsteinLine's for transitions between the lines from the given multiplet and by taking
            into account the particular selections and settings of for this computation; an Array{Einstein.Line,1} is
            returned. Apart from the level specification, all physical properties are set to zero during the initialization
            process.
    """
    function  determineLines(multiplet::Multiplet, settings::Einstein.Settings)
        if    settings.selectLines    selectLines   = true
            selectedLines = Basics.determineSelectedLines(settings.selectedLines, multiplet, multiplet)
        else                          selectLines   = false
        end

        lines = Einstein.Line[]
        for  i = 1:length(multiplet.levels)
            for  f = 1:length(multiplet.levels)
                if  selectLines  &&  !((i,f) in selectedLines )    continue   end
                omega    = multiplet.levels[i].energy - multiplet.levels[f].energy + settings.photonEnergyShift
                if  omega <= settings.mimimumPhotonEnergy  ||  omega > settings.maximumPhotonEnergy    continue   end

                channels = determineChannels(multiplet.levels[f], multiplet.levels[i], settings)
                push!( lines, Einstein.Line(multiplet.levels[i], multiplet.levels[f], omega, EmProperty(0., 0.), EmProperty(0., 0.),
                                            true, channels) )
            end
        end
        return( lines )
    end ==#

    """
    `Einstein.displayLifetimes(stream::IO, lines::Array{Einstein.Line,1})`
        ... to list all lifetimes of the selected lines as can be derived from the calculated data. A neat table
            is printed but nothing is returned otherwise.
    """
    function  displayLifetimes(stream::IO, lines::Array{Einstein.Line,1})
        # Determine all initial levels (and their level information) for printing the lifetimes
        ilevels = Int64[];   istr = String[]
        for  i = 1:length(lines)
            ii = lines[i].initialLevel.index;
            if  !(ii in ilevels)
               sa  = "  ";    sym = LevelSymmetry( lines[i].initialLevel.J, lines[i].initialLevel.parity )
               sa = sa * TableStrings.center(10, TableStrings.level(lines[i].initialLevel.index); na=2)
               sa = sa * TableStrings.center(10, string(sym); na=4)
               sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", lines[i].initialLevel.energy)) * "    "
               push!( ilevels, ii);    push!( istr, sa )
            end
        end
        # Determine the lifetime (in a.u.) of the selected initial levels
        irates = Basics.EmProperty[]
        for  ii in  ilevels
            waCoulomb = 0.;    waBabushkin = 0.
            for  i = 1:length(lines)
                if   lines[i].initialLevel.index == ii
                    waCoulomb   = waCoulomb   + lines[i].photonRate.Coulomb
                    waBabushkin = waBabushkin + lines[i].photonRate.Babushkin
                end
            end
            push!(irates, EmProperty( waCoulomb, waBabushkin) )
        end

        nx = 108
        println(stream, " ")
        println(stream, "  Radiative lifetimes (as derived from these computations):")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                              sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                              sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(12, "Level energy"   ; na=4);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(12, "Used Gauge"    ; na=6);                     sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(26, "Lifetime"; na=2);
        sb = sb * TableStrings.center(26, "[a.u.]"*"          "*TableStrings.inUnits("time"); na=5)
        sa = sa * TableStrings.center(12, "Decay widths"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  ii = 1:length(ilevels)
            sa = istr[ii]
            sa = sa * "Coulomb          " * @sprintf("%.6e",              1.0/irates[ii].Coulomb)               * "  "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("time: from atomic",   1.0/irates[ii].Coulomb) )   * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic",     irates[ii].Coulomb) )
            println(stream, sa)
            sa = repeat(" ", length(istr[ii]) )
            sa = sa * "Babushkin        " * @sprintf("%.6e",              1.0/irates[ii].Babushkin)             * "  "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("time: from atomic",   1.0/irates[ii].Babushkin) ) * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic",     irates[ii].Babushkin) )
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `Einstein.displayLines(lines::Array{Einstein.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of
            all selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{Einstein.Line,1})
        nx = 115
        println(" ")
        println("  Selected Einstein lines:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(30, "List of multipoles"; na=4);             sb = sb * TableStrings.hBlank(34)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.omega)) * "    "
            mpGaugeList = Tuple{Basics.EmMultipole,Basics.EmGauge}[]
            for  i in 1:length(line.channels)
                push!( mpGaugeList, (line.channels[i].multipole, line.channels[i].gauge) )
            end
            sa = sa * TableStrings.multipoleGaugeTupels(50, mpGaugeList)
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        println(" ")
        #
        return( nothing )
    end

    """
    `Einstein.displayRates(stream::IO, lines::Array{Einstein.Line,1})`
        ... to list all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayRates(stream::IO, lines::Array{Einstein.Line,1})
        nx = 145
        println(stream, " ")
        println(stream, "  Einstein coefficients, transition rates and oscillator strengths:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "Energy"   ; na=4);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center( 9, "Multipole"; na=1);                         sb = sb * TableStrings.hBlank(10)
        sa = sa * TableStrings.center(11, "Gauge"    ; na=4);                         sb = sb * TableStrings.hBlank(17)
        sa = sa * TableStrings.center(26, "A--Einstein--B"; na=2);
        sb = sb * TableStrings.center(26, TableStrings.inUnits("rate")*"           "*TableStrings.inUnits("rate"); na=2)
        sa = sa * TableStrings.center(11, "Osc. strength"    ; na=3);                 sb = sb * TableStrings.hBlank(17)
        sa = sa * TableStrings.center(12, "Decay widths"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            for  ch in line.channels
                sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                               fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
                sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
                sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.omega)) * "    "
                sa = sa * TableStrings.center(9,  string(ch.multipole); na=4)
                sa = sa * TableStrings.flushleft(11, string(ch.gauge);  na=2)
                chRate =  8pi * Defaults.getDefaults("alpha") * line.omega / (Basics.twice(line.initialLevel.J) + 1) * (abs(ch.amplitude)^2) *
                                                                 (Basics.twice(line.finalLevel.J) + 1)
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to Einstein A",  line, chRate)) * "  "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to Einstein B",  line, chRate)) * "    "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to g_f",         line, chRate)) * "    "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to decay width", line, chRate)) * "    "
                println(stream, sa)
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractDomain

An abstract type to aid with custom spatial domain types. Users
can subtype their domains from this type, and implement the methods
in `geotraits/domain.jl`.
"""
abstract type AbstractDomain end

"""
    view(domain, inds)

Return a view of `domain` at given indices `inds`.
"""
Base.view(domain::AbstractDomain, inds::AbstractVector{Int}) =
  DomainView(domain, inds)

#------------------
# IMPLEMENTATIONS
#------------------
"""
    SpatialDomain{T,N}

Spatial domain in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type SpatialDomain{T,N} <: AbstractDomain end

ncoords(::SpatialDomain{T,N}) where {T,N} = N
coordtype(::SpatialDomain{T,N}) where {T,N} = T

include("domains/point_set.jl")
include("domains/regular_grid.jl")
include("domains/structured_grid.jl")

"""
    DomainView(domain, inds)

Return a view of `domain` at `inds`.
"""
struct DomainView{𝒟,I} <: AbstractDomain
  domain::𝒟
  inds::I
end

nelms(dv::DomainView) = length(dv.inds)
ncoords(dv::DomainView) = ncoords(dv.domain)
coordtype(dv::DomainView) = coordtype(dv.domain)
coordinates!(buff::AbstractVector, dv::DomainView, ind::Int) =
  coordinates!(buff, dv.domain, dv.inds[ind])

==(dv₁::DomainView, dv₂::DomainView) =
  dv₁.domain == dv₂.domain && dv₁.inds == dv₂.inds

Base.collect(dv::DomainView) = PointSet(coordinates(dv))

function Base.show(io::IO, dv::DomainView)
  N = ncoords(dv)
  T = coordtype(dv)
  npts = nelms(dv)
  print(io, "$npts DomainView{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", dv::DomainView)
  println(io, dv)
  Base.print_array(io, coordinates(dv))
end


text:
function anyₛ(xs::AbstractArray{T, N}) where {T <: SoftBool, N}
  SoftBool(maximum([sb.logerr for sb in xs]))
end

function allₛ(xs::AbstractArray{T, N}) where {T <: SoftBool, N}
  SoftBool(sum([sb.logerr for sb in xs]))
end

function anyₛ(xs::AbstractArray{T, N}) where {T <: DualSoftBool, N}
  b0 = allₛ([x.b0 for x in xs])
  b1 = anyₛ([x.b1 for x in xs])
  DualSoftBool(b0, b1)
end

function allₛ(xs::AbstractArray{T, N}) where {T <: DualSoftBool, N}
  b0 = anyₛ([x.b0 for x in xs])
  b1 = allₛ([x.b1 for x in xs])
  DualSoftBool(b0, b1)
end

softeq(a::AbstractArray, b::AbstractArray) = allₛ(softeq.(a, b))