text:
struct StreamerRegistry
	streamers::Dict{MIME,Vector{FormatHandler}}
	favorites::Dict{MIME,FormatHandler}
	global_favorites::Vector{FormatHandler}
end

StreamerRegistry() = StreamerRegistry(
		Dict{MIME,Vector{FormatHandler}}(),
		Dict{MIME,FormatHandler}(),
		Vector{FormatHandler}())

Base.copy(registry::StreamerRegistry) = StreamerRegistry(
		copy(registry.streamers),
		copy(registry.favorites),
		copy(registry.global_favorites))

function Base.empty!(registry::StreamerRegistry)
	empty!(registry.streamers)
	empty!(registry.favorites)
	empty!(registry.global_favorites)
end

function Base.merge!(registry::StreamerRegistry, other::StreamerRegistry)
	merge!(registry.streamers, other.streamers)
	merge!(registry.favorites, other.favorites)
	append!(registry.global_favorites, other.global_favorites)
end

const registry = StreamerRegistry()

function newregistry(f::Function)
	tmp_registry = copy(registry)
	empty!(registry)
	f()
	merge!(registry, tmp_registry)
end

"""
	addstreamer(format::AbstractString, streamer::FormatHandler)

Register a `streamer` for the specified `format`.

# Example

```julia
struct MyIO <: Formats.FormatHandler end
Formats.addstreamer("image/png", MyIO())
```
"""
function addstreamer(format::AbstractString, streamer::FormatHandler)
	mime = MIME{Symbol(format)}()
	streamers = get!(registry.streamers, mime, [])
	if streamer in streamers
		error("streamer $(streamer) already registered for $(format)")
	else
		push!(streamers, streamer)
	end
	if length(streamers) > 1
		@info("$(format) has multiple registered streamers")
	end
end

"""
	preferstreamer(streamer::FormatHandler, [format::AbstractString])

Register `streamer` as the preferred alternative for the specified `format`.

If `format` is omitted, `streamer` becomes a globally preferred alternative: it
is used for all the formats it can stream.
```
"""
function preferstreamer(streamer::FormatHandler)
	if streamer in registry.global_favorites
		error("streamer $(streamer) is already globally preferred")
	end
	push!(registry.global_favorites, streamer)
end

function preferstreamer(streamer::FormatHandler, format::AbstractString)
	mime = MIME{Symbol(format)}()
	streamers = get(registry.streamers, mime, [])
	if ! (streamer in streamers)
		error("streamer $(streamer) is not registered for $(format)")
	end
	if haskey(registry.favorites, mime)
		@warn("replacing preferred streamer for $(format)")
	end
	registry.favorites[mime] = streamer
end

resolvestreamer(mime::MIME) =
		Formats.resolvehandler(mime, "streamer", registry.streamers,
		registry.favorites, registry.global_favorites)


text:
export PoolingDesc, pooling!, pooling, ∇pooling!
export CUDNN_POOLING_MAX
export CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING
export CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING

type PoolingDesc{N}
  ptr::Ptr{Void}
  window::NTuple{N,Int}
  padding::NTuple{N,Int}
  stride::NTuple{N,Int}
end

function PoolingDesc{N}(mode, window::NTuple{N,Int}, padding::NTuple{N,Int}, stride::NTuple{N,Int};
  maxpooling_nanopt=CUDNN_NOT_PROPAGATE_NAN)
  p = Ptr{Void}[0]
  cudnnCreatePoolingDescriptor(p)
  cwindow = Cint[window[i] for i=N:-1:1]
  cpadding = Cint[padding[i] for i=N:-1:1]
  cstride = Cint[stride[i] for i=N:-1:1]
  cudnnSetPoolingNdDescriptor(p[1], mode, maxpooling_nanopt, N, cwindow, cpadding, cstride)
  desc = PoolingDesc(p[1], window, padding, stride)
  finalizer(desc, cudnnDestroyPoolingDescriptor)
  desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::PoolingDesc) = desc.ptr

function pooling!{T}(x::CuArray{T}, y::CuArray{T}, poolingdesc; alpha=1.0, beta=0.0)
  h = gethandle(x.dev)
  xdesc = TensorDesc(x)
  ydesc = TensorDesc(y)
  cudnnPoolingForward(h, poolingdesc, T[alpha], xdesc, x, T[beta], ydesc, y)
end

function pooling(x::CuArray, poolingdesc; alpha=1.0)
  N = length(poolingdesc.window)
  w, p, s = poolingdesc.window, poolingdesc.padding, poolingdesc.stride
  outdims = Array(Int, N)
  for i = 1:N
    outdims[i] = (size(x,i) + 2*p[i] - w[i]) ÷ s[i] + 1
  end
  y = similar(x, outdims..., size(x,N+1), size(x,N+2))
  pooling!(x, y, poolingdesc, alpha)
end

function ∇pooling!{T}(poolingdesc, y::CuArray{T}, dy::CuArray{T},
  x::CuArray{T}, dx::CuArray{T}; alpha=1.0, beta=0.0)
  h = gethandle(x.dev)
  ydesc = TensorDesc(y)
  dydesc = TensorDesc(dy)
  xdesc = TensorDesc(x)
  dxdesc = TensorDesc(dx)
  cudnnPoolingBackward(h, poolingdesc, ydesc, y, dydesc, dy, xdesc, x,
    T[beta], dxdesc, dx)
end

ext: jl


text:
using SparseMerkleTree
using BitIntegers
using Test

@testset "SparseMerkleTree.jl" begin
    db = new_tree()
    @test db.root == SparseMerkleTree.Hash(parse(UInt256, "0x876422b7697ae7c337e2ee7727feb3db474adf7be1cf04b6b5857d82d610e88a"))

    k = UInt256(200)
    v = Vector{UInt8}(b"foo")
    db[k] = v
    @test db.root == Hash(parse(UInt256, "0x108b52f48b1d7f361359def8c726f1b6256de9e2b41e363e9604ce0556e3fb2d"))

    p1 = merkle_proof(db, k)
    @test verify_proof(p1, db.root, k, v)

    k = UInt256(0xabababab000ffffccccc)
    v = Vector{UInt8}(b"bar")
    db[k] = v
    @test db.root == Hash(parse(UInt256, "0x973de6aec749609fa159cb96f53e2b130f9280faefa16bf61ef97915da4401f6"))

    p2 = merkle_proof(db, k)
    @test verify_proof(p2, db.root, k, v)

    p3 = compress_proof(p2)
    p4 = decompress_proof(p3)
    @test p2 == p4
    @test verify_proof(p4, db.root, k, v)
end


text:
export localips

const port = Ref{Int}()

@init port[] = haskey(ENV, "BLINK_PORT") ? parse(Int, ENV["BLINK_PORT"]) : rand(2_000:10_000)

const ippat = r"([0-9]+\.){3}[0-9]+"

@static if Sys.isunix()
    localips() = map(IPv4, readlines(`ifconfig` |>
                      `grep -Eo $("inet (addr:)?$(ippat.pattern)")` |>
                      `grep -Eo $(ippat.pattern)` |>
                      `grep -v $("127.0.0.1")`))
end

# Browser Window

@static if Sys.isapple()
    launch(x) = run(`open $x`)
elseif Sys.islinux()
    launch(x) = run(`xdg-open $x`)
elseif Sys.iswindows()
    launch(x) = run(`cmd /C start $x`)
end

localurl(p::Page) = "http://127.0.0.1:$(port[])/$(id(p))"

launch(p::Page) = (launch(localurl(p)); p)


text:
using QuantileRegression
using DataFrames, CSV, LinearAlgebra
using Test

include("engel_test.jl")


text:
using LightGraphs
using MetaGraphs
using Tweaks.GraphUtil

@testset "edge_weight" begin
    g = MetaGraph(5, 0.5)
    @test get_weight(g, Edge(1, 2), -11) == -11
    @test get_weight(g, 1, 2, -11) == -11

    add_edge!(g, 2, 3)
    @test get_weight(g, Edge(2, 3), -11) == 0.5
    @test get_weight(g, 2, 3, -11) == 0.5

    add_edge!(g, 3, 4)
    set_prop!(g, 3, 4, :weight, 34)
    @test get_weight(g, Edge(3, 4), -11) == 34
    @test get_weight(g, 3, 4, -11) == 34
end

@testset "vertex_weight" begin
    g = MetaGraph(5, 0.5)
    @test get_weight(g, 2) == 0.5
    set_prop!(g, 3, :weight, 0.3)
    @test get_weight(g, 3) == 0.3
end

@testset "set_weight" begin
    g = MetaGraph(5, 0.5)
    add_edge!(g, 2, 3)
    set_weight!(g, 2, 3, 23)
    @test get_prop(g, 2, 3, :weight) == 23
    add_edge!(g, 3, 4)
    set_weight!(g, Edge(3, 4), 34)
    @test get_prop(g, 3, 4, :weight) == 34
    set_weight!(g, 2, 0.2)
    @test get_prop(g, 2, :weight) == 0.2
end

@testset "add_edge" begin
    g = MetaGraph(5, 0.5)
    add_edge!(g, 2, 3, 23)
    @test get_prop(g, 2, 3, :weight) == 23
    add_edge!(g, Edge(3, 4), 34)
    @test get_prop(g, 3, 4, :weight) == 34
end

@testset "iterate_edge" begin
    @test collect(Edge(1, 2)) == [1, 2]
end

@testset "show_light_graph" begin
    buf = IOBuffer()
    show(buf, MIME("image/svg+xml"), Graph(1))
    seek(buf, 0)
    @test String(read(buf, 5)) == "<?xml"
end

@testset "show_meta_graph" begin
    buf = IOBuffer()
    show(buf, MIME("image/svg+xml"), MetaGraph(1))
    seek(buf, 0)
    @test String(read(buf, 5)) == "<?xml"
end


text:
using Bokeh

# more or less the simplest plotting example
x = linspace(0, 2pi)
y = sin(x)#[sin(x) cos(x)]
# by default autoopen is set to true so `plot` will try to open a browser
# with the plot file (which by default is "bokeh_plot.html")
Bokeh.debug(true)
# plot(y, "k--")
hold(true)
autoopen(false)
ys = [sin(x) cos(x)]
plot(ys, "ko|bo")
plot(sin)
showplot()


text:
"""
    abstract type Quadrature{D,T} <: GridapType end

-[`get_coordinates(q::Quadrature)`](@ref)
-[`get_weights(q::Quadrature)`](@ref)
-[`test_quadrature`](@ref)
"""
abstract type Quadrature{D,T} <: GridapType end

# Abstract API

"""
    get_coordinates(q::Quadrature)
"""
function get_coordinates(q::Quadrature)
  @abstractmethod
end

"""
    get_weights(q::Quadrature)
"""
function get_weights(q::Quadrature)
  @abstractmethod
end

# Default API

"""
    num_points(q::Quadrature)
"""
num_points(q::Quadrature) = length(get_weights(q))

"""
    num_point_dims(::Quadrature{D}) where D
    num_point_dims(::Type{<:Quadrature{D}}) where D
"""
num_point_dims(::Quadrature{D}) where D = D

num_point_dims(::Type{<:Quadrature{D}}) where D = D

"""
    num_dims(::Quadrature{D}) where D where D
    num_dims(::Type{<:Quadrature{D}}) where D
"""
num_dims(::Quadrature{D}) where D = D

num_dims(::Type{<:Quadrature{D}}) where D = D

# Tester

