text:
# TODO:
# x midpoint rule
# x scattered/interpolation
# x GCV
# x L-curve
# - bounds on λ
# - initial guess for λ? will require mods to RegularizationTools
# - scaled λ? will need to work out equivalency with λ² formulation, and resolve with
#   derivative rather than difference matrix
# - optimize λ via standard deviation?
# - relative weights?
# - arbitrary weighting -- implemented but not yet tested
# - midpoint rule with scattered?
# x add argument types for `RegularizationSmooth` constructor methods (why isn't this done
#   for the other interpolaters?)
# - make use of `munge_data` features (allow for matrix rather than vector u & t arguments?)
# - validate data and t̂
# x unit tests

const LA = LinearAlgebra

### Regularization data smoothing and interpolation
struct RegularizationSmooth{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
    u::uType
    û::uType
    t::tType
    t̂::tType
    wls::uType
    wr::uType
    d::Int              # derivative degree used to calculate the roughness
    λ::AbstractFloat    # regularization parameter
    alg::Symbol  # how to determine λ: `:fixed`, `:gcv_svd`, `:gcv_tr`, `L_curve`
    Aitp::AbstractInterpolation{FT,T}
    RegularizationSmooth{FT}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp) where FT =
        new{typeof(u),typeof(t),FT,eltype(u)}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp)
end

"""
# Arguments
- `u::Vector`:  dependent data
- `t::Vector`:  independent data

# Optional Arguments
- `t̂::Vector`: t-values to use for the smooth curve (useful when data has missing values or
               is "scattered"); if not provided, then `t̂ = t`; must be monotonically
               increasing
- `wls::{Vector,Symbol}`: weights to use with the least-squares fitting term; if set to
                          `:midpoint`, then midpoint-rule integration weights are used for
                          _both_ `wls` and `wr`
- `wr::Vector`: weights to use with the roughness term
- `d::Int = 2`: derivative used to calculate roughness; e.g., when `d = 2`, the 2nd
                derivative (i.e. the curvature) of the data is used to calculate roughness.

# Keyword Arguments
- `λ::{Number,Tuple} = 1.0`: regulariation parameter; larger values result in a smoother
                             curve; the provided value is used directly when `alg = :fixed`;
                             otherwise it is used as an initial guess for the optimization
                             method, or as bounds if a 2-tuple is provided (TBD)
- `alg::Symbol = :gcv_svd`: algorithm for determining an optimal value for λ; the provided λ
                            value is used directly if `alg = :fixed`; otherwise `alg =
                            [:gcv_svd, :gcv_tr, :L_curve]` is passed to the
                            RegularizationTools solver

## Example Constructors
Smoothing using all arguments
```julia
A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
                              wls::AbstractVector, wr::AbstractVector, d::Int=2;
                              λ::Real=1.0, alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    M = _mapping_matrix(t̂, t)
    Wls½ = LA.diagm(sqrt.(wls))
    Wr½ = LA.diagm(sqrt.(wr))
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp)
end
"""
Direct smoothing, no `t̂` or weights
```julia
A = RegularizationSmooth(u, t, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int=2; λ::Real=1.0,
                              alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    t̂ = t
    N = length(t)
    M = Array{Float64}(LA.I, N, N)
    Wls½ = Array{Float64}(LA.I, N, N)
    Wr½ = Array{Float64}(LA.I, N-d, N-d)
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp)
end
"""
`t̂` provided, no weights
```julia
A = RegularizationSmooth(u, t, t̂, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
                              d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    N, N̂ = length(t), length(t̂)
    M = _mapping_matrix(t̂, t)
    Wls½ = Array{Float64}(LA.I, N, N)
    Wr½ = Array{Float64}(LA.I, N̂-d, N̂-d)
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp)
end
"""
`t̂` and `wls` provided
```julia
A = RegularizationSmooth(u, t, t̂, wls, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
                              wls::AbstractVector, d::Int=2; λ::Real=1.0,
                              alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    N, N̂ = length(t), length(t̂)
    M = _mapping_matrix(t̂, t)
    Wls½ = LA.diagm(sqrt.(wls))
    Wr½ = Array{Float64}(LA.I, N̂-d, N̂-d)
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp)
end
"""
`wls` provided, no `t̂`
```julia
A = RegularizationSmooth(u, t, nothing, wls,d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
                              wls::AbstractVector, d::Int=2; λ::Real=1.0,
                              alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    t̂ = t
    N = length(t)
    M = Array{Float64}(LA.I, N, N)
    Wls½ = LA.diagm(sqrt.(wls))
    Wr½ = Array{Float64}(LA.I, N-d, N-d)
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp)
end
"""
`wls` and `wr` provided, no `t̂`
```julia
A = RegularizationSmooth(u, t, nothing, wls, wr, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
                              wls::AbstractVector, wr::AbstractVector, d::Int=2;
                              λ::Real=1.0, alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    t̂ = t
    N = length(t)
    M = Array{Float64}(LA.I, N, N)
    Wls½ = LA.diagm(sqrt.(wls))
    Wr½ = LA.diagm(sqrt.(wr))
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp)
end
"""
Keyword provided for `wls`, no `t̂`
```julia
A = RegularizationSmooth(u, t, nothing, :midpoint, d; λ=[1.0], alg=[:gcv_svd])
```
"""
function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
                              wls::Symbol, d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd)
    u, t = munge_data(u, t)
    t̂ = t
    N = length(t)
    M = Array{Float64}(LA.I, N, N)
    wls, wr = _weighting_by_kw(t, d, wls)
    Wls½ = LA.diagm(sqrt.(wls))
    Wr½ = LA.diagm(sqrt.(wr))
    û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg)
    RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp)
end
# """ t̂ provided and keyword for wls  _TBD_ """
# function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
#                               wls::Symbol, d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd)

""" Solve for the smoothed depedent variables and create spline interpolator """
function _reg_smooth_solve(u::AbstractVector, t̂::AbstractVector, d::Int, M::AbstractMatrix,
                           Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol)
    λ = float(λ) # `float` expected by RT
    D = _derivative_matrix(t̂,d)
    Ψ = RT.setupRegularizationProblem(Wls½*M, Wr½*D)
    Wls½u = Wls½*u

    if alg == :fixed
        b̄ = RT.to_standard_form(Ψ, Wls½u) # via b̄
        ū = RT.solve(Ψ, b̄, λ)
        û = RT.to_general_form(Ψ, Wls½u, ū)
    else
        # note that the provided λ is currently ignored; using as initial guess or bounds
        # TBD, JJS 12/21/21
        result = RT.solve(Ψ, Wls½u, alg=alg)
        û = result.x
        λ = result.λ
    end
    Aitp = CubicSpline(û,t̂)
    # It seems logical to use B-Spline of order d+1, but I am unsure if theory supports the
    # extra computational cost, JJS 12/25/21
    #Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average)
    return û, λ, Aitp
end

"""
Order d derivative matrix for the provided t vector
"""
function _derivative_matrix(t::AbstractVector, d::Int)
    N = length(t)
    if d == 0
        return Array{Float64}(LA.I, (N, N))
    end
    dt = t[d+1:end] - t[1:end-d]
    V = LA.diagm(1 ./ dt)
    Ddm1 = diff(_derivative_matrix(t, d-1), dims=1)
    D = d * V*Ddm1
    return D
end

"""Linear interpolation mapping matrix, which maps `û` to `u`."""
function _mapping_matrix(t̂::AbstractVector,t::AbstractVector)
    N = length(t)
    N̂ = length(t̂)
    # map the scattered points to the appropriate index of the smoothed points
    idx = searchsortedlast.(Ref(t̂),t)
    # allow for "extrapolation"; i.e. for t̂ extremum that are interior to t
    idx[idx.==0] .+= 1
    idx[idx.==N̂] .+= -1
    # create the linear interpolation matrix
    m2 = @. (t - t̂[idx])/(t̂[idx+1] - t̂[idx])
    M = zeros(eltype(t), (N,N̂))
    for i in 1:N
        M[i,idx[i]] = 1 - m2[i]
        M[i,idx[i]+1] = m2[i]
    end
    return M
end