"""
    test_quadrature(q::Quadrature{D,T}) where {D,T}
"""
function test_quadrature(q::Quadrature{D,T}) where {D,T}
  x = get_coordinates(q)
  w = get_weights(q)
  @test isa(x,Vector{Point{D,T}})
  @test isa(w,Vector{T})
  @test length(x) == num_points(q)
  @test length(w) == num_points(q)
  @test D == num_dims(q)
  @test D == num_point_dims(q)
  @test D == num_dims(typeof(q))
  @test D == num_point_dims(typeof(q))
end

# Generic concrete implementation

"""
    struct GenericQuadrature{D,T} <: Quadrature{D,T}
      coordinates::Vector{Point{D,T}}
      weights::Vector{T}
    end
"""
struct GenericQuadrature{D,T} <: Quadrature{D,T}
  coordinates::Vector{Point{D,T}}
  weights::Vector{T}
end

get_coordinates(q::GenericQuadrature) = q.coordinates

get_weights(q::GenericQuadrature) = q.weights


text:
#---------------------------------------------------------------------#
#This function warps a 2D grid
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function warp_grid!(coord,Npoin,DFloat)

    for i=1:Npoin
        x, y = coord[1,i], coord[2,i]
        coord[1,i]=x + sin(π * x) * sin(2 * π * y) / 15
        coord[2,i]=y + sin(2 * π * x) * sin(π * y) / 15
    end

end


text:
using KUnet
using HDF5

@KUnet.useif CUDArt

if isdefined(:CUDArt)
const libkunet = find_library(["libkunet"], ["."])
typealias Mat CudaArray{Float32,2}
typealias Cmat Ptr{Float32}
end

if isdefined(:CUDArt)
function forward!(x::Matrix{Float32}, x1::Mat, w1::Mat, b1::Mat, x2::Mat, w2::Mat, b2::Mat, x3::Mat, y::Matrix{Float32})
    batch = 937
    xrows,xcols = size(x)
    yrows,ycols = size(y)
    for i=1:batch:xcols
        copy!(x1,1,x,(i-1)*xrows+1,length(x1))
        CUBLAS.gemm!('N','N',1.0f0,w1,x1,0.0f0,x2)
        ccall((:badd,libkunet),Void,(Cint,Cint,Cmat,Cmat),size(x2,1),size(x2,2),x2,b1)
        ccall((:reluforw,libkunet),Void,(Cint,Cmat),length(x2),x2)
        CUBLAS.gemm!('N','N',1.0f0,w2,x2,0.0f0,x3)
        ccall((:badd,libkunet),Void,(Cint,Cint,Cmat,Cmat),size(x3,1),size(x3,2),x3,b2)
        copy!(y,(i-1)*yrows+1,x3,1,length(x3))
    end
    y
end
end

if isdefined(:CUDArt)
function speedtest()
    x = h5read("devx.h5","/data")
    l1 = Op("dev1.h5")
    l2 = Op("dev2.h5")
    l1.xforw = KUnet.noop
    l2.xforw = KUnet.noop
    l1.w = CudaArray(l1.w)
    l1.b = CudaArray(l1.b)
    l2.w = CudaArray(l2.w)
    l2.b = CudaArray(l2.b)
    xrows,xcols = size(x)
    yrows,ycols = size(l2.w,1), xcols
    batch = 937
    y = similar(x, yrows, ycols)
    # xx = similar(l1.w,(xrows,batch))
    xx = CudaArray(x[:,1:batch])
    l1.y = similar(l1.w,(size(l1.w,1), batch))
    l2.y = similar(l2.w,(size(l2.w,1), batch))
    @time forward!(x, xx, l1.w, l1.b, l1.y, l2.w, l2.b, l2.y, y)
end
end

if isdefined(:CUDArt)
function speedtest2()
    batch = 937
    x = h5read("devx.h5","/data")
    l1 = Op("dev1.h5")
    l2 = Op("dev2.h5")
    l1.xforw = KUnet.noop
    l2.xforw = KUnet.noop
    l1.w = CudaArray(l1.w)
    l1.b = CudaArray(l1.b)
    l2.w = CudaArray(l2.w)
    l2.b = CudaArray(l2.b)
    # xrows,xcols = size(x)
    # yrows,ycols = size(l2.w,1), xcols
    # y = similar(x, yrows, ycols)
    # xx = similar(l1.w,(xrows,batch))
    # xx = CudaArray(x[:,1:batch])
    # l1.y = similar(l1.w,(size(l1.w,1), batch))
    # l2.y = similar(l2.w,(size(l2.w,1), batch))
    net = [l1,l2]
    @time KUnet.predict(net, x, batch)
    @time KUnet.predict(net, x, batch)
end
end

function speedtest3()
    # blas_set_num_threads(20)
    batch = 937
    x = h5read("devx.h5","/data")
    l1 = KUnet.Op("dev1.h5")
    l2 = KUnet.Op("dev2.h5")
    net = [l1,l2]
    @time KUnet.predict(net, x, batch=batch)
    @time KUnet.predict(net, x, batch=batch)
end

function speedtest4()
    # blas_set_num_threads(20)
    batch = 937
    x = h5read("devx.h5","/data")
    y = h5read("devy.h5","/data")
    l1 = KUnet.Op("dev1.h5")
    l2 = KUnet.Op("dev2.h5")
    net = [l1,l2]
    xx = x[:,1:batch]
    yy = y[:,1:batch]
    @time KUnet.backprop(net, xx, yy)
    @time KUnet.backprop(net, xx, yy)
    net
end

if isdefined(:CUDArt)
function speedtest5()
    # blas_set_num_threads(20)
    batch = 937
    x = h5read("devx.h5","/data")
    y = h5read("devy.h5","/data")
    l1 = KUnet.Op("dev1.h5")
    l2 = KUnet.Op("dev2.h5")
    l1.xforw = KUnet.noop
    l2.xforw = KUnet.noop
    l1.w = CudaArray(l1.w)
    l1.b = CudaArray(l1.b)
    l2.w = CudaArray(l2.w)
    l2.b = CudaArray(l2.b)
    net = [l1,l2]
    xx = CudaArray(x[:,1:batch])
    yy = CudaArray(y[:,1:batch])
    @time KUnet.backprop(net, xx, yy)
    @time KUnet.backprop(net, xx, yy)
    net
end
end

function speedtest6()
    # blas_set_num_threads(20)
    batch = 937
    x = h5read("devx.h5","/data")
    y = h5read("devy.h5","/data")
    l1 = KUnet.Op("rnd1.h5")
    l2 = KUnet.Op("rnd2.h5")
    net = [l1,l2]
    KUnet.setparam!(net, :l2reg, 0.5f0)
    @time KUnet.train(net, x, y; batch=937, iters=10)
    @time KUnet.train(net, x, y; batch=937, iters=10)
    net
end

if isdefined(:CUDArt)
function speedtest7()
    batch = 937
    x = h5read("devx.h5","/data")
    y = h5read("devy.h5","/data")
    l1 = KUnet.Op("rnd1.h5")
    l2 = KUnet.Op("rnd2.h5")
    l1.w = CudaArray(l1.w)
    l1.b = CudaArray(l1.b)
    l2.w = CudaArray(l2.w)
    l2.b = CudaArray(l2.b)
    net = [l1,l2]
    @time KUnet.train(net, x, y; batch=937, iters=1, l2reg=0.5f0)
    @time KUnet.train(net, x, y; batch=937, iters=1, l2reg=0.5f0)
    net
end
end


text:
using HydraulicFracturing
using Test

f(2,1)

@testset "HydraulicFracturing.jl" begin

    f(2,1)

end


text:
module TestWithSequential

using Transducers: DistributedEx, SequentialEx, ThreadedEx
using Folds.Testing: test_with_sequential, eval_call_data

executors = [
    SequentialEx(),
    ThreadedEx(),
    ThreadedEx(basesize = 1),
    DistributedEx(),
    DistributedEx(basesize = 1),
]

test_with_sequential(executors)

raw_extra_testdata = """
reduce(Mean(), 1:10)
reduce(Mean(), (x for x in 1:10 if isodd(x)))
reduce(KHist(25), 1.0:10.0)
"""

preamble = quote
    using OnlineStats
end

extra_testdata = eval_call_data(raw_extra_testdata, preamble)

test_with_sequential(extra_testdata, executors)

end


text:
using Test
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Density test" begin
    pde  = Helmholtz(;dim=3,k=1)
    Ω,M  = WaveProp.IO.gmsh_sphere(dim=2)
    Γ    = boundary(Ω)
    mesh = NystromMesh(view(M,Γ),order=1)
    σ     = γ₀(x->x[1],mesh)
    @test eltype(σ) == Float64
    σ     = γ₀(x->exp(im*x[3]),mesh)
    @test eltype(σ) == ComplexF64
    σ     = γ₁((x,n)->sum(x.*n),mesh)
    @test eltype(σ) == Float64
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
UTriangular(a::Matrix) = UpperTriangular(a)
if VERSION < v"0.5.0-dev"
    Base.cholfact(s::Symmetric) = cholfact(full(s))
end

## NumericExtensions is no longer supported, underoptimized implementation:
function logsumexp{T<:AbstractFloat}(x::AbstractVector{T})
    m = maximum(x)
    log(sum(exp.(x .- m))) + m
end
logsumexp{T<:AbstractFloat}(x::Matrix{T}, dim::Integer) = mapslices(logsumexp, x, dim)

## Also NumericExtensions' semantics of dot() is no longer supported.
function Base.dot{T<:AbstractFloat}(x::Matrix{T}, y::Matrix{T})
    size(x) == size(y) || error("Matrix sizes must match")
    dot(vec(x), vec(y))
end
function Base.dot{T<:AbstractFloat}(x::Matrix{T}, y::Matrix{T}, dim::Integer)
    size(x) == size(y) || error("Matrix sizes must match")
    if dim==1
        r = zeros(T, 1, size(x,2))
        for j in 1:size(x,2)
            for i in 1:size(x,1)
                r[j] += x[i,j]*y[i,j]
            end
        end
    else
        r = zeros(T, size(x,1), 1)
        for j in 1:size(x,2)
            for i in 1:size(x,1)
                r[i] += x[i,j]*y[i,j]
            end
        end
    end
    r
end

if VERSION < v"0.5.0-dev+2023"
    displaysize(io::IO) = Base.tty_size()
end

## this we need for xμTΛxμ!
#Base.A_mul_Bc!(A::StridedMatrix{Float64}, B::AbstractTriangular{Float32}) = A_mul_Bc!(A, convert(AbstractMatrix{Float64}, B))
#Base.A_mul_Bc!(A::Matrix{Float32}, B::AbstractTriangular{Float64}) = A_mul_Bc!(A, convert(AbstractMatrix{Float32}, B))
## this for diagstats
#Base.BLAS.gemm!(a::Char, b::Char, alpha::Float64, A::Matrix{Float32}, B::Matrix{Float64}, beta::Float64, C::Matrix{Float64}) = Base.BLAS.gemm!(a, b, alpha, float64(A), B, beta, C)


text:
using DataFrames
using Unitful
using StatsBase
using ATools

import Unitful:
    nm, μm, mm, cm, ns, μs, ms, ps, s

"""
    distance_lor_to_point(x1::Real,y1::Real,x2::Real,y2::Real,
                          x0::Real=0.0,y0::Real=0.0)

Compute the distance between lor defined by (x1,y1) and (x2, y2) and point (x0,y0)
"""
function distance_lor_to_point(x1::Real,y1::Real,x2::Real,y2::Real,
                               x0::Real=0.0,y0::Real=0.0)

    num = abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1))
    den = sqrt((x2 -x1)^2 + (y2 - y1)^2)
    num/den