""" Common-use weighting, currently only `:midpoint` for midpoint-rule integration """
function _weighting_by_kw(t::AbstractVector, d::Int, wls::Symbol)
    # `:midpoint` only for now, but plan to add functionality for `:relative` weighting
    N = length(t)
    if wls == :midpoint
        bmp = zeros(N)
        bmp[1] = -t[1] + t[2]
        for i in 2:N-1
            bmp[i] = -t[i-1] + t[i+1]
        end
        bmp[N] = -t[N-1] + t[N]
        # divide by 2 doesn't matter in the minimize step, but keeping for correctness if
        # used as a template elsewhere
        bmp = bmp/2
        start = floor(Int, d/2) + 1
        final = iseven(d) ? N - (start - 1)  : N - start
        b̃mp = bmp[start:final]
        return bmp, b̃mp
    else
        throw("Unknown `$(wls)` keyword used for weighting, use `:midpoint`.")
    end
end

_interpolate(A::RegularizationSmooth{<:AbstractVector{<:Number}}, t::Number) =
    _interpolate(A.Aitp, t)


text:
# xn = a11xn_1 + a12xn_2 + a13xn_3 (mod m1)
# Yn = a21Yn-i + a22Yn-2 + a23 (mod m2)
# zn = xn - Yn (mod m1)
# un = zn/m1
# Period : 2^191

import Dates

a11 = 0
a12 = 1403580
a13 = -810728
m1  = 2^32 - 209
a21 = 527612
a22 = 0
a23 = -1370589
m2  = 2^32 - 22853

function usage()
    println("Usage: julia " * PROGRAM_FILE * " <lower-limit> <upper-limit> <how-many>")
    exit(1)
end

macro assert(ex, msgs...)
    msg_body = isempty(msgs) ? ex : msgs[1]
    msg      = string(msg_body)
    return :($ex ? nothing : throw(AssertionError($msg)))
end

function validate_args()
    ll = nothing
    ul = nothing
    n  = nothing
    try
        ll = parse(Int64,ARGS[1])
        ul = parse(Int64,ARGS[2])
        n  = parse(Int64,ARGS[3])
        if n > ul-ll+1
            error("LimitsMismatch")
        end
        if ll > ul
            error("LimitsMismatch")
        end
    catch e
        println(e)
        exit(1)
    end
    return (ll,ul,n)
end

ticks() = (return convert(UInt64,Dates.datetime2unix(Dates.unix2datetime(time()))*1000);)

function generate_initial_xs()
    x1 = 2096367628361 + Int(ticks())
    x2 = 1638267496885 + Int(ticks())
    x3 = 6377396958484 + Int(ticks())
    y1 = 9542607296802 + Int(ticks())
    y2 = 7716177293560 + Int(ticks())
    return (x1,x2,x3,y1,y2)
end

function generate_next(xn_1,xn_2,xn_3,Yn_1,Yn_2,m1,m2)
    xn = mod(a11 * xn_1 + a12 * xn_2 + a13 * xn_3, m1)
    Yn = mod(a21 * Yn_1 + a22 * Yn_2 + a23, m2)
    zn = mod(xn - Yn, m1)
    un = zn/m1
    return (xn,Yn,un)
end

function squeeze_range(ll,ul,r,float_mode)
    val = ll + (ul - ll + 1) * r
    if float_mode
        return val
    end
    return trunc(Int,val)
end

function main()
    try
        @assert(length(ARGS) == 3 , "")
    catch
        usage()
    end
    ll, ul, n          = validate_args()
    x0, x1, x2, y1, y2 = generate_initial_xs()
    i = 0
    d = Dict()
    while i<n
        xn,Yn,r = generate_next(x0,x1,x2,y1,y2,m1,m2)
        ran_num = squeeze_range(ll,ul,r,false)
        x0 = x1
        x1 = x2
        x2 = xn
        y1 = y2
        y2 = Yn
        if !(string(ran_num) in keys(d))
            println(ran_num)
            d[string(ran_num)] = nothing
            i=i+1
        end
    end
end

main()

ext: jl

max_stars_repo_name: mannasoumya/RALOTA

max_issues_repo_name: mannasoumya/RALOTA

max_forks_repo_name: mannasoumya/RALOTA


text:
module JuliaProjectFormatterTests

using JuliaProjectFormatter
using Test

function withfiles(f, tree)
    mktempdir() do dir
        for (path, content) in tree
            fullpath = joinpath(dir, path)
            mkpath(dirname(fullpath))
            open(fullpath; write = true) do io
                write(io, content)
            end
        end
        f(dir)
    end
end

function test_ignore_git()
    bad = """
    c = 3
    b = 2
    a = 1
    """

    good = """
    a = 1
    b = 2
    c = 3
    """

    tree = [
        # path => content
        ".git/Project.toml" => bad,
        "Project.toml" => bad,
        "sub/Project.toml" => bad,
    ]
    withfiles(tree) do dir
        format_projects(dir)
        @test Text(read(joinpath(dir, ".git/Project.toml"), String)) == Text(bad)
        @test Text(read(joinpath(dir, "Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/Project.toml"), String)) == Text(good)
    end
end

function test_juliaproject()
    bad = """
    c = 3
    b = 2
    a = 1
    """

    good = """
    a = 1
    b = 2
    c = 3
    """

    tree = [
        # path => content
        "JuliaProject.toml" => bad,
        "sub/Project.toml" => bad,
    ]
    withfiles(tree) do dir
        format_projects(dir)
        @test Text(read(joinpath(dir, "JuliaProject.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/Project.toml"), String)) == Text(good)
    end
end

function test_subsubsub()
    bad = """
    c = 3
    b = 2
    a = 1
    """

    good = """
    a = 1
    b = 2
    c = 3
    """

    tree = [
        # path => content
        "Project.toml" => bad,
        "sub/Project.toml" => bad,
        "sub/sub/Project.toml" => bad,
        "sub/sub/sub/Project.toml" => bad,
    ]
    withfiles(tree) do dir
        format_projects(dir)
        @test Text(read(joinpath(dir, "Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/sub/Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/sub/sub/Project.toml"), String)) == Text(good)
    end
end

function test_other_files()
    bad = """
    c = 3
    b = 2
    a = 1
    """

    good = """
    a = 1
    b = 2
    c = 3
    """

    tree = [
        # path => content
        "Project.toml" => bad,
        "sub/Project.toml" => bad,
        "other_file.txt" => "other file with invalid TOML",
        "sub/other_file.txt" => "other file with invalid TOML",
    ]
    withfiles(tree) do dir
        format_projects(dir)
        @test Text(read(joinpath(dir, "Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/Project.toml"), String)) == Text(good)
    end
end

function test_real_project()
    bad = """
    uuid = "139e0cf6-64b2-48ae-acd2-4285c0bb47c3"
    authors = ["Takafumi Arakaki <aka.tkf@gmail.com> and contributors"]
    version = "0.1.0-DEV"
    name = "JuliaProjectFormatter"

    [deps]
    TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
    [compat]
    julia = "1.6"
    """

    good = """
    name = "JuliaProjectFormatter"
    uuid = "139e0cf6-64b2-48ae-acd2-4285c0bb47c3"
    authors = ["Takafumi Arakaki <aka.tkf@gmail.com> and contributors"]
    version = "0.1.0-DEV"

    [deps]
    TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

    [compat]
    julia = "1.6"
    """

    tree = [
        # path => content
        "Project.toml" => bad,
        "sub/Project.toml" => bad,
        "other_file.txt" => "other file with invalid TOML",
        "sub/other_file.txt" => "other file with invalid TOML",
    ]
    withfiles(tree) do dir
        format_projects(dir)
        @test Text(read(joinpath(dir, "Project.toml"), String)) == Text(good)
        @test Text(read(joinpath(dir, "sub/Project.toml"), String)) == Text(good)
    end
end

end  # module JuliaProjectFormatterTests


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

# using Markdown
# using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ f5486f67-7bfc-44e2-91b9-9401d81666da
#=╠═╡ notebook_exclusive
begin
	using PlutoDevMacros
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ e3d5c718-d98c-4d53-8fc9-911be34c9f2d
#=╠═╡ notebook_exclusive
using BenchmarkTools
  ╠═╡ notebook_exclusive =#

# ╔═╡ 47c42d27-88f1-4a27-bda9-54a2439b09a1
#=╠═╡ notebook_exclusive
include("basics.jl")
  ╠═╡ notebook_exclusive =#

# ╔═╡ fcbd82ae-c04d-4f87-bbb7-5f73bdbf8bd0
html"""

Disclaimer

The code in this notebook is made to be viewed with a fork of Pluto that provides functionality to make cells exclusive to the notebook (meaning that they are commented out in the .jl file).

This is a custom feature I use to clean up notebooks and only execute the relevant cells of a notebook when this is included from normal julia (I use notebooks as building blocks for packages)

This is heavily inspired by the cell disabling that exists in Pluto and the source code to save notebook exclusivity on the file is copied/adapted from pull request #1209.
The actual modifications to achieve this functionalities
are shown here

When opening this notebook without that functionality, all cells after the macro and functions definition are notebook_exclusive and are thus surrounded by block comments.

To try out the exclusive parts of the notebook, press this <button>button</button> toggle between commenting in or out the cells by removing (or adding) the leading and trailing block comments from the cells that are marked as notebook_exclusive.

You will then have to use Ctrl-S to execute all modified cells (where the block comments were removed)

You still need to use at least version 0.17 of Pluto as the @plutoinclude macro only works properly with the macro analysis functionality that was added in that version (PlutoHooks)

The automatic reload of the macro when re-executing the cell is broken with CM6 so the whole cell should add/delete empty spaces after the macro before re-executing

"""

# ╔═╡ 5089d8dd-6587-4172-9ffd-13cf43e8c341
#=╠═╡ notebook_exclusive
md"""
## Main Functions
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ b87d12be-a37b-4202-9426-3eef14d8253c
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol("#plutoinclude_",basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
			 :(using PlutoDevMacros: @plutoinclude), # This is needed for nested @plutoinclude calls
             :(include($path))))
	m
end

# ╔═╡ 57efc195-6a2f-4ad3-94fd-53e884838789
#=╠═╡ notebook_exclusive
md"""
# Other Ingredients Helpers
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ aa28b5d8-e0d7-4b97-9220-b61a0c5f4fc4
html_reload_button() = html"""

	Reload @pluto_include

"""

# ╔═╡ 98b1fa0d-fad1-4c4f-88a0-9452d492c4cb
function include_expr(from::Module,kwargstrs::String...; to::Module)
	modname = Symbol("#plutoincluded_module")
	ex = Expr(:block, :(const $modname = $from))
	kwargs = (Symbol(s) => true for s ∈ kwargstrs if s ∈ ("all","imported"))
	varnames = names(from;kwargs...)
	# Remove the symbols that start with a '#' (still to check what is the impact)
	filter!(!Base.isgensym,varnames)
	# Symbols to always exclude from imports
	exclude_names = (
			nameof(from),
			:eval,
			:include,
			Symbol("@bind"),
			Symbol("@plutoinclude"), # Since we included this in the module
		)
	for s ∈ varnames
		if s ∉ exclude_names
			if getfield(from,s) isa Function
				# _copymethods!(ex, s; to, from, importedlist = varnames, fromname = modname)
				ret_types = Base.return_types(getfield(from,s))
				candidate_type = ret_types[1]
				# Get the eventual docstring
				docstring = Base.doc(Base.Docs.Binding(from, s))
				if all(x -> x === candidate_type, ret_types) && Base.isconcretetype(candidate_type)
					push!(ex.args, :(@doc ($docstring) $s(args...; kwargs...)::$candidate_type = $modname.$s(args...; kwargs...)))
				else
					push!(ex.args, :(@doc ($docstring) $s(args...; kwargs...) = $modname.$s(args...; kwargs...)))
				end
			else
				push!(ex.args,:(const $s = $modname.$s))
			end
		end
	end
	# Add the html to re-run the cell
	push!(ex.args,:($(html_reload_button())))
	ex
end

# ╔═╡ 872bd88e-dded-4789-85ef-145f16003351
"""
	@plutoinclude path nameskwargs...
	@plutoinclude modname=path namekwargs...

This macro is used to include external julia files inside a pluto notebook and is inspired by the discussion on [this Pluto issue](https://github.com/fonsp/Pluto.jl/issues/1101).

It requires Pluto >= v0.17.0 and includes and external file, taking care of putting in the caller namespace all varnames that are tagged with `export varname` inside the included file.

The macro relies on the use of [`names`](@ref) to get the variable names to be exported, and support providing the names of the keyword arguments of `names` to be set to true as additional strings

When called from outside Pluto, it simply returns nothing
"""
macro plutoinclude(ex,kwargstrs...)
	path = ex isa String ? ex : Base.eval(__module__,ex)
	if is_notebook_local(__source__.file::Symbol |> String)
		# If this is called directly from the notebook, do the hack to export the various variables from the module
		m = ingredients(path)
		esc(include_expr(m,kwargstrs...; to = __module__))
	elseif first(nameof(__module__) |> String, 13) == "#plutoinclude"
		# We are in a chained plutoinclude, simply include the subfile
		:(include($path)) |> esc
	else
		# We are not in the notebook and not in a chained include, so do nothing
		nothing
	end
end

# ╔═╡ 748b8eab-2f3d-4afd-bfb4-fee3240d391b
export @plutoinclude

# ╔═╡ 1f291bd2-9ab1-4fd2-bf50-49253726058f
#=╠═╡ notebook_exclusive
md"""
## Example Use
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ cf0d13ea-7562-4b8c-b7e6-fb2f1de119a7
#=╠═╡ notebook_exclusive
md"""
The cells below assume to also have the test notebook `ingredients_include_test.jl` from PlutoUtils in the same folder, download it and put it in the same folder in case you didn't already
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ bd3b021f-db44-4aa1-97b2-04002f76aeff
#=╠═╡ notebook_exclusive
notebook_path = "./plutoinclude_test.jl"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0e3eb73f-091a-4683-8ccb-592b8ccb1bee
#=╠═╡ notebook_exclusive
md"""
Try changing the content of the included notebook by removing some exported variables and re-execute (**using Shift-Enter**) the cell below containing the @plutoinclude call to see that variables are correctly updated.

You can also try leaving some variable unexported and still export all that is defined in the notebook by using
```julia
@plutoinclude notebook_path "all"
```

Finally, you can also assign the full imported module in a specific variable by doing
```julia
@plutoinclude varname = notebook_path
```
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ d2ac4955-d2a0-48b5-afcb-32baa59ade21
#=╠═╡ notebook_exclusive
@plutoinclude notebook_path "all"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0d1f5079-a886-4a07-9e99-d73e0b8a2eec
#=╠═╡ notebook_exclusive
@macroexpand @plutoinclude notebook_path "all"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 61924e22-f052-43a5-84b1-5512d222af26
#=╠═╡ notebook_exclusive
@benchmark asd(TestStruct())
  ╠═╡ notebook_exclusive =#

# ╔═╡ 50759ca2-45ca-4005-9182-058a5cb68359
#=╠═╡ notebook_exclusive
const mm = ingredients(notebook_path)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 4cec781b-c6d7-4fd7-bbe3-f7db0f973698
#=╠═╡ notebook_exclusive
a
  ╠═╡ notebook_exclusive =#

# ╔═╡ a7e7123f-0e7a-4771-9b9b-d0da97fefcef
#=╠═╡ notebook_exclusive
b
  ╠═╡ notebook_exclusive =#

# ╔═╡ 2c41234e-e1b8-4ad8-9134-85cd65a75a2d
#=╠═╡ notebook_exclusive
c
  ╠═╡ notebook_exclusive =#

# ╔═╡ ce2a2025-a6e0-44ab-8631-8d308be734a9
#=╠═╡ notebook_exclusive
d
  ╠═╡ notebook_exclusive =#

# ╔═╡ d8be6b4c-a02b-43ec-b176-de6f64fefd87
#=╠═╡ notebook_exclusive
# Extending the method
asd(s::String) = "STRING"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8090dd72-a47b-4d9d-85df-ceb0c1bcedf5
#=╠═╡ notebook_exclusive
asd(3)
  ╠═╡ notebook_exclusive =#

# ╔═╡ d1fbe484-dcd0-456e-8ec1-c68acd708a08
#=╠═╡ notebook_exclusive
asd(TestStruct())
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8df0f262-faf2-4f99-98e2-6b2a47e5ca31
#=╠═╡ notebook_exclusive
asd(TestStruct(),3,4)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 1754fdcf-de3d-4d49-a2f0-9e3f4aa3498e
#=╠═╡ notebook_exclusive
asd("S")
  ╠═╡ notebook_exclusive =#

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
PlutoDevMacros = "a0499f29-c39b-4c5c-807c-88074221b949"

[compat]
BenchmarkTools = "~1.2.0"
PlutoDevMacros = "~0.3.6"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[deps.PlutoDevMacros]]
deps = ["MacroTools", "PlutoHooks"]
git-tree-sha1 = "5d6d3c0f37bd5a635e0795943e33fc863b430035"
uuid = "a0499f29-c39b-4c5c-807c-88074221b949"
version = "0.3.6"

[[deps.PlutoHooks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "f297787f7d7507dada25f6769fe3f08f6b9b8b12"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
"""