end

"""
    function sinogramdf(lors)

Return a DataFrame from which sinograms can be computed
"""
function sinogramdf(lors)

    ldf = DataFrame(x1=[lors[i].x1 for i in 1:length(lors)],
                    y1=[lors[i].y1 for i in 1:length(lors)],
                    z1=[lors[i].z1 for i in 1:length(lors)],
                    x2=[lors[i].x2 for i in 1:length(lors)],
                    y2=[lors[i].y2 for i in 1:length(lors)],
                    z2=[lors[i].z2 for i in 1:length(lors)])

    # add z, θ and r columns
    ldf[!, "zl"] = ldf.z1 - ldf.z2;
    ldf[!, "rl"] = distance_lor_to_point.(ldf.x1, ldf.y1, ldf.x2, ldf.y2,
                                      zeros(nrow(ldf)), zeros(nrow(ldf)) )
    ldf[!, "tl"] = atan.(ldf.x1, ldf.y1)
    return ldf
end

function zxdf(sdf::DataFrame;  nbinz::Integer=3, zmin::Float64=-350.0, zmax::Float64=350.0)

    if nbinz == 1
        ZDF = [sdf]
    else
        # bins in z
        hz  = hist1d(sdf.zl, nbinz, zmin, zmax)
        hze = hz.edges[1]
        ZDF = [select_by_column_value_interval(sdf, "zl",
        hze[i], hze[i+1]) for i in 1:length(hze)-1]
    end
    return ZDF
end

"""
    ztsinogram(sdf; nbinz=3, zmin=-350.0, zmax=350.0,
                    nbint=3, tmin=-Float64(π), tmax=Float64(π))

Return a dictionary binning in z and θ for sinogram calculations
"""
function ztsinogram(sdf::DataFrame;
                    nbinz::Integer=3, zmin::Float64=-350.0, zmax::Float64=350.0,
                    nbint::Integer=3, tmin::Float64=-Float64(π), tmax::Float64=Float64(π))

    if nbinz == 1
        RL = [sdf]
    else
        # bins in z
        hz  = hist1d(sdf.zl, nbinz, zmin, zmax)
        hze = hz.edges[1]
        RL = [select_by_column_value_interval(sdf, "zl",
        hze[i], hze[i+1]) for i in 1:length(hze)-1]
    end

    # select in theta bins
    dRL = Dict()
    for (i, rl) in enumerate(RL)
        if nbint == 1
            RT =[rl]
        else
            ht  = hist1d(rl.zl, nbint, tmin, tmax)
            hte = ht.edges[1]
            RT = [select_by_column_value_interval(sdf,
                  "tl", hte[i], hte[i+1]) for i in 1:length(hte)-1]
        end
        dRL[i] = RT
    end
    return dRL
end

"""
    thetasgrm(drl; nbint, nbin=20, rmin=0.0, rmax=75.0)

Take the dictionary drl and return projected θ  sinograms
that is histograms in θ for each theta in the dictionary
"""
function thetasgrm(drl; nbint, nbin=20, rmin=0.0, rmax=75.0)
    nbinz = 1
    PR = []
    HR = []
    for i in 1:nbint
        sng = drl[1][i]
        h,p = hist1d(sng.rl, "rl", nbin, rmin, rmax)
        push!(PR, p)
        push!(HR, h)
    end
    return HR, PR
end

"""
    thetasgrm(df::DataFrame; ntproj::Integer, nbinr=20, rmin=0.0, rmax=75.0)

"""
function thetasgrm(df::DataFrame; ntproj::Integer, nbinr=20, rmin=0.0, rmax=75.0)

    ht  = hist1d(df.tl, ntproj, -Float64(π), Float64(π))
    hte = ht.edges[1]
    RT = [select_by_column_value_interval(df,
                                         "tl",
                                         hte[i], hte[i+1]) for i in 1:length(hte)-1]
    rH = [hist1d(dfx.rl, nbinr, rmin, rmax) for dfx in RT]
end

"""
    zsgrm(drl; nbinz, nbin=20, rmin=0.0, rmax=75.0)
Take the dictionary drl and return projected z sinograms
that is histograms in r for each z in the dictionary
"""
function zsgrm(drl; nbinz, nbin=20, rmin=0.0, rmax=75.0)
    nbint = 1
    PR = []
    HR = []
    for i in 1:nbinz
        sng = drl[i][1]
        h,p = hist1d(sng.rl, "rl", nbin, rmin, rmax)
        push!(PR, p)
        push!(HR, h)
    end
    return HR, PR
end

"""
    stp(hrl::Histogram; i0::Integer = 5, i1::Integer = 15)

Return Scatters, Trues and Prompts from an Histogram of nema4 radial distance
to source
    i0 :  defines the bin of background (left)
    i1 :  defines the bin of background (right)

Background is computed as linear interpolation between i0 and i1
"""
function stp(hrl::Histogram; i0::Integer = 5, i1::Integer = 15)
    w = Float64.(hrl.weights)
    e = hrl.edges[1]
    w0 = w[i0]                        # contents of bin
    w1 = w[i1]
    wm, im, em = find_max_xy(w, e)
    spb = w[i0:i1]
    fxy = gline2p(i0,w0,i1, w1)
    bkg = [fxy(i) for i in i0:i1]
    sgn = spb - bkg
    sgn = [s >0 ? ceil(s) : 0.0 for s in sgn]
    T = sum(sgn)
    P = sum(w)
    S = P - T

    @info "i0 = %d, im= %d, i1 = %d" i0 im i1
    @info "w0 = %5.1f, wm= %5.1f, w1 =%5.1f" w0 wm w1

    return S,T,P
end

"""
    pok1(λ::Float64)
    Poisson probability for k=0, 1 (k! = 1)

"""
pok1(λ::Float64) = exp(-λ) * (1.0 + λ)

"""
    pdt(rkcps::Float64, wmus::Float64)

Given an activity rkcps (in kcps) and a dead time wmus (in mus), the average number
of events in that window is kcps/(1/wmus) = kcps * wmus.
"""
function pdt(rkcps::Float64, wmus::Float64)
    λ = rkcps * wmus * 1e-3  # λ = rate/(1/w) = rate * w where rate in kcps, w in mus
    return pok1(λ)
end

"""
    cnec(n_kcps::Float64,
         eff_gp::Float64,
         eff_sp::Float64,
         SOT::Float64       =1.6,
         eff_cwns::Float64  = 1e-9,
        deadt_mus::Float64 = 0.6)

Computes NEC curve
"""
function cnec(n_kcps::Float64,
    eff_gp::Float64,
    eff_sp::Float64,
    SOT::Float64       =1.6,
    eff_cwns::Float64  = 1e-9,
    deadt_mus::Float64 = 0.6)

# n_kcps is the number of counts in kcps
# eff_gp is the efficiency for good prompts
# eff_sp is the efficiency for single prompts
# SOT = S/T, the ratio of plots to trues
# eff_cwns is the coincidence window in ns. Since rates come in kcps multiply by 10^3

eff_w =  eff_cwns *1e+3
nT = n_kcps * eff_gp               # number of good coincidences in scanner
nS = nT * SOT                      # number of scatters in scanner
n1 = n_kcps * eff_sp                # number of singles in scanner
nR = 2 * n1^2 * eff_w              # number of randoms
ntot = nT + nS + nR

prob = pdt(ntot, deadt_mus)
nT   *= prob
nS   *= prob
#nR   *= prob
ntot *= prob

nec = nT / (1.0 + nS/nT + nR/nT)   # nec

return nT, nS, nR, ntot, nec, prob
end

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
ENV["PYTHON"] = ""
using PyCall
using Conda

pipcmd = joinpath(Conda.PYTHONDIR,"pip")

#Conda.add("astroquery.jplhorizons")
Conda.add("astroquery", channel="conda-forge")
#run(`$pipcmd install --pre astroquery`)
#Conda.add("astroquery.jplhorizons", channel="conda-forge")

run(`$pipcmd install barycorrpy`)

Conda.add("astropy")
#Conda.add("astropy.coordinates")
#Conda.add("astropy.time")

Conda.add("scipy")
#Conda.add("scipy.spatial.transform")

#run(`$pipcmd install pyneid`)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the Jacchia-Bowman 2008 model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] http://sol.spacenvironment.net/~JB2008/
#   [2] http://sol.spacenvironment.net/~JB2006/
#   [3] https://github.com/JuliaSpace/JB2008_Test
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/atmospheric_models/jb2008
# ===========================================

# Functions: jb2008
# -----------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
#   Compare the results with the tests in files:
#
#       JB_AUTO_OUTPUT_01.DAT
#       JB_AUTO_OUTPUT_02.DAT
#       JB_AUTO_OUTPUT_03.DAT
#       JB_AUTO_OUTPUT_04.DAT
#
#   which were created using JuliaSpace/JB2008_Test@1552551 .
#
################################################################################

@testset "Function jb2008" begin
    # Scenario 01
    # ===========

    # Initialize the space indices with the local files.
    init_space_indices(enabled_files  = [:dtcfile, :fluxtable, :solfsmy],
                       fluxtable_path = "./fluxtable_demo.txt",
                       dtcfile_path   = "./DTCFILE_DEMO.TXT",
                       solfsmy_path   = "./SOLFSMY_DEMO.TXT")
    # TODO: The last line is a workaround to do not download the WDC files,
    # which we do not need for this test. This must be modified once a way to
    # avoid those unnecessary downloads is implemented.

    # Files with the test results.
    #
    # TODO: Two tests are disable due to the lack of data in
    # `fluxtable_demo.txt`.
    test_list = ["JB2008_AUTO_OUTPUT_01.DAT",
                 # "JB2008_AUTO_OUTPUT_02.DAT",
                 "JB2008_AUTO_OUTPUT_03.DAT",
                 #"JB2008_AUTO_OUTPUT_04.DAT"
                ]

    # Execute the tests.
    for filename in test_list
        open(filename) do file
            line_num = 0
            year     = 0
            doy      = 0
            hour     = 0
            min      = 0
            sec      = 0.0

            for line in eachline(file)
                line_num += 1

                # Ignore 5 lines to skip the header.
                (line_num <= 5) && continue

                # Ignore others non important lines.
                (line_num in [7, 8, 9]) && continue

                if line_num == 6
                    # Read the next line to obtain the input data related to the
                    # time.
                    tokens = split(line)

                    year = parse(Int,     tokens[1])
                    doy  = parse(Int,     tokens[2])
                    hour = parse(Int,     tokens[3])
                    min  = parse(Int,     tokens[4])
                    sec  = parse(Float64, tokens[5])
                else
                    tokens = split(line)

                    # Read the information about the location.
                    h    = parse(Float64, tokens[1])*1000
                    glon = parse(Float64, tokens[2])*pi/180
                    glat = parse(Float64, tokens[3])*pi/180

                    # Read the model output.
                    T_exo = parse(Float64, tokens[4])
                    Tz    = parse(Float64, tokens[5])
                    ρ     = parse(Float64, tokens[6])

                    # Run the model.
                    JD  = date_to_jd(year, 1, 1, hour, min, sec) - 1 + doy
                    out = jb2008(JD, glat, glon, h)

                    # Compare the results.
                    @test out.T_exo ≈ T_exo atol = 0.6 rtol = 0.0
                    @test out.Tz    ≈ Tz    atol = 0.6 rtol = 0.0
                    @test out.rho   ≈ ρ     atol = 0.0 rtol = 5e-3
                end
            end
        end
    end