# ╔═╡ Cell order:
# ╠═f5486f67-7bfc-44e2-91b9-9401d81666da
# ╠═e3d5c718-d98c-4d53-8fc9-911be34c9f2d
# ╠═748b8eab-2f3d-4afd-bfb4-fee3240d391b
# ╠═47c42d27-88f1-4a27-bda9-54a2439b09a1
# ╟─fcbd82ae-c04d-4f87-bbb7-5f73bdbf8bd0
# ╟─5089d8dd-6587-4172-9ffd-13cf43e8c341
# ╠═b87d12be-a37b-4202-9426-3eef14d8253c
# ╟─57efc195-6a2f-4ad3-94fd-53e884838789
# ╠═98b1fa0d-fad1-4c4f-88a0-9452d492c4cb
# ╠═872bd88e-dded-4789-85ef-145f16003351
# ╠═aa28b5d8-e0d7-4b97-9220-b61a0c5f4fc4
# ╟─1f291bd2-9ab1-4fd2-bf50-49253726058f
# ╟─cf0d13ea-7562-4b8c-b7e6-fb2f1de119a7
# ╠═bd3b021f-db44-4aa1-97b2-04002f76aeff
# ╟─0e3eb73f-091a-4683-8ccb-592b8ccb1bee
# ╠═d2ac4955-d2a0-48b5-afcb-32baa59ade21
# ╠═0d1f5079-a886-4a07-9e99-d73e0b8a2eec
# ╠═8090dd72-a47b-4d9d-85df-ceb0c1bcedf5
# ╠═61924e22-f052-43a5-84b1-5512d222af26
# ╠═50759ca2-45ca-4005-9182-058a5cb68359
# ╠═4cec781b-c6d7-4fd7-bbe3-f7db0f973698
# ╠═a7e7123f-0e7a-4771-9b9b-d0da97fefcef
# ╠═2c41234e-e1b8-4ad8-9134-85cd65a75a2d
# ╠═ce2a2025-a6e0-44ab-8631-8d308be734a9
# ╠═d1fbe484-dcd0-456e-8ec1-c68acd708a08
# ╠═8df0f262-faf2-4f99-98e2-6b2a47e5ca31
# ╠═d8be6b4c-a02b-43ec-b176-de6f64fefd87
# ╠═1754fdcf-de3d-4d49-a2f0-9e3f4aa3498e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
Meta.parse("""
    { "Hello": "World" }
""")


text:
module TestAdd

module Wrapper
    using CxxCall
    dir = mktempdir()
    lib = joinpath(dir, "libAddCxx")
    filepath = joinpath(dir, "AddCxx.cxx")
    eval(cxxsetup())
    eval(cxxnewfile(filepath,
    """
    #include <iostream>
    """))

    @cxx lib function add(x::Cint, y::Cint)::Cfloat
        """
        float ret = x + y;
        std::cout << "Welcome to libAddCxx" << std::endl;
        std::cout << "x=" << x << " y=" << y << " ret=" << ret << std::endl;
        return ret;
        """
    end

    @cxx lib function add_sloppy(x::Convert(Any=>Cint), y::Convert(Int64=>Cint))::Cfloat
        """
        float ret = x + y;
        return ret;
        """
    end
end#module Wrapper

using Test
import .Wrapper
using Libdl

@testset "AddCxx" begin
    @test !ispath(Wrapper.filepath)
    Wrapper.cxx_write_code!()
    @test isfile(Wrapper.filepath)
    libpath = "$(Wrapper.lib).$dlext"
    run(`g++ -shared -fPIC $(Wrapper.filepath) -o $libpath`)
    @test isfile(libpath)
    dlopen(libpath)
    @test Wrapper.add(Cint(1),Cint(2)) === 3f0
    @test_throws MethodError Wrapper.add(Int64(1),Cint(2)) === 3f0
    @test Wrapper.add_sloppy(Int128(1), Int64(2)) === 3f0
    @test Wrapper.add_sloppy(Float16(1), Int64(2)) === 3f0
    @test_throws MethodError Wrapper.add_sloppy(Float16(1), Int32(2)) === 3f0
end

end#module


text:
include("../src/SoapySDR.jl")
using Printf
using FFTW
using PyPlot
using DSP

include("fmDemod.jl")

# enumerate devices
(kwargs, sz) = SoapySDR.SoapySDRDevice_enumerate()

t = unsafe_load(kwargs)
for i = 1:Int(sz[])
    @printf "\nnumber of results in device = %i\n" t.size
    @printf "Found device #%d: " i
    keys = unsafe_string.(unsafe_wrap(Array, t.keys, t.size))
    vals = unsafe_string.(unsafe_wrap(Array, t.vals, t.size))
    for j = 1:t.size
        @printf "%s=%s, " keys[j] vals[j]
    end
    @printf "\n"
end

@printf "\nnumber of devices = %i\n" Int(sz[])
SoapySDR.SoapySDRKwargs_clear(kwargs)

# create device instance
# args can be user defined or from the enumeration result
sdr = SoapySDR.SoapySDRDevice_make(kwargs)

if (unsafe_load(sdr) == C_NULL)
    @printf "SoapySDRDevice_make fail: %s\n" unsafe_string(SoapySDR.SoapySDRDevice_lastError())
end

# query device info
(name, sz) = SoapySDR.SoapySDRDevice_listAntennas(sdr, SoapySDR.SOAPY_SDR_RX, 0)
@printf "Rx antennas: "
for i=1:Int(sz[])
    @printf "%s, " unsafe_string.(unsafe_wrap(Array, name, Int(sz[])))[i]
end
@printf "\n"

(name2, sz2) = SoapySDR.SoapySDRDevice_listGains(sdr, SoapySDR.SOAPY_SDR_RX, 0)
@printf "Rx gains: "
for i=1:Int(sz[])
    @printf "%s, " unsafe_string.(unsafe_wrap(Array, name2, Int(sz2[])))[i]
end
@printf "\n"

(ranges, sz) = SoapySDR.SoapySDRDevice_getFrequencyRange(sdr, SoapySDR.SOAPY_SDR_RX, 0)

@printf "Rx freq ranges: "
for i = 1:Int(sz[])
    range =  unsafe_wrap(Array, ranges, Int(sz[]))[i]
    @printf "[%g Hz -> %g Hz], " range.minimum range.maximum
end
@printf "\n"

# apply settings
#sampRate = 1024e3
sampRate = 2048e3
#sampRate = 512e3
if (SoapySDR.SoapySDRDevice_setSampleRate(sdr, SoapySDR.SOAPY_SDR_RX, 0, sampRate) != 0)
    @printf "setSampleRate fail: %s\n" unsafe_string(SoapySDR.SoapySDRDevice_lastError())
end

#f0 = 103.3e6
f0 = 104.1e6
#f0 = 938.2e6
if (SoapySDR.SoapySDRDevice_setFrequency(sdr, SoapySDR.SOAPY_SDR_RX, 0, f0) != 0)
    @printf "setFrequency fail: %s\n" unsafe_string(SoapySDR.SoapySDRDevice_lastError())