end


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
# File   : sphere_approx.jl
# License: MIT
# Author : Andrei Leonard Nicusan <a.l.nicusan@bham.ac.uk>
# Date   : 21.05.2021

using GLMakie
using Polyhedra
using ColorSchemes
using LinearAlgebra

import QHull
lib = QHull.Library()

# Return points on a sphere, generated with the method of Markus Deserno from
# https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
function sphere(;centre = [0., 0., 0.], radius = 1., n = 16)
    nc = 0
    a = 4π * radius^2 / n
    d = √a

    mθ = floor(π / d)

    dθ = π / mθ
    dϕ = a / dθ

    points = zeros(n, 3)

    for m in 0:mθ - 1
        θ = π * (m + 0.5) / mθ
        mϕ = floor(2π * sin(θ) / dϕ)
        for p in 0:mϕ - 1
            ϕ = 2π * p / mϕ
            nc += 1

            points[nc, 1] = centre[1] + radius * sin(θ) * cos(ϕ)
            points[nc, 2] = centre[2] + radius * sin(θ) * sin(ϕ)
            points[nc, 3] = centre[3] + radius * cos(θ)
        end
    end

    return points[1:nc, :]
end

# Very dense circle points
sphere_points = sphere(n = 512)
sphere_poly = polyhedron(vrep(sphere_points), lib)

# Plot polygonal approximations and area ratio
cm = ColorSchemes.magma
# GLMakie.mesh(Polyhedra.Mesh{3}(sphere_poly), color = get(cm, 0.))

nverts = collect(10:25:500)   # [16, 32, 64, 128, 256]
ratios = zeros(length(nverts))
colors = [get(cm, i) for i in LinRange(0.1, 0.9, length(nverts))]

for (i, n) in enumerate(nverts)
    local sphere_points = sphere(centre = [0., 0. - i / 10, 0.], n = n)

    local sphere_poly = polyhedron(vrep(sphere_points), lib)
    # GLMakie.mesh!(Polyhedra.Mesh{3}(sphere_poly), color = colors[i])

    ratios[i] = Polyhedra.volume(sphere_poly) / (4 / 3 * pi)
end

current_figure()

# Ploting
GLMakie.lines(nverts, ratios)

max_stars_repo_name: anicusan/KonigCell

max_issues_repo_name: anicusan/KonigCell

max_forks_repo_name: anicusan/KonigCell


text:
export LongstepWrapper

include("saveplanes.jl")

mutable struct LongstepWrapper{T<:FOSAlgorithm} <: FOSAlgorithm
    longinterval::Int64
    nsave::Int64
    alg::T
    options
end

mutable struct LongstepWrapperData{T<:FOSSolverData} <: FOSSolverData
    longinterval::Int64
    nsave::Int64
    saved::SavedPlanes{Float64}
    saveineq::Bool
    savepos::Int64
    eqi::Int64
    uneqi::Int64
    tmp::Array{Float64,1}
    algdata::T
end

function LongstepWrapper(alg::T; longinterval=100, nsave=10, kwargs...) where T
    LongstepWrapper{T}(longinterval, nsave, alg, [kwargs...,alg.options...])
end

function init_algorithm!(long::LongstepWrapper, model::AbstractFOSModel)
    alg, longinterval, nsave = long.alg, long.longinterval, long.nsave
    !support_longstep(alg) && @error "Algorithm alg does not support longstep"
    data, status_generator =  init_algorithm!(alg, model)
    neq, nineq = projections_per_step(alg)
    #TODO x
    x = getinitialvalue(model)
    saved = SavedPlanes(x, nsave, neq, nineq)
    # TODO Fix case with saveineq=false
    data = LongstepWrapperData(longinterval, nsave, saved, true, 0, 1, 1, similar(x), data)
    return data, status_generator
end

getmn(data::LongstepWrapperData) = getmn(data.algdata)

function Base.step(wrap::LongstepWrapper, longstep::LongstepWrapperData, x, i, status::AbstractStatus)
    nsave, longinterval = longstep.nsave, longstep.longinterval
    #println("At iteration $i")
    savepos = (i-1)%longinterval - longinterval + nsave + 2 # Index of collection to save
    #println("savepos: $savepos")
    if 0 < savepos  #Should we start saving collections?
        longstep.savepos = savepos
        longstep.eqi, longstep.uneqi = 0, 0
    end
    step(wrap.alg, longstep.algdata, x, i, status, longstep)
    #If time to do longstep
    if longstep.savepos == nsave+1
        #println("Doing projection!")
        fail = projectonnormals!(longstep.saved, x, longstep.tmp)
        longstep.savepos = -1
        x .= longstep.tmp
    end
end

function getsol(alg::LongstepWrapper, data::LongstepWrapperData, x)
    getsol(alg.alg, data.algdata, x)
end

getcgiter(data::LongstepWrapperData) = getcgiter(data.algdata)

addprojeq(::Nothing, ::Any, ::Any) = nothing
addprojineq(::Nothing, ::Any, ::Any) = nothing

function addprojeq(long::T, y, x) where T<:LongstepWrapperData
    if long.savepos > 0
        #println("adding EQ, nsave: $(long.nsave) i:$(long.i) eqi:$(long.eqi) uneqi:$(long.uneqi)")
        s = long.saved
        i = (long.savepos-1)*(s.neq + s.nineq*long.saveineq) + long.eqi + 1
        s.A[i,:] .= x .- y
        b = 0.0
        for j = 1:length(x)
            b += (x[j]-y[j])*y[j]
        end
        s.b[i] = b
        long.eqi += 1
    end
    return
end

"""Saves inequalities as equality for now"""
function addprojineq(long::T, y, x) where T<:LongstepWrapperData
    if long.savepos > 0 && (long.saveineq || long.savepos == long.nsave+1)
        #println("adding inEQ, nsave: $(long.nsave) i:$(long.i) eqi:$(long.eqi) uneqi:$(long.uneqi)")
        s = long.saved
        i = (long.savepos-1)*(s.neq + s.nineq*long.saveineq) + long.eqi + long.uneqi + 1
        s.A[i,:] .= x .- y
        b = 0.0
        for j = 1:length(x)
            b += (x[j]-y[j])*y[j]
        end
        s.b[i] = b
        long.uneqi += 1
    end
end


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

mutable struct AuthorizationPropertiesFormat <: SwaggerModel
    authorizationKey::Any # spec type: Union{ Nothing, String } # spec name: authorizationKey
    authorizationUseStatus::Any # spec type: Union{ Nothing, String } # spec name: authorizationUseStatus
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function AuthorizationPropertiesFormat(;authorizationKey=nothing, authorizationUseStatus=nothing, provisioningState=nothing)
        o = new()
        validate_property(AuthorizationPropertiesFormat, Symbol("authorizationKey"), authorizationKey)
        setfield!(o, Symbol("authorizationKey"), authorizationKey)
        validate_property(AuthorizationPropertiesFormat, Symbol("authorizationUseStatus"), authorizationUseStatus)
        setfield!(o, Symbol("authorizationUseStatus"), authorizationUseStatus)
        validate_property(AuthorizationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type AuthorizationPropertiesFormat

const _property_map_AuthorizationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("authorizationKey")=>Symbol("authorizationKey"), Symbol("authorizationUseStatus")=>Symbol("authorizationUseStatus"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_AuthorizationPropertiesFormat = Dict{Symbol,String}(Symbol("authorizationKey")=>"String", Symbol("authorizationUseStatus")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ AuthorizationPropertiesFormat }) = collect(keys(_property_map_AuthorizationPropertiesFormat))
Swagger.property_type(::Type{ AuthorizationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AuthorizationPropertiesFormat[name]))}
Swagger.field_name(::Type{ AuthorizationPropertiesFormat }, property_name::Symbol) =  _property_map_AuthorizationPropertiesFormat[property_name]

const _allowed_AuthorizationPropertiesFormat_authorizationUseStatus = ["Available", "InUse"]

function check_required(o::AuthorizationPropertiesFormat)
    true
end

function validate_property(::Type{ AuthorizationPropertiesFormat }, name::Symbol, val)
    if name === Symbol("authorizationUseStatus")
        Swagger.validate_param(name, "AuthorizationPropertiesFormat", :enum, val, _allowed_AuthorizationPropertiesFormat_authorizationUseStatus)
    end
end


text:
using ImageMethodReverb, Random
Fs = 8e3                        # sampling frequency
L = Lx, Ly, Lz =  4.,5.,3.      # room dimensions in meters
T60 = 0.7                       # reverberation time
# or alternatively
#β = (0.9,0.9,0.7,0.7,0.8,0.8)
#(  βx1     ,    βx2   ,    βy1   ,    βy2    ,  βz1 ,   βz2  )
#(front wall, rear wall, left wall, right wall, floor, ceiling)

xs = (0.5, 0.5, 0.5)          #src pos (in meters)
xr = (Lx-0.1, Ly-0.3, Lz-0.2) #mic pos
Nt = round(Int,Fs/2)          #time samples (1/5 sec)

#h, = rim(xs,xr,L,β,Nt,Fs)
h, = rim(xs,xr,L,T60,Nt,Fs)
t = range(0; length = Nt, step = 1/Fs )

# more mics
h, = rim(xs,[(1,1,1),(1,2,1)],L,T60,Nt,Fs)

#with source signal
s = randn(Nt)   #src signal
p, = rim(s,xs,[(1,1,1),(1,1.5,1)],L,T60,Nt,Fs)

using Plots
p = plot(t,h, xlabel="Time", ylabel="Sound Pressure")


text:
"""
    DevicePtr{T,A}
A memory address that refers to data of type `T` that is accessible from the GPU. It is the
on-device counterpart of `MtlDAdrv.MtlPtr`, additionally keeping track of the address space
`A` where the data resides (shared, global, constant, etc). This information is used to
provide optimized implementations of operations such as `unsafe_load` and `unsafe_store!.`
"""
DevicePtr

# constructors
DevicePtr{T,A}(x::Union{Int,UInt,DevicePtr}) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, x)
#DevicePtr{T}(ptr::MtlPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)
#DevicePtr(ptr::MtlPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

## conversions
Base.convert(::Type{DevicePtr{T,A}}, x::Union{Int,UInt}) where {T,A<:AddressSpace} = DevicePtr{T,A}(x)

# between host and device pointers
#Base.convert(::Type{MtlPtr{T}},  p::DevicePtr)  where {T}                 = Base.bitcast(MtlPtr{T}, p)
Base.convert(::Type{MtlBuffer{T}},   p::DevicePtr)  where {T}               = MtlBuffer{T}(Base.bitcast(MTL.MTLBuffer, p))
#Base.convert(::Type{DevicePtr{T,A}}, p::MtlPtr) where {T,A<:AddressSpace}   = Base.bitcast(DevicePtr{T,A}, p)
#Base.convert(::Type{DevicePtr{T}},   p::MtlPtr)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p)
#Base.convert(::Type{DevicePtr{T,A}}, p::MtlBuffer) where {T,A<:AddressSpace}   = Base.bitcast(DevicePtr{T,A}, p.handle)
#Base.convert(::Type{DevicePtr{T}},   p::MtlBuffer)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p.handle)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{Ptr{T}}, x::DevicePtr{T}) where {T} = reinterpret(Ptr{T}, x)

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = Base.bitcast(DevicePtr{T,A}, p)
## identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
## convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

## memory operations