end

# set up a stream (complex floats)
rxStream = SoapySDR.SoapySDRStream()
if (SoapySDR.SoapySDRDevice_setupStream(sdr, rxStream, SoapySDR.SOAPY_SDR_RX, SoapySDR.SOAPY_SDR_CF32, C_NULL, 0) != 0)
    @printf "setupStream fail: %s\n" unsafe_string(SoapySDR.SoapySDRDevice_lastError())
end

# start streaming
SoapySDR.SoapySDRDevice_activateStream(sdr, rxStream, 0, 0, 0)

# create a re-usable buffer for rx samples
buffsz = 1024
buff = Array{ComplexF32}(undef, buffsz)

# receive some samples
timeS = 15
timeSamp = Int(floor(timeS * sampRate / buffsz))
storeIq = zeros(ComplexF32, buffsz*timeSamp)

flags = Ref{Cint}()
timeNs = Ref{Clonglong}()
buffs = [buff]

#storeFft = zeros(timeSamp, buffsz)
storeBuff = zeros(ComplexF32,timeSamp, buffsz)
for i=1:timeSamp
    SoapySDR.SoapySDRDevice_readStream(sdr, rxStream, buffs, buffsz, flags, timeNs, 100000)
    local storeBuff[i,:] = buff
end

b = blackman(20)
storeFft = zeros(timeSamp, buffsz)
for i = 1:size(storeBuff)[1]
    local storeFft[i,:] = 20 .*log10.(abs.(fftshift(fft(storeBuff[i,:]))))
end

# get IQ array
storeIq = Array(reshape(storeBuff', :, size(storeBuff)[1]*size(storeBuff)[2])')[:]

# shutdown the stream
SoapySDR.SoapySDRDevice_deactivateStream(sdr, rxStream, 0, 0)  # stop streaming
SoapySDR.SoapySDRDevice_closeStream(sdr, rxStream)
SoapySDR.SoapySDRDevice_unmake(sdr)

function plotTimeFreq(storeFft, fs, f0)
    w, h = figaspect(0.5)
    figure(figsize=[w,h])
    minF = (f0 - fs/2 )./ 1e6
    maxF = (f0 + fs/2 )./ 1e6
    maxT = 1/fs * size(storeFft)[1] * size(storeFft)[2]
    imshow(storeFft, extent=[minF,maxF,0,maxT], aspect="auto")
    xlabel("Frequnecy (MHz)")
    ylabel("Time (s)")
    cb = colorbar()
    cb[:set_label]("Magnitude (dB)")
end

function plotTime(data, fs)
    w, h = figaspect(0.25)
    figure(figsize=[w,h])
    maxT = 1/fs * size(data)[1]
    t = range(0, maxT, length=size(data)[1])
    plot(t, data, linewidth=0.1)
    xlim([0, maxT])
    xlabel("Time (s)")
    ylabel("Amplitude")
end

function plotIq(data)
    figure()
    scatter(real.(data), imag.(data), s=10)
    xlabel("In-phase")
    ylabel("Quadrature")
    title("IQ data")
end
plotIq(storeIq[1:10000])
plotTimeFreq(storeFft, sampRate, f0)

(data, fs) =  fmDemod(storeIq, sampRate)
plotTime(data, fs)

wavwrite(data, "demod.wav", Fs=fs)


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

# --------------------------------------------------------
# Pupil
# --------------------------------------------------------
mutable struct Pupil
    _tr::Transform
    _size::Float64

    function Pupil(tr::Transform, size::Float64)
        return new(tr, size)
    end
end

tr(p::Pupil) = p._tr
text(p::Pupil) = p._size
pupilsize(p::Pupil) = p._size

# --------------------------------------------------------
# Eye
# --------------------------------------------------------
mutable struct Eye
    _tr::Transform
    _text::String
    _pupil::Pupil

    function Eye(tr::Transform, text::String)
        pupil = Pupil(Transform(SVector(0.0, 0.0, 8)), 4.0)

        return new(tr, text, pupil)
    end
end

tr(e::Eye) =  e._tr
text(e::Eye) =  e._text
pupil(e::Eye) =  e._pupil
function tr!(e::Eye, t::Transform)
    e._tr = t
end
function lookAt!(e::Eye, dir::Vec3)
    tr!(e, Transform(origin(tr(e)), dir))
end

# --------------------------------------------------------
# Head
# --------------------------------------------------------
mutable struct Head
    _tr::Transform
    _eyes::Vector{Eye}

    function Head(tr::Transform)
        x_dist = 30.0
        y_dist = 3.0
        z_dist = 4.0
        left_eye = Eye(Transform(Vec3(-x_dist, y_dist, z_dist)), "Left Eye")
        right_eye = Eye(Transform(Vec3(+x_dist, y_dist, z_dist)), "Right Eye")
        return new(tr, [left_eye, right_eye])
    end
end

tr(h::Head) =  h._tr
eyes(h::Head) =  h._eyes
eye(h::Head, lr::Symbol) = h._eyes[eye_index(lr)]
eye_index(lr::Symbol) = lr == :left ? 1 : 2

function tr(h::Head, part::Symbol)
    if (part == :head)
        return tr(h)
    elseif (part == :left_eye)
        return tr(h, :head) * tr(eyes(h)[1])
    elseif (part == :right_eye)
        return tr(h, :head) * tr(eyes(h)[2])
    elseif (part == :left_pupil)
        return return tr(h, :left_eye) * tr(pupil(eyes(h)[1]))
    elseif (part == :right_pupil)
        return tr(h, :right_eye) * tr(pupil(eyes(h)[2]))
    else
        @error "unrecognized part in Head::tr [$part]"
    end
end

# -------------------------------------------------------

__face_model = nothing
function get_face_model()
    global __face_model
    if (__face_model === nothing)
        fn = joinpath(dirname(@__FILE__), "face_model.obj")
        __face_model = FileIO.load(fn)
    end
    return __face_model
end

function draw_arrow(scene::Makie.LScene, from::Vec3, to::Vec3; scale=1.0, color=Makie.MakieCore.Automatic())
    scale = scale * 0.5
    to2 = to - from
    Makie.arrows!(scene,
        [from[1]], [from[2]], [from[3]],
        [to2[1]], [to2[2]], [to2[3]],
        linewidth=scale*0.25,
        arrowsize=scale,
        arrowcolor=color,
        linecolor=color,
        transparency=true,
    )
end

function quat(tr::Transform)
    rot, _, _  =  decomposeRTS(tr)
	w = sqrt(1.0 + tr[1,1] + tr[2,2] + tr[3,3]) / 2.0
	w4 = (4.0 * w)
	x = (tr[3,2] - tr[2,3]) / w4
	y = (tr[1,3] - tr[3,1]) / w4
	z = (tr[2,1] - tr[1,2]) / w4
    return Makie.Quaternionf0(x, y, z, w)
end

function draw_local_frame(scene::Makie.LScene, tr::Transform, caption::String = "???"; scale=5.0)
    o = zero(Vec3)
    len = 6
    ux = unitX3() * len
    uy = unitY3() * len
    uz = unitZ3() * len
    l2w = local2world(tr)

    wo = l2w * o
    wux = l2w * ux
    wuy = l2w * uy
    wuz = l2w * uz

    Makie.scatter!(scene, [wo[1]], [wo[2]], [wo[3]], color=RGBA(0.9, 0.9, 0.1, 0.2), markersize=20, transparency=true)
    draw_arrow(scene, wo, wux, scale=scale, color=RGBA(0.9, 0.1, 0.1, 0.2))
    draw_arrow(scene, wo, wuy, scale=scale, color=color=RGBA(0.1, 0.9, 0.1, 0.2))
    draw_arrow(scene, wo, wuz, scale=scale, color=color=RGBA(0.1, 0.1, 0.9, 0.2))

    q = quat(tr)
    text_pos = l2w * (o + Vec3(0.0, 0.01, 0.0))
    Makie.text!( scene,
        caption,
        rotation = q,
        position = text_pos,
        color = RGBA(0.1, 0.1, 0.1, 0.2),
        align = (:left, :baseline),
        textsize = scale *2.0/5.0,
        space = :data
    )

end

# transform a mesh coordinates
function transform_mesh(mesh::GeometryBasics.Mesh, tr::Transform)
    tr = convert(Transform{Float32}, tr)

    vertices = GeometryBasics.coordinates(mesh)
    faces = GeometryBasics.faces(mesh)
    vertices = [GeometryBasics.Point3f0(tr * Vec3(p...)) for p in vertices]
    return GeometryBasics.Mesh(vertices, faces)
end


text:
function ck_generate_test()
    X, y, mask = load_image_dataset(:ck_max)
    rbm = RBM(Degenerate, Bernoulli, 1, 1)
    h5open(joinpath(MODEL_DIR, "rbm_512.h5")) do h5
        load_params(h5, rbm, "model")
    end
    Xg = generate(rbm, X)
    imgs = [to_image(Xg[:,i], mask) for i in 1:size(Xg, 2)]
    nviewall(imgs[1:36])
end

function lfw_rbm()
    imgs =
end


text:
@testset "KNNClassifier" begin
    # create something rather silly: 3 massive clusters very separated
    rng = StableRNG(10)
    n, p = 50, 3
    x, w = make_x_and_w(n, p, rng=rng)
    w1 = rand(rng, n-5)

    y1 = fill("A", n)
    y2 = fill("B", n)
    y3 = fill("C", n)
    y = categorical(vcat(y1, y2, y3))

    # generate test data
    ntest =  5
    xtest = make_xtest(ntest, p, rng=rng)
    ytest1 = fill("A", ntest)
    ytest2 = fill("B", ntest)
    ytest3 = fill("C", ntest)
    ytest = categorical(vcat(ytest1, ytest2, ytest3))

    # Test model fit
    knnc = KNNClassifier(weights=Inverse())
    f, _, _ = fit(knnc, 1, x, y)
    f2, _ , _ = fit(knnc, 1, x, y, w)
    @test_throws DimensionMismatch fit(knnc, 1, x, y, w1)
    @test f2 isa KNNResult
    @test f2.sample_weights == w
    @test fitted_params(knnc, f) == (tree=f.tree,)

    # Check predictions gotten from the model
    p = predict(knnc, f, xtest)
    p2 = predict(knnc, f2, xtest)
    @test p[1] isa UnivariateFinite
    @test p2[1] isa UnivariateFinite
    p = predict_mode(knnc, f, xtest)
    p2 = predict_mode(knnc, f2, xtest)
    @test accuracy(p, ytest) == 1.0
    @test accuracy(p2, ytest) == 1.0
end

@testset "MultitargetKNNClassifier" begin
    # create something rather silly: 3 massive clusters very separated
    rng = StableRNG(10)
    n, p = 50, 3
    x, w = make_x_and_w(n, p, rng=rng)
    y11 = fill("A", n)
    y21 = fill("B", n)
    y31 = fill("C", n)
    y12 = fill("D", n)
    y22 = fill("E", n)
    y32 = fill("F", n)
    y1 = categorical(vcat(y11, y21, y31))
    y2 = categorical(vcat(y12, y22, y32), ordered=true)
    y = (a=y1, b=y2)

    # generate test data
    ntest =  5
    xtest = make_xtest(ntest, p, rng=rng)
    ytest11 = fill("A", ntest)
    ytest21 = fill("B", ntest)
    ytest31 = fill("C", ntest)
    ytest12 = fill("D", ntest)
    ytest22 = fill("E", ntest)
    ytest32 = fill("F", ntest)
    ytest1 = categorical(vcat(ytest11, ytest21, ytest31))
    ytest2 = categorical(vcat(ytest12, ytest22, ytest32), ordered=true)
    ytest = (a=ytest1, b=ytest2)

    # Create two models, the first has an `output_type` of `DictTable` while the second
    # has an `output_type` of `ColumnTable`
    multi_knnc = MultitargetKNNClassifier(weights=Inverse())
    multi_knnc2 = MultitargetKNNClassifier(
        weights=Inverse(), output_type=ColumnTable
    )

    # Test model fit
    f, _, _ = fit(multi_knnc, 1, x, y)
    f2, _ , _ = fit(multi_knnc, 1, x, y, w)
    f3, _, _ = fit(multi_knnc2, 1, x, y)
    f4, _ , _ = fit(multi_knnc2, 1, x, y, w)
    @test f2 isa KNNResult
    @test f2.sample_weights == w
    @test f4.sample_weights == w
    @test fitted_params(multi_knnc, f) == (tree=f.tree,)

    # Check predictions gotten from the model
    p = predict(multi_knnc, f, xtest)
    p2 = predict(multi_knnc, f2, xtest)
    p3 = predict(multi_knnc2, f3, xtest)
    p4 = predict(multi_knnc2, f4, xtest)
    p5 = predict_mode(multi_knnc, f, xtest)
    p6 = predict_mode(multi_knnc, f2, xtest)
    p7 = predict_mode(multi_knnc2, f3, xtest)
    p8 = predict_mode(multi_knnc2, f4, xtest)
    for col in [:a, :b]
        @test p[col][1] isa UnivariateFinite
        @test p2[col][1] isa UnivariateFinite
        @test p3[col][1] isa UnivariateFinite
        @test p4[col][1] isa UnivariateFinite

        @test accuracy(p5[col], ytest[col]) == 1.0
        @test accuracy(p6[col], ytest[col]) == 1.0
        @test accuracy(p7[col], ytest[col]) == 1.0
        @test accuracy(p8[col], ytest[col]) == 1.0
    end
end

# the following test is a little more rigorous:
@testset "classifier sample weights" begin
    rng = StableRNG(200)
    # assign classes a, b and c randomly to 10N points on the interval:
    N = 80
    X = (x = rand(rng, 10N), );
    y = categorical(rand(rng, "abc", 10N));
    model = KNNClassifier(K=N)

    # define sample weights corresponding to class weights 2:4:1 for
    # a:b:c:
    w = map(y) do η
        if η == 'a'
            return 2
        elseif η == 'b'
            return 4
        else
            return 1
        end
    end

    f, _, _ = MLJBase.fit(model, 1, X, y, w)
    posterior3 = average([predict(model, f, X)...])

    # skewed weights gives similarly skewed posterior:
    @test abs(pdf(posterior3, 'b')/(2*pdf(posterior3, 'a'))  - 1) < 0.1
    @test abs(pdf(posterior3, 'b')/(4*pdf(posterior3, 'c'))  - 1) < 0.1
end

@testset "KNNRegressor" begin
    # create something rather silly: 3 massive clusters very separated
    rng = StableRNG(50)
    n, p = 50, 3
    x, w = make_x_and_w(n, p, rng=rng)
    y1 = fill( 0.0, n)
    y2 = fill( 2.0, n)
    y3 = fill(-2.0, n)
    y = vcat(y1, y2, y3)

    # Test model fit
    knnr = KNNRegressor(weights=Inverse())
    f, _, _ = fit(knnr, 1, x, y)
    f2, _, _ = fit(knnr, 1, x, y, w)
    @test f2 isa KNNResult
    @test fitted_params(knnr, f) == (tree=f.tree,)

    # Create test data
    ntest =  5
    xtest = make_xtest(ntest, p, rng=rng)

    # Check predictions gotten from the model
    p = predict(knnr, f, xtest)
    @test p isa AbstractVector
    p2 = predict(knnr, f2, xtest)
    @test all(p[1:ntest] .≈ 0.0)
    @test all(p[ntest+1:2*ntest] .≈ 2.0)
    @test all(p[2*ntest+1:end] .≈ -2.0)
end

@testset "MultitargetKNNRegressor" begin
    # create something rather silly: 3 massive clusters very separated
    rng = StableRNG(50)
    n, p = 50, 3
    x, w = make_x_and_w(n, p, rng=rng)
    ymat = vcat(fill( 0.0, n, 2), fill(2.0, n, 2), fill(-2.0, n, 2))
    Y = table(ymat; names = [:a, :b])

    # Test model fit
    multi_knnr = MultitargetKNNRegressor(weights=Inverse())
    f, _ , _ = fit(multi_knnr, 1, x, Y)
    f2, _, _ = fit(multi_knnr, 1, x, Y, w)
    @test fitted_params(multi_knnr, f) == (tree=f.tree,)
    @test f2 isa KNNResult
    @test f2.sample_weights == w

    # Create test data
    ntest = 5
    xtest = make_xtest(ntest, p, rng=rng)

    # Check predictions gotten from the model
    pr = predict(multi_knnr, f, xtest)
    for col in [:a, :b]
        @test all(pr[col][1:ntest] .≈ [0.0])
        @test all(pr[col][ntest+1:2*ntest] .≈ [2.0])
        @test all(pr[col][2*ntest+1:end] .≈ [-2.0])
    end
end


text:
function normalize(x)
    x / sqrt(x' * x)
end

function findPerpendicularPart(x, y)
    y - x' * y * x / (x' * x)
end

function gramSchmdit(xs...)
    qs = Array{Array{Float64,1}}(undef, length(xs))
    i = 1
    for x ∈ xs
        setindex!(qs, x,  i)
        i += 1
    end
    i = 1
    qsLen = length(qs)
    for x ∈ qs
        if i != 1
            for j in range(i, qsLen, step = 1)
                setindex!(qs, findPerpendicularPart(getindex(qs, i - 1), getindex(qs, j)), j)
            end
        end
        i += 1
    end
    i = 1
    for x ∈ qs
        setindex!(qs, normalize(x),  i)
        i += 1
    end
    return qs
end

a = [1;2;3;0]
b = [2;3;1;0]
c = [3;2;1;0]
d = [0;0;0;1]

qs = gramSchmdit(a, b, c, d)
println(qs)


text:
# GPU memory management and pooling

using Printf
using Logging

## allocation statistics

mutable struct AllocStats
  alloc_count::Int
  alloc_bytes::Int

  free_count::Int
  free_bytes::Int

  total_time::Float64
end

const alloc_stats = AllocStats(0, 0, 0, 0, 0.0)

Base.copy(alloc_stats::AllocStats) =
  AllocStats((getfield(alloc_stats, field) for field in fieldnames(AllocStats))...)

AllocStats(b::AllocStats, a::AllocStats) =
  AllocStats(
    b.alloc_count - a.alloc_count,
    b.alloc_bytes - a.alloc_bytes,
    b.free_count - a.free_count,
    b.free_bytes - a.free_bytes,
    b.total_time - a.total_time)

## CUDA allocator

@timeit_ci function actual_alloc(bytes::Integer;
                                 stream_ordered::Bool=false,
                                 stream::Union{CuStream,Nothing}=nothing)
  # try the actual allocation
  buf = try
    time = Base.@elapsed begin
      buf = @timeit_ci "Mem.alloc" begin
        Mem.alloc(Mem.Device, bytes; async=true, stream_ordered, stream)
      end
    end

    buf
  catch err
    isa(err, OutOfGPUMemoryError) || rethrow()
    return nothing
  end

  return buf
end

@timeit_ci function actual_free(buf::Mem.DeviceBuffer; stream_ordered::Bool=false,
                                stream::Union{CuStream,Nothing}=nothing)
  # free the memory
  time = Base.@elapsed begin
    @timeit_ci "Mem.free" Mem.free(buf; async=true, stream_ordered, stream)
  end

  return
end

## stream-ordered memory pool

const __stream_ordered = LazyInitialized{Vector{Bool}}()
function stream_ordered(dev::CuDevice)
  flags = get!(__stream_ordered) do
    val = Vector{Bool}(undef, ndevices())
    if version() < v"11.2" || haskey(ENV, "CUDA_MEMCHECK")
      fill!(val, false)
    else
      for dev in devices()
        val[deviceid(dev)+1] = attribute(dev, DEVICE_ATTRIBUTE_MEMORY_POOLS_SUPPORTED) == 1
      end
    end
    val
  end
  @inbounds flags[deviceid(dev)+1]
end

function allocatable_memory(dev::CuDevice)
  # NOTE: this function queries available memory, which obviously changes after we allocate.
  device!(dev) do
    available_memory()
  end
end

function reserved_memory(dev::CuDevice)
  # taken from TensorFlow's `MinSystemMemory`
  #
  # if the available memory is < 2GiB, we allocate 225MiB to system memory.
  # otherwise, depending on the capability version assign
  #  500MiB (for cuda_compute_capability <= 6.x) or
  # 1050MiB (for cuda_compute_capability <= 7.x) or
  # 1536MiB (for cuda_compute_capability >= 8.x)
  available = allocatable_memory(dev)
  if available <= 1<<31
    225 * 1024 * 1024
  else
    cap = capability(dev)
    if cap <= v"6"
      500 * 1024 * 1024
    elseif cap <= v"7"
      1050 * 1024 * 1024
    else
      1536 * 1024 * 1024
    end
  end
end

# per-device flag indicating the status of a pool
const _pool_status = PerDevice{Base.RefValue{Union{Nothing,Bool}}}()
pool_status(dev::CuDevice) = get!(_pool_status, dev) do
  # nothing=uninitialized, false=idle, true=active
  Ref{Union{Nothing,Bool}}(nothing)
end
function pool_mark(dev::CuDevice)
  status = pool_status(dev)
  if status[] === nothing
      pool = memory_pool(dev)

      # first time on this context, so configure the pool
      attribute!(memory_pool(dev), MEMPOOL_ATTR_RELEASE_THRESHOLD,
                 UInt64(reserved_memory(dev)))

      # also launch a task to periodically trim the pool
      if isinteractive() && !isassigned(__pool_cleanup)
        __pool_cleanup[] = @async pool_cleanup()
      end
  end
  status[] = true
  return
end

# reclaim unused pool memory after a certain time
const __pool_cleanup = Ref{Task}()
function pool_cleanup()
  idle_counters = fill(0, ndevices())
  while true
    for (i, dev) in enumerate(devices())
      stream_ordered(dev) || continue

      status = pool_status(dev)
      status[] === nothing && continue

      if status[]
        idle_counters[i] = 0
      else
        idle_counters[i] += 1
      end
      status[] = 0

      if idle_counters[i] == 5
        # the pool hasn't been used for a while, so reclaim unused buffers
        device!(dev) do
          reclaim()
        end
      end
    end

    sleep(60)
  end
end

## interface

export OutOfGPUMemoryError

"""
    OutOfGPUMemoryError()

An operation allocated too much GPU memory for either the system or the memory pool to
handle properly.
"""
struct OutOfGPUMemoryError <: Exception
  sz::Int

  OutOfGPUMemoryError(sz::Integer=0) = new(sz)
end

function Base.showerror(io::IO, err::OutOfGPUMemoryError)
    print(io, "Out of GPU memory")
    if err.sz > 0
      print(io, " trying to allocate $(Base.format_bytes(err.sz))")
    end
    println(io)
    memory_status(io)
end

"""
    alloc(sz)

Allocate a number of bytes `sz` from the memory pool. Returns a buffer object; may throw
an [`OutOfGPUMemoryError`](@ref) if the allocation request cannot be satisfied.
"""
@inline @timeit_ci function alloc(sz; stream::Union{Nothing,CuStream}=nothing)
  # 0-byte allocations shouldn't hit the pool
  sz == 0 && return Mem.DeviceBuffer(CU_NULL, 0)

  state = active_state()

  buf = nothing
  gctime = 0.0  # using Base.@timed/gc_num is too expensive
  time = Base.@elapsed begin
    if stream_ordered(state.device)
      # mark the pool as active
      pool_mark(state.device)

      for phase in 1:4
          if phase == 2
              gctime += Base.@elapsed GC.gc(false)
          elseif phase == 3
              gctime += Base.@elapsed GC.gc(true)
          elseif phase == 4
              device_synchronize()
          end

          buf = actual_alloc(sz; stream_ordered=true, stream=something(stream, state.stream))
          buf === nothing || break
      end
    else
      for phase in 1:4
          if phase == 2
              gctime += Base.@elapsed GC.gc(false)
          elseif phase == 3
              gctime += Base.@elapsed GC.gc(true)
          end

          buf = actual_alloc(sz; stream_ordered=false)
          buf === nothing || break
      end
    end
  end
  buf === nothing && throw(OutOfGPUMemoryError(sz))

  alloc_stats.alloc_count += 1
  alloc_stats.alloc_bytes += sz
  alloc_stats.total_time += time - gctime
  # NOTE: total_time might be an over-estimation if we trigger GC somewhere else

  return buf
end

"""
    free(buf)

Releases a buffer `buf` to the memory pool.
"""
@inline @timeit_ci function free(buf::Mem.DeviceBuffer;
                                 stream::Union{Nothing,CuStream}=nothing)
  # XXX: have @timeit use the root timer, since we may be called from a finalizer

  # 0-byte allocations shouldn't hit the pool
  sizeof(buf) == 0 && return

  state = active_state()

  # this function is typically called from a finalizer, where we can't switch tasks,
  # so perform our own error handling.
  try

    time = Base.@elapsed begin
      if stream_ordered(state.device)
        # mark the pool as active
        pool_mark(state.device)

        actual_free(buf; stream_ordered=true, stream=something(stream, state.stream))
      else
        actual_free(buf; stream_ordered=false)
      end
    end

    alloc_stats.free_count += 1
    alloc_stats.free_bytes += sizeof(buf)
    alloc_stats.total_time += time
  catch ex
    Base.showerror_nostdio(ex, "WARNING: Error while freeing $buf")
    Base.show_backtrace(Core.stdout, catch_backtrace())
    Core.println()
  end

  return
end

"""
    reclaim([sz=typemax(Int)])

Reclaims `sz` bytes of cached memory. Use this to free GPU memory before calling into
functionality that does not use the CUDA memory pool. Returns the number of bytes
actually reclaimed.
"""
function reclaim(sz::Int=typemax(Int))
  dev = device()
  if stream_ordered(dev)
      # TODO: respect sz
      device_synchronize()
      trim(memory_pool(dev))
  else
    0
  end
end

"""
    @retry_reclaim isfailed(ret) ex

Run a block of code `ex` repeatedly until it successfully allocates the memory it needs.
Retries are only attempted when calling `isfailed` with the current return value is true.
At each try, more and more memory is freed from the CUDA memory pool. When that is not
possible anymore, the latest returned value will be returned.

This macro is intended for use with CUDA APIs, which sometimes allocate (outside of the
CUDA memory pool) and return a specific error code when failing to.
"""
macro retry_reclaim(isfailed, ex)
  quote
    ret = nothing
    phase = 0
    while true
      phase += 1
      ret = $(esc(ex))
      $(esc(isfailed))(ret) || break

      # incrementally more costly reclaim of cached memory
      if stream_ordered(device())
        if phase == 1
          # synchronizing streams forces asynchronous free operations to finish.
          device_synchronize()
        elseif phase == 2
          GC.gc(false)
          device_synchronize()
        elseif phase == 3
          GC.gc(true)
          device_synchronize()
        elseif phase == 4
          # maybe this allocation doesn't use the pool, so trim it.
          reclaim()
          device_synchronize()
        else
          break
        end
      else
        if phase == 1
          GC.gc(false)
        elseif phase == 2
          GC.gc(true)
        else
          break
        end
      end
    end
    ret
  end
end

# XXX: function version for use in CUDAdrv where we haven't loaded pool.jl yet
function retry_reclaim(f, check)
  @retry_reclaim check f()
end

## utilities

used_memory(ctx=context()) = @lock allocated_lock begin
    mapreduce(sizeof∘first, +, values(allocated(ctx)); init=0)
end

"""
    @allocated

A macro to evaluate an expression, discarding the resulting value, instead returning the
total number of bytes allocated during evaluation of the expression.
"""
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = alloc_stats.alloc_bytes
                $(esc(ex))
                alloc_stats.alloc_bytes - b0
            end
            f()
        end
    end
end

"""
    @time ex

Run expression `ex` and report on execution time and GPU/CPU memory behavior. The GPU is
synchronized right before and after executing `ex` to exclude any external effects.

"""
macro time(ex)
    quote
        local val, cpu_time,
            cpu_alloc_size, cpu_gc_time, cpu_mem_stats,
            gpu_alloc_size, gpu_mem_time, gpu_mem_stats = @timed $(esc(ex))

        local cpu_alloc_count = Base.gc_alloc_count(cpu_mem_stats)
        local gpu_alloc_count = gpu_mem_stats.alloc_count

        Printf.@printf("%10.6f seconds", cpu_time)
        for (typ, gctime, memtime, bytes, allocs) in
            (("CPU", cpu_gc_time, 0, cpu_alloc_size, cpu_alloc_count),
             ("GPU", 0, gpu_mem_time, gpu_alloc_size, gpu_alloc_count))
          if bytes != 0 || allocs != 0
              allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
              if ma == 1
                  Printf.@printf(" (%d%s %s allocation%s: ", allocs, Base._cnt_units[ma], typ, allocs==1 ? "" : "s")
              else
                  Printf.@printf(" (%.2f%s %s allocations: ", allocs, Base._cnt_units[ma], typ)
              end
              print(Base.format_bytes(bytes))
              if gctime > 0
                  Printf.@printf(", %.2f%% gc time", 100*gctime/cpu_time)
              end
              if memtime > 0
                  Printf.@printf(", %.2f%% memmgmt time", 100*memtime/cpu_time)
              end
              print(")")
          else
              if gctime > 0
                  Printf.@printf(", %.2f%% %s gc time", 100*gctime/cpu_time, typ)
              end
              if memtime > 0
                  Printf.@printf(", %.2f%% %s memmgmt time", 100*memtime/cpu_time, typ)
              end
          end
        end
        println()

        val
    end
end

macro timed(ex)
    quote
        while false; end # compiler heuristic: compile this block (alter this if the heuristic changes)

        # @time(d) might surround an application, so be sure to initialize CUDA before that
        CUDA.prepare_cuda_state()

        # coarse synchronization to exclude effects from previously-executed code
        synchronize()

        local gpu_mem_stats0 = copy(alloc_stats)
        local cpu_mem_stats0 = Base.gc_num()
        local cpu_time0 = time_ns()

        # fine-grained synchronization of the code under analysis
        local val = @sync blocking=false $(esc(ex))

        local cpu_time1 = time_ns()
        local cpu_mem_stats1 = Base.gc_num()
        local gpu_mem_stats1 = copy(alloc_stats)

        local cpu_time = (cpu_time1 - cpu_time0) / 1e9

        local cpu_mem_stats = Base.GC_Diff(cpu_mem_stats1, cpu_mem_stats0)
        local gpu_mem_stats = AllocStats(gpu_mem_stats1, gpu_mem_stats0)

        (value=val, time=cpu_time,
         cpu_bytes=cpu_mem_stats.allocd, cpu_gctime=cpu_mem_stats.total_time / 1e9, cpu_gcstats=cpu_mem_stats,
         gpu_bytes=gpu_mem_stats.alloc_bytes, gpu_memtime=gpu_mem_stats.total_time, gpu_memstats=gpu_mem_stats)
    end
end

"""
    memory_status([io=stdout])

Report to `io` on the memory status of the current GPU and the active memory pool.
"""
function memory_status(io::IO=stdout)
  state = active_state()
  ctx = context()

  free_bytes, total_bytes = Mem.info()
  used_bytes = total_bytes - free_bytes
  used_ratio = used_bytes / total_bytes
  @printf(io, "Effective GPU memory usage: %.2f%% (%s/%s)\n",
              100*used_ratio, Base.format_bytes(used_bytes),
              Base.format_bytes(total_bytes))

  if stream_ordered(state.device)
    pool = memory_pool(state.device)
    if version() >= v"11.3"
      pool_reserved_bytes = attribute(UInt64, pool, MEMPOOL_ATTR_RESERVED_MEM_CURRENT)
      pool_used_bytes = attribute(UInt64, pool, MEMPOOL_ATTR_USED_MEM_CURRENT)
      @printf(io, "Memory pool usage: %s (%s reserved)",
                  Base.format_bytes(pool_used_bytes),
                  Base.format_bytes(pool_reserved_bytes))
    else
      @printf(io, "Memory pool statistics require CUDA 11.3.")
    end
  else
    @printf(io, "No memory pool is in use.")
  end
end