@generated function pointerref(p::DevicePtr{T,A}, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return T.instance
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ld = load!(builder, ptr)

        if A != AS.Generic
            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(ld, align)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,A}, Int}, :((p, Int(i-one(i)))))
end

@generated function pointerset(p::DevicePtr{T,A}, x::T, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, eltyp, T_int]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[3]])
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)

        if A != AS.Generic
            metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(st, align)

        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{DevicePtr{T,A}, T, Int},
                  :((p, convert(T,x), Int(i-one(i)))))
end

# operand types supported by llvm.nvvm.ldg.global
# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet
const LDGTypes = Union{UInt8, UInt16, UInt32, UInt64,
                       Int8, Int16, Int32, Int64,
                       Float32, Float64}

# interface

export unsafe_cached_load

Base.unsafe_load(p::DevicePtr{T}, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerref(p, Int(i), align)

Base.unsafe_store!(p::DevicePtr{T}, x, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerset(p, convert(T, x), Int(i), align)

# NOTE: fall back to normal pointerref for unsupported types. we could be smarter here,
#       e.g. destruct/load/reconstruct, but that's too complicated for what it's worth.
unsafe_cached_load(p::DevicePtr, i::Integer=1, align::Val=Val(1)) =
    pointerref(p, Int(i), align)


text:
using QuasiMonteCarlo, Distributions, Test
using Test

#1D
lb = 0.0
ub = 5.0
n = 5
d = 1
QuasiMonteCarlo.sample(n,lb,ub,GridSample(0.1))
QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
QuasiMonteCarlo.sample(n,lb,ub,LatticeRuleSample())
QuasiMonteCarlo.sample(20,lb,ub,LowDiscrepancySample(10))
QuasiMonteCarlo.sample(5,d,Cauchy())
QuasiMonteCarlo.sample(5,d,Normal(0,4))

#ND
lb = [0.1,-0.5]
ub = [1.0,20.0]
n = 5
d = 2

@testset "GridSample" begin
    #GridSample{T}
    s = QuasiMonteCarlo.sample(n,lb,ub,GridSample([0.1,0.5]))
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "UniformSample" begin
    #UniformSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "SobolSample" begin
    #SobolSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LHS" begin
    #LHS
    s = QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LatticeRuleSample" begin
    #LatticeRuleSample()
    s = QuasiMonteCarlo.sample(n,lb,ub,LatticeRuleSample())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "LDS" begin
    #LDS
    s = QuasiMonteCarlo.sample(n,lb,ub,LowDiscrepancySample([10,3]))
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "Distribution 1" begin
    #Distribution 1
    s = QuasiMonteCarlo.sample(n,d,Cauchy())
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "Distribution 2" begin
    #Distribution 2
    s = QuasiMonteCarlo.sample(n,d,Normal(3,5))
    @test isa(s,Matrix{typeof(s[1][1])}) == true
    @test size(s) == (d, n)
end

@testset "generate_design_matrices" begin
    num_mat = 3
    Ms = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,UniformSample(), num_mat)
    @test length(Ms) == num_mat
    A = Ms[1]
    B = Ms[2]
    @test isa(A, Matrix{typeof(A[1][1])}) == true
    @test size(A) == (d, n)
    @test isa(B, Matrix{typeof(B[1][1])}) == true
    @test size(B) == (d, n)
end


text:
#=
$ docker-compose run --rm julia
julia> pwd()
julia> /work
julia> include("experiments/test/runtests.jl")
=#

using Test
using Glob
using Base.Threads

ignore_files = ["wav_example.md", "clang.md", "linear_regression.md"]

@testset "MyWorkflow.jl" begin
    files = glob("*.md", joinpath(@__DIR__, "..", "notebook")) |> sort
    @threads for f in files
        basename(f) in ignore_files && continue
        @info "Running $f"
        proc = run(`jupytext --to ipynb --execute $f`)
        @test proc.exitcode == 0
    end
end


text:
function triangular(n)
  sum = 0
  for i in 1:n
    sum = sum + i
  end
  sum
end

function countDivisors(n)
  factors = 0
  for i in 1:floor(√n)
    if mod(n,i) == 0
      factors = factors + 1
    end
  end
  factors = factors * 2
  if √n == floor(√n)
    factors = factors - 1
  end
  factors
end

for i in 1:100000000000000000000000000
  if countDivisors(triangular(i)) >= 500
    println(triangular(i))
    quit()
  end
end

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
using Mocha

include("VAE.jl")

backend = DefaultBackend()
init(backend)

# Number of latent variables, size of encoding, decoding layers
# The Kingma & Welling 2013 paper uses 50,500,500 when optimizing likelihood.
# Here we use low dimensional latent space for exploration
N_Z = 50
net = make_vae(backend, N_Z, 400, 400)

############# Train the model ############

base_dir = "snapshots_mnist_VAE_exp_50_400"

method = Adam()
params = make_solver_parameters(method, max_iter=20000, regu_coef=0.0,
                                lr_policy=LRPolicy.Inv(0.01, 0.0001, 0.5),
                                #lr_policy=LRPolicy.Fixed(0.002),
                                load_from=base_dir)
solver = Solver(method, params)

setup_coffee_lounge(solver, save_into="$base_dir/statistics.jld", every_n_iter=5000)
add_coffee_break(solver, TrainingSummary(:iter, :obj_val, :learning_rate, "kl-loss", "bce-loss"), every_n_iter=100)
add_coffee_break(solver, Snapshot(base_dir), every_n_iter=5000)

solve(solver, net)

############# Play with it ##############

# We use three sliders to choose a value of Z, the latent variable.
# Then we use the network weights to decode this to an image.

using Interact, Winston
# (Need a recent Interact.jl on Ipython 4/Jupyter - https://github.com/JuliaLang/Interact.jl/issues/73)
# Pkg.checkout("Interact.jl"), and make sure ipywidgets Python package is installed.

latent_to_output = make_latent_to_output(net, 12, 13)

xx = -1:0.01:1
@Interact.manipulate for x in xx, y in xx, z in xx
  imagesc(latent_to_output([x,y,z]))
end

destroy(net)
shutdown(backend)


text:
# Test for gausschebyshev().

n = 10

# Test x.^2:
x, w = gausschebyshev(n,1)
@test dot(x.^2,w) ≈ pi/2
x, w = gausschebyshev(n,2)
@test dot(x.^2,w) ≈ pi/8
x, w = gausschebyshev(n,3)
@test dot(x.^2,w) ≈ pi/2
x, w = gausschebyshev(n,4)
@test dot(x.^2,w) ≈ pi/2

# Test x^3:
x, w = gausschebyshev(n,1)
@test abs(dot(x.^3,w))<1e-15
x, w = gausschebyshev(n,2)
@test abs(dot(x.^3,w))<1e-15
x, w = gausschebyshev(n,3)
@test dot(x.^3,w) ≈ 3*pi/8
x, w = gausschebyshev(n,4)
@test dot(x.^3,w) ≈ -3*pi/8


text:
## Comparison of exogenous grid search and endogenous grid search
include("Step1Grid.jl")
include("Step1EndGrid.jl")

using BenchmarkTools

# Set parameters
β = 0.9932
σ = 1.5

# Exogenous shocks
eH = 1.0
eL = 0.1
E = [eL; eH]

# Transition matrix: Pr(j|i) = Π(i,j) (columns sum to 1)
πHH = 0.925 # e_H | e_H
πHL = 0.5   # e_H | e_L
Π = [1 - πHL  πHL
    1 - πHH πHH]

# Assets
a_lb = -2.0
a_max = 5
a_size = 350
A = range(a_lb, a_max, length = a_size)
A = [A; ];

# Guess initial price
q = 1.1

# Check performance simple grid search method
@btime step1(A, E, Π, q, β, σ)

# Check performance endogenous grid search method
@btime EndGridSearch(A, E, Π, β, σ, q)


text:
"""
	gmtread(fname::String; kwargs...)

Read GMT object from file. The object is one of "grid" or "grd", "image" or "img",
"data" or "table", "cmap" or "cpt" and "ps" (for postscript), and OGR formats (shp, kml, json).
Use a type specificatin to force a certain reading path (e.g. `grd=true` to read grids) or take
the chance of letting the data type be guessed via the file extension. Known extensions are:

- Grids:      .grd, .nc
- Images:     .jpg, .png, .tif, .tiff, .bmp, .webp
- Datasets:   .dat, .txt, .csv
- Datasets:   .shp, .kml, .json, .geojson, .gmt, .gpkg
- CPT:        .cpt
- PostScript: .ps, .eps

Parameters
----------

Specify data type (with *type*=true, e.g. `img=true`).  Choose among:
- **grd** | **grid** :: [Type => Any]

    Load a grid.
- **img** | **image** :: [Type => Any]

    Load an image.
- **cpt** | **cmap** :: [Type => Any]

    Load a GMT color palette.
- **data** | **dataset** | **table** :: [Type => Any]

    Load a dataset (a table of numbers).
- **ogr** :: [Type => Any]

    Load a dataset via GDAL OGR (a table of numbers). Many things can happen here.
- **ps** :: [Type => Any]

    Load a PostScript file
- **gdal** :: [Type => Any]

    Force reading the file via GDAL. Should only be used to read grids.
- **varname** :: [Type => Str]

    When netCDF files have more than one 2D (or higher) variables use *varname* to pick the wished variable.
    e.g. varname=:slp to read the variable named 'slp'. This option defaults data type to 'grid'
- **layer** | **band** :: [Type => Str, Number, Array]

When files are multiband or nc files with 3D or 4D arrays, we access them via these keywords.
`layer=4` reads the fourth layer (or band) of the file. But the file can be a grid or an image. If it is a
grid, layer can be a scalar (to read 3D arrays) or an array of two elements (to read a 4D array).

If file is an image `layer` can be a 1 or a 1x3 array (to read a RGB image). Not that in this later case
bands do not need to be contiguous. A `band=[1,5,2]`` composes an RGB out of those bands. See more at
$(GMTdoc)/GMT_Docs.html#modifiers-for-coards-compliant-netcdf-files) but note that we use **1 based** indexing here.

- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_f)

Example: to read a nc called 'lixo.grd'

    G = gmtread("lixo.grd");

to read a jpg image with the bands reversed (this example is currently broken in GMT5. Needs GMT6dev)

    I = gmtread("image.jpg", band=[2,1,0]);
"""
function gmtread(fname::String; kwargs...)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, opt_R = parse_R(d, "")
	cmd = parse_common_opts(d, cmd, [:V_params :f :i :h])[1]
	cmd, opt_bi = parse_bi(d, cmd)

	# Process these first so they may take precedence over defaults set below
	opt_T = add_opt(d, "", "Tg", [:grd :grid])
	if (opt_T != "")		# Force read via GDAL
		if ((val = find_in_dict(d, [:gdal])[1]) !== nothing)  fname = fname * "=gd"  end
	else
		opt_T = add_opt(d, "", "Ti", [:img :image])
	end
	if (opt_T == "")  opt_T = add_opt(d, "", "Td", [:data :dataset :table])  end
	if (opt_T == "")  opt_T = add_opt(d, "", "Tc", [:cpt :cmap])  end
	if (opt_T == "")  opt_T = add_opt(d, "", "Tp", [:ps])   end
	if (opt_T == "")  opt_T = add_opt(d, "", "To", [:ogr])  end

	ogr_layer = Int32(0)			# Used only with ogrread. Means by default read only the first layer
	if ((varname = find_in_dict(d, [:varname])[1]) !== nothing) # See if we have a nc varname / layer request
		if (opt_T == "")			# Force read via GDAL
			if ((val = find_in_dict(d, [:gdal])[1]) !== nothing)  fname = fname * "=gd"  end
			opt_T = " -Tg"
		end
		fname = fname * "?" * arg2str(varname)
		if ((val = find_in_dict(d, [:layer :band])[1]) !== nothing)
			if (isa(val, Number))     fname *= @sprintf("[%d]", val-1)
			elseif (isa(val, Array))  fname *= @sprintf("[%d,%d]", val[1]-1, val[2]-1)	# A 4D array
			end
		end
	else									# See if we have a bands request
		if ((val = find_in_dict(d, [:layer :layers :band :bands])[1]) !== nothing)
			if ((lix = guess_T_from_ext(fname)) == " -To")		# See if it's a OGR layer request
				ogr_layer = Int32(val - 1)	# -1 because it's going to be sent to C (zero based)
			else
				fname = fname * "+b"
				if (isa(val, String) || isa(val, Symbol) || isa(val, Number))
					fname = string(fname, parse(Int, string(val))-1)
				elseif (isa(val, Array) || isa(val, Tuple))
					# Replacement for the annoying fact that one cannot do @sprintf(repeat("%d,", n), val...)
					fname  *= @sprintf("%d", val[1]-1)
					[fname *= @sprintf(",%d", val[k]-1) for k = 2:length(val)]
				end
				(opt_T == "") && (opt_T = " -Ti")
			end
		end
	end

	if (opt_T == "")
		((opt_T = guess_T_from_ext(fname)) == "") && error("Must select one input data type (grid, image, dataset, cmap or ps)")
		(opt_T == " -Tg" && haskey(d, :ignore_grd)) && return nothing	# contourf uses this
		if (opt_T == " -To" && fname[1] == '@')	# Because GMT ogrread does not know the '@' mechanism.
			fname = joinpath(readlines(`gmt --show-userdir`)[1], "cache", fname[2:end])
			!isfile(fname) && error("File $fname does not exist in cache.")
		end
	else
		opt_T = opt_T[1:4]      				# Remove whatever was given as argument to type kwarg
	end

	if (opt_T == " -Ti" || opt_T == " -Tg")		# See if we have a mem layout request
		if ((val = find_in_dict(d, [:layout :mem_layout])[1]) !== nothing)
			(opt_T == " -Ti" && startswith(string(val), "TRB")) && return gdaltranslate(fname)
			# MUST TAKE SOME ACTION HERE. FOR IMAGES I THINK ONLY THE "I" FOR IMAGES.JL IS REALLY POSSIBLE
			cmd = (opt_T == " -Ti") ? cmd * " -%" * arg2str(val) : cmd * " -&" * arg2str(val)
		end
	end

	if (opt_T != " -To")			# All others but OGR
		if (opt_T == " -Td" && opt_bi != "")  cmd *= opt_bi  end		# Read from binary file
		cmd *= opt_T
		return (dbg_print_cmd(d, cmd) !== nothing) ? "gmtread " * cmd : gmt("read " * fname * cmd)
	end

	(dbg_print_cmd(d, cmd) !== nothing) && return "ogrread " * fname * " " * cmd
	# Because of the certificates shits on Windows. But for some reason the set in gmtlib_check_url_name() is not visible
	(Sys.iswindows()) && run(`cmd /c set GDAL_HTTP_UNSAFESSL=YES`)
	API2 = GMT_Create_Session("GMT", 2, GMT_SESSION_NOEXIT + GMT_SESSION_EXTERNAL + GMT_SESSION_NOGDALCLOSE + GMT_SESSION_COLMAJOR);
	if (GMTver >= v"6.1")
		x = opt_R2num(opt_R)		# See if we have a limits request
		if (GMTver > v"6.1.1")
			lims = (x === nothing) ? (0.0, 0.0, 0.0, 0.0, 0.0, 0.0) : tuple(vcat(x,[0.0, 0.0])...)
			ctrl = OGRREAD_CTRL(Int32(0), ogr_layer, pointer(fname), lims)
			O = ogr2GMTdataset(gmt_ogrread(API2, pointer([ctrl])))
		else
			O = ogr2GMTdataset(gmt_ogrread(API2, fname, (x === nothing) ? C_NULL : x))
		end
	end
	(GMTver == v"6.0") && (O = ogr2GMTdataset(gmt_ogrread(API2, fname)))
	ressurectGDAL()				# Because GMT called GDALDestroyDriverManager()
	GMT_Destroy_Session(API2)
	return O
end

# ---------------------------------------------------------------------------------
function guess_T_from_ext(fname::String)::String
	# Guess the -T option from a couple of known extensions
	ext = splitext(fname)[2]
	(length(ext) > 8 || occursin("?", ext)) && return (occursin("?", ext)) ? " -Tg" : "" # A SUBDATASET encoded fname?
	ext = lowercase(ext[2:end])
	((ext == "jp2" || ext == "tif" || ext == "tiff") && (!isfile(fname) && !startswith(fname, "/vsi") &&
		!occursin("https:", fname) && !occursin("http:", fname) && !occursin("ftps:", fname) && !occursin("ftp:", fname))) &&
		error("File $fname does not exist.")
	if     (findfirst(isequal(ext), ["grd", "nc", "nc=gd"])  !== nothing)  out = " -Tg";
	elseif (findfirst(isequal(ext), ["dat", "txt", "csv"])   !== nothing)  out = " -Td";
	elseif (findfirst(isequal(ext), ["jpg", "png", "bmp", "webp"]) 	!== nothing)  out = " -Ti";
	elseif (findfirst(isequal(ext), ["shp", "kml", "json", "geojson", "gmt", "gpkg"]) !== nothing)  out = " -To";
	elseif (ext == "jp2") ressurectGDAL(); out = (findfirst("Type=UInt", gdalinfo(fname)) !== nothing) ? " -Ti" : " -Tg"
	elseif (ext == "cpt")  out = " -Tc";
	elseif (ext == "ps" || ext == "eps")  out = " -Tp";
	elseif (startswith(ext, "tif"))
		ressurectGDAL();
		gdinfo = gdalinfo(fname)
		(gdinfo === nothing) && error("gdalinfo failed - unable to open $fname")
		out = (findfirst("Type=UInt", gdinfo) !== nothing || findfirst("Type=Byte", gdinfo) !== nothing) ? " -Ti" : " -Tg"
	else
		out = ""
	end
end

"""
	gmtwrite(fname::AbstractString, data; kwargs...)

Write a GMT object to file. The object is one of "grid" or "grd", "image" or "img",
"dataset" or "table", "cmap" or "cpt" and "ps" (for postscript).

When saving grids we have a panoply of formats at our disposal.

Parameters
----------

- **id** ::  [Type => Str]

    Use an ``id`` code when not not saving a grid into a standard COARDS-compliant netCDF grid. This ``id``
    is made up of two characters like ``ef`` to save in ESRI Arc/Info ASCII Grid Interchange format (ASCII float).
    See the full list of ids at $(GMTdoc)grdconvert.html#format-identifier.

    ($(GMTdoc)grdconvert.html#g)
- **scale** | **offset** :: [Type => Number]

    You may optionally ask to scale the data and then offset them with the specified amounts.
    These modifiers are particularly practical when storing the data as integers, by first removing an offset
    and then scaling down the values.
- **nan** | **novalue** | **invalid** | **missing** :: [Type => Number]

    Lets you supply a value that represents an invalid grid entry, i.e., ‘Not-a-Number’.
- **gdal** :: [Type => Bool]

    Force the use of the GDAL library to write the grid (to be used only with grids).
    ($(GMTdoc)GMT_Docs.html#grid-file-format-specifications)
- **driver** :: [Type => Str]

    When saving in other than the netCDF format we must tell the GDAL library what is wished format.
    That is done by specifying the driver name used by GDAL itself (e.g., netCDF, GTiFF, etc...).
- **datatype** :: [Type => Str] 		``Arg = u8|u16|i16|u32|i32|float32``

    When saving with GDAL we can specify the data type from u8|u16|i16|u32|i32|float32 where ‘i’ and ‘u’ denote
    signed and unsigned integers respectively.
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bo)
- $(GMT.opt_f)

Example: write the GMTgrid 'G' object into a nc file called 'lixo.grd'

	gmtwrite("lixo.grd", G);
"""
gmtwrite(data, fname::AbstractString=""; kwargs...) = gmtwrite(fname, data; kwargs...)
function gmtwrite(fname::AbstractString, data; kwargs...)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_R(d, "")
	cmd = parse_V_params(d, cmd)
	if (fname == "")
		opt_G = parse_G(d, "")[1]
		(opt_G != "") && (fname = opt_G[4:end])	# opt_G comes with the " -G" prefix
	end
	(fname == "") && error("Output file name cannot be empty.")

	if (isa(data, GMTgrid))
		opt_T = " -Tg"
		fname = fname * parse_grd_format(d)		# If we have format requests
		cmd, = parse_f(d, cmd)
		CTRL.proj_linear[1] = true				# To force pad=0 and julia memory (no dup)
	elseif (isa(data, GMTimage))
		opt_T = " -Ti"
		fname *= parse_grd_format(d)			# If we have format requests
		CTRL.proj_linear[1] = true				# To force pad=0 and julia memory (no dup) in image_init()
		transpcmap!(data, true)
	elseif (isa(data, GMTdataset) || isa(data, Vector{<:GMTdataset}))
		opt_T = " -Td"
		cmd, = parse_bo(d, cmd)					# Write to binary file
	elseif (isa(data, GMTcpt))
		opt_T = " -Tc"
	elseif (isa(data, GMTps))
		opt_T = " -Tp"
	elseif (isa(data, Array{UInt8}))
		fmt = parse_grd_format(d)				# See if we have format requests
		if (fmt == "")							# If no format, write a dataset
			opt_T = " -Td"
			cmd, = parse_bo(d, cmd)				# Write to binary file
		else
			data = mat2img(data)
			fname *= fmt
			opt_T = " -Ti"
		end
	elseif (isa(data, AbstractArray))
		fmt = parse_grd_format(d)				# See if we have format requests
		if (fmt == "")							# If no format, write a dataset
			opt_T = " -Td"
			cmd, = parse_bo(d, cmd)				# Write to binary file
		else
			data = mat2grid(data)
			fname *= fmt
			opt_T = " -Tg"
		end
	else
		error("Input data of unknown data type $(typeof(data))")
	end
	cmd = cmd * opt_T

	if (opt_T == " -Ti" || opt_T == " -Tg")		# See if we have a mem layout request
		if ((val = find_in_dict(d, [:layout :mem_layout])[1]) !== nothing)
			cmd = (opt_T == " -Ti") ? cmd * " -%" * arg2str(val) : cmd * " -&" * arg2str(val)
		end
	end

	if (dbg_print_cmd(d, cmd) !== nothing)  return "gmtwrite " * fname * cmd  end

	(opt_T == " -Tg" && isa(data, GMTgrid) && (data.scale != 1 || data.offset != 0)) && (fname *= "+s$(data.scale)+o$(data.offset)")
	gmt("write " * fname * cmd, data)
	(opt_T == " -Ti") && transpcmap!(data, false)		# Reset original cmap (in case it was changed)
	return nothing
end

# -----------------------------------------------------------------------------------------------
function parse_grd_format(d::Dict)::String
# Scan options to fill any of  [=<id>][+s<scale>][+o<offset>][+n<nan>][:<driver>[/<dataType>]
# that control the grid/image output format
	out = ""
	for sym in [:id :gdal :driver]
		if (haskey(d, sym))
			if (sym == :gdal || sym == :driver)
				out = "=gd"
			else
				t = arg2str(d[sym])
				(length(t) != 2) && error(@sprintf("Format code MUST have 2 characters and not %s", t))
				out = "=" * t
			end
			break
		end
	end
	if ((val = find_in_dict(d, [:scale])[1]) !== nothing)   out *= "+s" * arg2str(val)  end
	if ((val = find_in_dict(d, [:offset])[1]) !== nothing)  out *= "+o" * arg2str(val)  end
	if ((val = find_in_dict(d, [:nan :novalue :invalid :missing])[1]) !== nothing)
		out *= "+n" * arg2str(val)
	end
	if ((val = find_in_dict(d, [:driver])[1]) !== nothing)
		out *= ":" * arg2str(val)
		((val = find_in_dict(d, [:datatype])[1]) !== nothing) && (out *= "/" * arg2str(val))
	end
	del_from_dict(d, [:id :gdal])
	return out
end

# -----------------------------------------------------------------------------------------------
function transpcmap!(I::GMTimage, toGMT::Bool=true)
	# So, the shit is in GMT when making plots it expects the colormap to be column-major (MEX inheritance)
	# but gdalwrite it expects it to be row-major. So we must do some transpose work here.
	# Not anymore (GMT6.3), now gdalwrite expects column-mjor. Shit.
	(ndims(I) != 2 || I.n_colors <= 1 || eltype(I) == UInt16) && return nothing 	# Nothing to do

	n_colors = I.n_colors
	(n_colors > 2000) && (n_colors = Int(floor(n_colors / 1000)))
	if (toGMT)
		n_cols = Int(length(I.colormap) / n_colors)
		(GMTver <= v"6.2.0") && (I.colormap = vec(collect(reshape(reshape(I.colormap, n_cols, n_colors)', n_colors * n_cols, 1))) )
		(I.n_colors < 2000 || n_cols == 4) && (I.n_colors *= 1000)	#  Because gmtwrite uses a trick to know if cmap is Mx2 or Mx4
	else
		#(I.n_colors > 2000) && (I.n_colors = div(I.n_colors, 1000))		# Revert the trick
		n_cols = Int(length(I.colormap) / n_colors)
		(GMTver <= v"6.2.0") && (I.colormap = vec(collect(reshape(reshape(I.colormap, n_colors, n_cols)', n_colors * n_cols, 1))) )
	end
	return nothing
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SoapyLoopback"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaTelecom/SoapyLoopback.git", "da9930830f6e7a97360fc959665c513a2d7aa4ea")
]

dependencies = [
    Dependency("soapysdr_jll"; compat="0.8.0")
]

# Bash recipe for building across all platforms
script = raw"""
cd SoapyLoopback
mkdir build
cd build
# This package is only used in testing so we want a Debug build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Debug \
      ..
make -j${nproc}
make install
if [[ "${target}" == *-apple-* ]]; then
    mv ${libdir}/SoapySDR/modules0.8/libsoapyloopback.so  ${libdir}/SoapySDR/modules0.8/libsoapyloopback.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(;experimental=true))
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libsoapyloopback", :librtlsdrSupport, ["lib/SoapySDR/modules0.8/"])
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
export Measurement

struct Measurement
    name::String
    fields::Dict{String,Union{String,Int,Float64}}
    tags::Dict{String,String}
    timestamp::Float64

    function Measurement(name::AbstractString, fields::Dict, tags::Dict, timestamp::Float64 = time())
        return new(
            String(name),
            Dict(String(k) => v for (k, v) in fields),
            Dict(String(k) => String(v) for (k, v) in tags),
            timestamp,
        )
    end
end

function payload(m::Measurement)
    buff = IOBuffer()
    payload!(buff, m)
    return String(take!(buff))
end

writify(x) = string(x)
writify(x::AbstractString) = string("\"", x, "\"")

function payload!(io::IO, m::Measurement)
    write(io, m.name)
    if !isempty(m.tags)
        join(io, [",$(k)=$(v)" for (k, v) in m.tags])
    end
    write(io, " ")
    join(io, ["$(k)=$(writify(v))" for (k, v) in m.fields], ",")
    write(io, " ")
    write(io, string(round(Int64, m.timestamp*1e9)))
end


text:
module TestIO

using Julio
using Test

function test_io_read()
    text = """
    line 1
    line 2
    """
    iobuf = IOBuffer(text)
    Julio.open(iobuf) do io
        @test read(io, String) == text
    end

    iobuf = IOBuffer(text)
    Julio.open(iobuf) do io
        @test readline(io) == "line 1"
        @test readline(io; keep = true) == "line 2\n"
    end
end

function test_io_write()
    iobuf = IOBuffer()
    Julio.open(iobuf) do io
        write(io, "hello")
    end
    @test String(take!(iobuf)) == "hello"
end

function test_io_read_cat()
    try
        read(`cat`)
    catch
        return @test_skip false
    end
    read_finished = Ref(false)
    output = Pipe()
    input = Pipe()
    proc = run(pipeline(`cat`; stdout = output, stdin = input); wait = false)
    close(output.in)
    ans1 = ans2 = :__notset__
    try
        Julio.open(output) do io
            write(input, "hello")
            flush(input)
            ans1 = Julio.withtimeout(0.1) do
                readline(io)
                read_finished[] = true
            end
            write(input, "\nworld")
            close(input)
            ans2 = Julio.withtimeout(60) do
                read(io, String)
            end
        end
    finally
        close(input)
        close(output)
        wait(proc)
    end
    @test ans1 === nothing
    @test !read_finished[]
    @test something(ans2) == "hello\nworld"
end

function test_open_cmd()
    Julio.open(`echo "hello"`) do output
        @test read(output, String) == "hello\n"
    end
    Julio.open(`echo "hello"`) do output
        @test readline(output) == "hello"
    end
    Julio.open(`echo "hello"`) do output
        @test readline(output; keep = true) == "hello\n"
    end
end

end  # module


text:
using ReinforcementLearningBase
using ReinforcementLearningCore
using Random
using Test
using StatsBase
using Distributions: probs
using ReinforcementLearningEnvironments
using Flux
using CUDAapi

if has_cuda()
    using CuArrays
end

@testset "ReinforcementLearningCore.jl" begin
    include("core/core.jl")
    include("components/components.jl")
    include("utils/utils.jl")
end


text:
@testset "Integer" begin
	for T in [Int, Int8, BigInt]
		@test lr_eq(T(0), (), ())
		@test lr_eq(T(1), (0,), ())
		@test lr_eq(T(-1), (), (0,))

		@test issurreal(T(1))

		@test isordinal(T(1))
		@test !isordinal(T(0))
		@test !isordinal(T(-1))
	end
end


text:
"""
    Block(inds...)

A `Block` is simply a wrapper around a set of indices or enums so that it can be used to dispatch on. By
indexing a `AbstractBlockArray` with a `Block` the a block at that block index will be returned instead of
a single element.

```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(ones(2,3), [1, 1], [2, 1])
2×2-blocked 2×3 BlockArray{Float64,2}:
 1.0  1.0  │  1.0
 ──────────┼─────
 1.0  1.0  │  1.0

julia> A[Block(1, 1)]
1×2 Array{Float64,2}:
 1.0  1.0
```
"""
struct Block{N, T}
    n::NTuple{N, T}
    Block{N, T}(n::NTuple{N, T}) where {N, T} = new{N, T}(n)
end

Block{N, T}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
Block{N}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
Block() = Block{0,Int}()
Block(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
Block{1}(n::Tuple{T}) where {T} = Block{1, T}(n)
Block{N}(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
Block(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)

@inline Block(blocks::NTuple{N, Block{1, T}}) where {N,T} = Block{N, T}(ntuple(i -> blocks[i].n[1], Val(N)))
@inline Block(::Tuple{}) = Block{0,Int}(())

# iterate and broadcast like Number
length(b::Block) = 1
size(b::Block) = ()
iterate(x::Block) = (x, nothing)
iterate(x::Block, ::Any) = nothing
isempty(x::Block) = false
broadcastable(x::Block) = x
ndims(::Type{<:Block}) = 0
ndims(::Block) = 0
eltype(::Type{B}) where B<:Block = B
getindex(B::Block, ::CartesianIndex{0}) = B

# The following code is taken from CartesianIndex
@inline (+)(index::Block{N}) where {N} = Block{N}(map(+, index.n))
@inline (-)(index::Block{N}) where {N} = Block{N}(map(-, index.n))

@inline (+)(index1::Block{N}, index2::Block{N}) where {N} =
    Block{N}(map(+, index1.n, index2.n))
@inline (-)(index1::Block{N}, index2::Block{N}) where {N} =
    Block{N}(map(-, index1.n, index2.n))
@inline min(index1::Block{N}, index2::Block{N}) where {N} =
    Block{N}(map(min, index1.n, index2.n))
@inline max(index1::Block{N}, index2::Block{N}) where {N} =
    Block{N}(map(max, index1.n, index2.n))

@inline (+)(i::Integer, index::Block) = index+i
@inline (+)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x+i, index.n))
@inline (-)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x-i, index.n))
@inline (-)(i::Integer, index::Block{N}) where {N} = Block{N}(map(x->i-x, index.n))
@inline (*)(a::Integer, index::Block{N}) where {N} = Block{N}(map(x->a*x, index.n))
@inline (*)(index::Block, a::Integer) = *(a,index)

# comparison
@inline isless(I1::Block{N}, I2::Block{N}) where {N} = Base.IteratorsMD._isless(0, I1.n, I2.n)

# conversions
convert(::Type{T}, index::Block{1}) where {T<:Number} = convert(T, index.n[1])
convert(::Type{T}, index::Block) where {T<:Tuple} = convert(T, index.n)

Int(index::Block{1}) = Int(index.n[1])
Integer(index::Block{1}) = index.n[1]
Number(index::Block{1}) = index.n[1]

# print
Base.show(io::IO, B::Block{0,Int}) = print(io, "Block()")
function Base.show(io::IO, B::Block{N,Int}) where N
    print(io, "Block($(B.n[1])")
    for n in Base.tail(B.n)
        print(io, ", $n")
    end
    print(io, ")")
end

"""
    BlockIndex{N}

A `BlockIndex` is an index which stores a global index in two parts: the block
and the offset index into the block.

It can be used to index into `BlockArrays` in the following manner:

```jldoctest; setup = quote using BlockArrays end
julia> arr = Array(reshape(1:25, (5,5)));

julia> a = PseudoBlockArray(arr, [3,2], [1,4])
2×2-blocked 5×5 PseudoBlockArray{Int64,2}:
 1  │   6  11  16  21
 2  │   7  12  17  22
 3  │   8  13  18  23
 ───┼────────────────
 4  │   9  14  19  24
 5  │  10  15  20  25

julia> a[BlockIndex((1,2), (1,2))]
11

julia> a[BlockIndex((2,2), (2,3))]
20
```
"""
struct BlockIndex{N}
    I::NTuple{N, Int}
    α::NTuple{N, Int}
end

@inline BlockIndex(a::NTuple{N,Block{1}}, b::Tuple) where N = BlockIndex(Int.(a), b)

@inline BlockIndex(a::Int, b::Int) = BlockIndex((a,), (b,))
@inline BlockIndex(a::Tuple, b::Int) = BlockIndex(a, (b,))
@inline BlockIndex(a::Int, b::Tuple) = BlockIndex((a,), b)

@inline BlockIndex(a::Block, b::Tuple) = BlockIndex(a.n, b)
@inline BlockIndex(a::Block, b::Int) = BlockIndex(a, (b,))

@generated function BlockIndex(I::NTuple{N, Int}, α::NTuple{M, Int}) where {M,N}
    @assert M < N
    α_ex = Expr(:tuple, [k <= M ? :(α[$k]) : :(1) for k = 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds α2 = $α_ex
        BlockIndex(I, α2)
    end
end

block(b::BlockIndex) = Block(b.I...)
blockindex(b::BlockIndex{1}) = b.α[1]

BlockIndex(indcs::NTuple{N,BlockIndex{1}}) where N = BlockIndex(block.(indcs), blockindex.(indcs))

function Base.show(io::IO, B::BlockIndex)
    show(io, Block(B.I...))
    print(io, "[$(B.α[1])")
    for α in Base.tail(B.α)
        print(io, ", $α")
    end
    print(io, "]")
end

##
# checkindex
##

@inline checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::Block{N}) where N = blockcheckbounds(Bool, A, I.n...)
@inline function checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::BlockIndex{N}) where N
    bl = block(I)
    checkbounds(Bool, A, bl) || return false
    B = A[bl]
    checkbounds(Bool, B, blockindex(I)...)
end

checkbounds(::Type{Bool}, A::AbstractArray{<:Any,N}, I::AbstractVector{BlockIndex{N}}) where N =
    all(checkbounds.(Bool, Ref(A), I))

struct BlockIndexRange{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{BlockIndex{N},N}
    block::Block{N,Int}
    indices::R
end

"""
    BlockIndexRange(block, startind:stopind)

represents a cartesian range inside a block.
"""
BlockIndexRange

BlockIndexRange(block::Block{N}, inds::NTuple{N,AbstractUnitRange{Int}}) where {N} =
    BlockIndexRange{N,typeof(inds)}(inds)
BlockIndexRange(block::Block{N}, inds::Vararg{AbstractUnitRange{Int},N}) where {N} =
    BlockIndexRange(block,inds)

block(R::BlockIndexRange) = R.block

getindex(B::Block{N}, inds::Vararg{Int,N}) where N = BlockIndex(B,inds)
getindex(B::Block{N}, inds::Vararg{AbstractUnitRange{Int},N}) where N = BlockIndexRange(B,inds)
getindex(B::Block{1}, inds::Colon) = B
getindex(B::Block{1}, inds::Base.Slice) = B

getindex(B::BlockIndexRange{1}, kr::AbstractUnitRange{Int}) = BlockIndexRange(B.block, B.indices[1][kr])
getindex(B::BlockIndexRange{N}, inds::Vararg{Int,N}) where N = B.block[Base.reindex(B.indices, inds)...]

eltype(R::BlockIndexRange) = eltype(typeof(R))
eltype(::Type{BlockIndexRange{N}}) where {N} = BlockIndex{N}
eltype(::Type{BlockIndexRange{N,R}}) where {N,R} = BlockIndex{N}
IteratorSize(::Type{<:BlockIndexRange}) = Base.HasShape{1}()

first(iter::BlockIndexRange) = BlockIndex(iter.block.n, map(first, iter.indices))
last(iter::BlockIndexRange)  = BlockIndex(iter.block.n, map(last, iter.indices))

@inline function iterate(iter::BlockIndexRange)
    iterfirst, iterlast = first(iter), last(iter)
    if any(map(>, iterfirst.α, iterlast.α))
        return nothing
    end
    iterfirst, iterfirst
end
@inline function iterate(iter::BlockIndexRange, state)
    nextstate = BlockIndex(state.I, inc(state.α, first(iter).α, last(iter).α))
    nextstate.α[end] > last(iter.indices[end]) && return nothing
    nextstate, nextstate
end

size(iter::BlockIndexRange) = map(dimlength, first(iter).α, last(iter).α)
length(iter::BlockIndexRange) = prod(size(iter))

Block(bs::BlockIndexRange) = bs.block

function Base.show(io::IO, B::BlockIndexRange)
    show(io, Block(B))
    print(io, "[")
    show(io, B.indices[1])
    for α in Base.tail(B.indices)
        print(io, ", ")
        show(io, α)
    end
    print(io, "]")
end

# #################
# # support for pointers
# #################
#
# function unsafe_convert(::Type{Ptr{T}},
#                         V::SubArray{T, N, BlockArray{T, N, AT}, NTuple{N, BlockSlice{Block{1,Int}}}}) where AT <: AbstractArray{T, N} where {T,N}
#     unsafe_convert(Ptr{T}, parent(V).blocks[Int.(Block.(parentindices(V)))...])
# end

"""
    BlockSlice(indices)

Represent an AbstractUnitRange of indices that attaches a block.

Upon calling `to_indices()`, Blocks are converted to BlockSlice objects to represent
the indices over which the Block spans.

This mimics the relationship between `Colon` and `Base.Slice`.
"""
struct BlockSlice{BB,INDS<:AbstractUnitRange{Int}} <: AbstractUnitRange{Int}
    block::BB
    indices::INDS
end

Block(bs::BlockSlice{<:Block}) = bs.block

for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length,
          :unsafe_length, :start)
    @eval $f(S::BlockSlice) = $f(S.indices)
end

getindex(S::BlockSlice, i::Integer) = getindex(S.indices, i)
getindex(S::BlockSlice{<:Block}, k::AbstractUnitRange{Int}) = BlockSlice(S.block[k],S.indices[k])
getindex(S::BlockSlice{<:BlockIndexRange}, k::AbstractUnitRange{Int}) = BlockSlice(S.block[k],S.indices[k])
show(io::IO, r::BlockSlice) = print(io, "BlockSlice(", r.block, ",", r.indices, ")")

Block(bs::BlockSlice{<:BlockIndexRange}) = Block(bs.block)

struct BlockRange{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{Block{N,Int},N}
    indices::R
    BlockRange{N,R}(inds::R) where {N,R} = new{N,R}(inds)
end

# The following is adapted from Julia v0.7 base/multidimensional.jl
# definition of CartesianRange

# deleted code that isn't used, such as 0-dimensional case
"""
    BlockRange(startblock, stopblock)

represents a cartesian range of blocks.

The relationship between `Block` and `BlockRange` mimicks the relationship between
`CartesianIndex` and `CartesianRange`.
"""
BlockRange

BlockRange(inds::NTuple{N,AbstractUnitRange{Int}}) where {N} =
    BlockRange{N,typeof(inds)}(inds)
BlockRange(inds::Vararg{AbstractUnitRange{Int},N}) where {N} =
    BlockRange(inds)

(:)(start::Block{1}, stop::Block{1}) = BlockRange((first(start.n):first(stop.n),))
(:)(start::Block, stop::Block) = throw(ArgumentError("Use `BlockRange` to construct a cartesian range of blocks"))
Base.BroadcastStyle(::Type{<:BlockRange{1}}) = DefaultArrayStyle{1}()
broadcasted(::DefaultArrayStyle{1}, ::Type{Block}, r::AbstractUnitRange) = BlockRange((r,))
broadcasted(::DefaultArrayStyle{1}, ::Type{Int}, block_range::BlockRange{1}) = first(block_range.indices)
broadcasted(::DefaultArrayStyle{0}, ::Type{Int}, block::Block{1}) = Int(block)

# AbstractArray implementation
axes(iter::BlockRange{N,R}) where {N,R} = map(axes1, iter.indices)
Base.IndexStyle(::Type{BlockRange{N,R}}) where {N,R} = IndexCartesian()
@inline function Base.getindex(iter::BlockRange{N,<:NTuple{N,Base.OneTo}}, I::Vararg{Integer, N}) where {N}
    @boundscheck checkbounds(iter, I...)
    Block(I)
end
@inline function Base.getindex(iter::BlockRange{N,R}, I::Vararg{Integer, N}) where {N,R}
    @boundscheck checkbounds(iter, I...)
    Block(I .- first.(axes1.(iter.indices)) .+ first.(iter.indices))
end

@inline function iterate(iter::BlockRange)
    iterfirst, iterlast = first(iter), last(iter)
    if any(map(>, iterfirst.n, iterlast.n))
        return nothing
    end
    iterfirst, iterfirst
end
@inline function iterate(iter::BlockRange, state)
    nextstate = Block(inc(state.n, first(iter).n, last(iter).n))
    nextstate.n[end] > last(iter.indices[end]) && return nothing
    nextstate, nextstate
end

# increment & carry
@inline inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline inc(state::Tuple{Integer}, start::Tuple{Integer}, stop::Tuple{Integer}) = (state[1]+1,)
@inline function inc(state, start, stop)
    if state[1] < stop[1]
        return (state[1]+1,tail(state)...)
    end
    newtail = inc(tail(state), tail(start), tail(stop))
    (start[1], newtail...)
end

# 0-d cartesian ranges are special-cased to iterate once and only once
iterate(iter::BlockRange{0}, done=false) = done ? nothing : (Block(), true)

size(iter::BlockRange) = map(dimlength, first(iter).n, last(iter).n)
dimlength(start, stop) = stop-start+1

length(iter::BlockRange) = prod(size(iter))

first(iter::BlockRange) = Block(map(first, iter.indices))
last(iter::BlockRange)  = Block(map(last, iter.indices))

@inline function in(i::Block{N}, r::BlockRange{N}) where {N}
    _in(true, i.n, first(r).n, last(r).n)
end
_in(b, ::Tuple{}, ::Tuple{}, ::Tuple{}) = b
@inline _in(b, i, start, stop) = _in(b & (start[1] <= i[1] <= stop[1]), tail(i), tail(start), tail(stop))

# We sometimes need intersection of BlockRange to return a BlockRange
intersect(a::BlockRange{1}, b::BlockRange{1}) = BlockRange(intersect(a.indices[1], b.indices[1]))

# needed for scalar-like broadcasting

BlockSlice{Block{1,BT},RT}(a::Base.OneTo) where {BT,RT<:AbstractUnitRange} =
    BlockSlice(Block(convert(BT, 1)), convert(RT, a))::BlockSlice{Block{1,BT},RT}
BlockSlice{BlockRange{1,Tuple{BT}},RT}(a::Base.OneTo) where {BT<:AbstractUnitRange,RT<:AbstractUnitRange} =
    BlockSlice(BlockRange(convert(BT, Base.OneTo(1))), convert(RT, a))::BlockSlice{BlockRange{1,Tuple{BT}},RT}
BlockSlice{BlockIndexRange{1,Tuple{BT}},RT}(a::Base.OneTo) where {BT<:AbstractUnitRange,RT<:AbstractUnitRange} =
    BlockSlice(BlockIndexRange(Block(1), convert(BT, Base.OneTo(1))), convert(RT, a))::BlockSlice{BlockIndexRange{1,Tuple{BT}},RT}


text:
import Random: rand

using Random, MultiFloats
using Random: AbstractRNG, SamplerType
using Base: IEEEFloat

Random.rand(rng::AbstractRNG, ::SamplerType{MultiFloat{T,N}}) where {T<:IEEEFloat,N} =
    renormalize(MultiFloat(ntuple(i -> rand(rng, T) * eps(T)^(i-1), Val{N}())))


text:
dice5() = rand(1:5)

function dice7()
    r = 5*dice5() + dice5() - 6
    r < 21 ? (r%7 + 1) : dice7()
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode