text:
module TestUtilities

# using Revise
using Test
using MLJ
using ScientificTypes
using Tables
using OrderedCollections

d = Dict(
    :x => Set([1, 2]),
    :y => Set([2, 3, 5]),
    :z => Set([4, 7]),
    :a => Set([8, 1]),
    :b => Set([4,]),
    :w => Set([3, 1, 2]),
    :t => Set([0,]))

dinv = Dict(
    0 => Set([:t,]),
    1 => Set([:x, :a, :w]),
    2 => Set([:x, :y, :w]),
    3 => Set([:y, :w]),
    4 => Set([:z, :b]),
    5 => Set([:y,]),
    7 => Set([:z,]),
    8 => Set([:a,]))

@test MLJ.inverse(d) == dinv

d = LittleDict(
    :x => Set([1, 2]),
    :y => Set([2, 3, 5]),
    :z => Set([4, 7]),
    :a => Set([8, 1]),
    :b => Set([4,]),
    :w => Set([3, 1, 2]),
    :t => Set([0,]))

dinv = LittleDict(
    0 => Set([:t,]),
    1 => Set([:x, :a, :w]),
    2 => Set([:x, :y, :w]),
    3 => Set([:y, :w]),
    4 => Set([:z, :b]),
    5 => Set([:y,]),
    7 => Set([:z,]),
    8 => Set([:a,]))

@test MLJ.inverse(d) == dinv

d = LittleDict()
d[:test] = Tuple{Union{Continuous,Missing}, Finite}
d["junk"] = LittleDict{Any,Any}("H" => Missing, :c => Table(Finite,Continuous),
                                :cross => "lemon", :t => :w, "r" => "r")
d["a"] = "b"
d[:f] = true
d["j"] = :post

@test MLJ.decode_dic(MLJ.encode_dic(d)) == d

end # module
true


text:
 # module ComplexDPR1

import Base: *, getindex, size, Matrix
import LinearAlgebra: eigen, eigvals, inv, transpose, mul!

T=Complex{Float64}

# export DPR1

# Define general Complex Diagonal+Rank-One Type D+ρ u*transpose(u)
struct CSymDPR1{T} <: AbstractMatrix{T}
    D::Vector{T} # diagonal
    u::Vector{T} # rank one, length n
    ρ::T  # rho
end

# Define its size
size(A::CSymDPR1, dim::Integer) = length(A.D)
size(A::CSymDPR1) = size(A,1), size(A,1)

# Index into a CSymDPR1
function getindex(A::CSymDPR1,i::Integer,j::Integer)
    Aij = A.ρ * (A.u[i] * A.u[j])
    return i == j ? A.D[i] + Aij : Aij
end # getindex

Matrix(A::CSymDPR1)= Matrix(Symmetric(BLAS.syr!('U',A.ρ,A.u,diagm(0=>A.D))))
transpose(A::CSymDPR1{T}) where T = A

# Define CSymDPR1 x vector in O(n) operations
function *(A::CSymDPR1{T}, x::Vector{T}) where T
    dotd= T==ComplexF64 ? BLAS.dotu : BLAS.dot
    y=Diagonal(A.D)*x
    BLAS.axpy!(A.ρ*dotd(A.u, x),A.u,y)
    # return y
end
function mul!(y::Vector{T}, A::CSymDPR1{T}, x::Vector{T}) where T
    dotd= T==ComplexF64 ? BLAS.dotu : BLAS.dot
    mul!(y,Diagonal(A.D),x)
    BLAS.axpy!(A.ρ*dotd(A.u, x),A.u,y)
end

# CSymDPR1 x Matrix
function *(A::CSymDPR1{T}, x::Matrix{T}) where T
    y=Diagonal(A.D)*x
    BLAS.ger!(A.ρ,A.u,conj(transpose(x)*A.u),y)
end
function mul!(y::Matrix{T}, A::CSymDPR1{T}, x::Matrix{T}) where T
    mul!(y,Diagonal(A.D),x)
    BLAS.ger!(A.ρ,A.u,conj(transpose(x)*A.u),y)
end

function inv(A::CSymDPR1{T}) where {T<:T}
    xₕ=A.u ./A.D
    Δₕ=one(T)./A.D
    ρₕ=-A.ρ/(one(T)+BLAS.dotu(xₕ,A.u)*A.ρ)
    return CSymDPR1(Δₕ,xₕ,ρₕ)
end

function inv(A::CSymDPR1{T}) where {T<:Float64}
    xₕ=A.u ./A.D
    Δₕ=one(T)./A.D
    ρₕ=-A.ρ/(one(T)+BLAS.dot(xₕ,A.u)*A.ρ)
    return CSymDPR1(Δₕ,xₕ,ρₕ)
end

function xtAx(A::CSymDPR1{T}, x::Vector{T}) where {T}
    xdx = zero(T)
    xu  = zero(T)
    @inbounds for i in 1:length(x)
        xdx += A.D[i]*x[i]^2
        xu  += x[i]*A.u[i]
    end
    return xdx + A.ρ*xu^2
end

function xtAx(A::CSymDPR1{T}, x::Vector{T}, μ::T) where {T}
    xdx = zero(T)
    xu  = zero(T)
    @inbounds for i in 1:length(x)
        xdx += (A.D[i]-μ)*x[i]^2
        xu  += x[i]*A.u[i]
    end
    return xdx + A.ρ*xu^2
end

function xtAxc(A::CSymDPR1{T}, x::Vector{T}, μ::T) where {T}
    xdx = zero(T)
    xu  = zero(T)
    xuc = zero(T)
    @inbounds for i in 1:length(x)
        xdx += (A.D[i]-μ)*x[i]*conj(x[i])
        xu  += A.u[i]*x[i]
        xuc += A.u[i]*conj(x[i])
    end
    return xdx + A.ρ*xu*xuc
end

function xtAxc(A::CSymDPR1{T}, x::Vector{T}) where {T}
    xdx = zero(T)
    xu  = zero(T)
    xuc = zero(T)
    @inbounds for i in 1:length(x)
        xdx += A.D[i]*x[i]*conj(x[i])
        xu  += A.u[i]*x[i]
        xuc += A.u[i]*conj(x[i])
    end
    return xdx + A.ρ*xu*xuc
end

function dpr1vector!(v::Vector{T}, A::CSymDPR1{T}, λ::T) where  T
    n = length(A.D)
    @inbounds for l in 1:n
        τ=A.D[l] - λ
        if τ==zero(T)
            for k=1:n
                v[k]=zero(T)
            end
            v[l]=one(T)
            return v
        end
        v[l] = A.u[l]/τ #  (A.D[l] - λ)
    end
    normalize!(v)
    # v
end

function mrqi!(z::Vector{T}, Du::Vector{T}, A::CSymDPR1{T}, λ::T) where T
    # Modified Rayleigh Quotient Iteration
    # Returns z=inv(A.D-λ)*A.u if γ is small or too large
    # and z=inv(A.D-λI)*z otherwise
    n=length(z)
    @inbounds for l in 1:n
        τ=A.D[l] - λ
        if τ==zero(T)
            for k=1:n
                z[k]=zero(T)
            end
            z[l]=one(T)
            return
        end
        Du[l] = A.u[l]/τ #  (A.D[l] - λ)
    end
    dotd= (T==Float64) ? BLAS.dot : BLAS.dotu
    γ=A.ρ*dotd(z,Du)/(one(T)+A.ρ*dotd(A.u,Du))
    if abs(γ)<1e-4 || abs(γ)>1e4 || isnan(γ)
        @inbounds for l=1:n
            z[l]=Du[l]
        end
    else
        @inbounds for l=1:n
            z[l]=z[l]/(A.D[l]-λ)-γ*Du[l]
        end
    end
    normalize!(z)
end

function rqi(z::Vector{T}, z₁::Vector{T}, τ::Vector{T}, A::CSymDPR1{T} ,λ::T) where T
    # Rayleigh Quotient Iteration
    # Returns λ=z'*A*z/z'*z
    # and z=inv(A-λI)*z otherwise
    dotd= (T==Float64) ? BLAS.dot : BLAS.dotu
    if count(isequal(λ),A.D)>0
        return λ
    end
    n=length(A.D)
    dpr1vector!(z,A,λ)
    δ = xtAxc(A, z)
    if count(isequal(δ),A.D)>0
        λ=δ
        return λ
    else
        copy!(z₁,A.u)
        copy!(τ,A.D)
        τ.-=δ
        z₁./=τ
        μ=dotd(A.u,z₁)
        γ=-A.ρ/(one(T)+μ*A.ρ)
        ξ=dotd(z₁,z)
        γ*=ξ
        z./=τ
        BLAS.axpy!(γ,z₁,z)
        return xtAxc(A, z)/BLAS.dotc(z,z)
    end
end

function eigvals(A₁::CSymDPR1)
    # A is complex and unreduced, that is, A.u.!=0 and
    # all elements of A.D are different.
    # It is better to use z^T A z / z^T z (modified RQI) instead
    # of standard RQI using z^*
    # This Standard Rayleigh Quotient Iteration, but solving the system
    # (multiplication by the inverse) is replaced by the CSymDPR1()
    # eigenvector formula

    # Preliminaries
    A=deepcopy(A₁)
    n = length(A.D)
    m = n

    T = typeof(A.D[1])

    # Eigenvalues
    λ = Vector{T}(undef,n)
    # Auxiliary arrays
    d = Vector{T}(undef,n)
    z = Vector{T}(undef,n)
    Du = Vector{T}(undef,n)
    AD=Vector{T}(undef,n)
    dotd= (T==Float64) ? BLAS.dot : BLAS.dotu
    # Tolerance
    tol = 100.0*eps()
    # Maxiter
    maxiter = 101

    for l in 1:n
        copyto!(AD,A.D)
        # Shift to the absolutely largest pole
        t=1
        μ = A.D[t]
        for i in 1:m
            A.D[i] -= μ
        end
        for i=1:m
            z[i]=zero(T)
        end
        z[t]=one(T)
        k = 1
        δ = one(T)
        γ = zero(T)
        while k < maxiter && abs(δ)>tol
            # use z^T
            ν = dotd(z, z)
            δ = xtAx(A, z)/ν
            γ += δ
            mrqi!(z, Du, A, δ)
            for i in 1:m
                A.D[i] -= δ
            end

            k += 1
        end
        if k==maxiter
            println(l," ",abs(δ)," ",tol)
        end
        λ[l] = μ + γ

        # Deflation - eliminating absolutely largest pole
        @inbounds for i in 1:(t - 1)
            d[i] = A.u[i]/sqrt(one(T) - γ/(AD[i] - μ))
            if isnan(d[i])
                println(l," nan")
                d[i]=zero(T)
            end
        end
        @inbounds for i in (t + 1):m
            d[i-1] = A.u[i]/sqrt(one(T) - γ/(AD[i] - μ))
            if isnan(d[i-1])
                d[i-1]=zero(T)
            end
        end
        m -= 1
        deleteat!(A.D,t)
        deleteat!(AD,t)
        deleteat!(z,t)
        deleteat!(Du,t)
        copyto!(A.D,AD)
        deleteat!(A.u,t)
        for i=1:m
            A.u[i]=d[i]
        end
    end
    # Two steps of RQI on the original CSymDPR1 matrix
    # Generate workspace for each thread

    nthreads=Threads.nthreads()
    # nthreads=1
    z₀ = [Vector{T}(undef,n) for i = 1:nthreads]
    z₁ = [Vector{T}(undef,n) for i = 1:nthreads]
    τ = [Vector{T}(undef,n) for i = 1:nthreads]
    # Threaded loop
    Threads.@threads for l=1:n
        tid=Threads.threadid()
        μ=λ[l]
        λ[l]=rqi(z₀[tid],z₁[tid],τ[tid],A₁,μ)
    end
    return λ
end

# Scaling
function scale!(X::CauchyLike{T}) where T
    n=length(X.x)
    dotd= (T==Float64) ? BLAS.dot : BLAS.dotu
    a=[Vector{T}(undef,n) for i=1:Threads.nthreads()]
    Ψ=Vector{T}(undef,Threads.nthreads())
    rzero=count(iszero,X.r)
    Threads.@threads for i=1:n
        tid=Threads.threadid()
        getindex!(a[tid],X,:,i)
        Ψ[tid]=dotd(a[tid],a[tid])
        if !isfinite(Ψ[tid]) || count(iszero,a[tid])>rzero
            X.s[i]=zero(T)
        else
            X.s[i]/=conj(sqrt(Ψ[tid]))
        end
    end
end

function eigen(A::CSymDPR1{T}) where T
    # A is complex and unreduced, that is, A.u.!=0 and
    # all elements of A.D are different.
    # eig() is just eigvals with normalized eigenvectors returned
    # as Cauchy-like matrix
    Λ=eigvals(A)
    n=length(A.D)
    # Boley-Golub Lemma
    z=Vector{T}(undef,n)
    Threads.@threads for i=1:n
        δ=A.D[i]
        τ=Λ[i]-δ
        for j=1:i-1
            τ*=(Λ[j]-δ)/(A.D[j]-δ)
        end
        for j=i+1:n
            τ*=(Λ[j]-δ)/(A.D[j]-δ)
        end
        τ/=A.ρ
        τ=sqrt(τ)
        if isnan(τ)
            # println(" nan ")
            z[i]=A.u[i]
        else
            if T==Float64
                z[i]=copysign(τ,A.u[i])
            else
                z[i]=copysign(real(τ),real(A.u[i]))+im*copysign(imag(τ),imag(A.u[i]))
            end
        end
    end

    nz=norm(A.u-z,Inf)
    if nz>100*eps()
        # println("warning: nz=",nz)
    end
    X=CauchyLike(A.D,Λ,z,ones(T,length(A.u)))
    scale!(X)
    return Λ,X
end

function traceX(S::CauchyLike{T},s::Int) where T
    n=div(length(S.x),2)
    G=Matrix{T}(undef,2n,2s)
    a=[Vector{T}(undef,2n) for i=1:Threads.nthreads()]
    b=[Vector{T}(undef,2n) for i=1:Threads.nthreads()]
    Threads.@threads for i=1:2n
        tid=Threads.threadid()
        getindex!(a[tid],S,:,i)
        G[i,:].=view(a[tid],[1:s;n+1:n+s])
    end
    Y=CauchyLikeS(S.y,G,conj(G))
    # trace(X)=trace(SYS')=trace(S'(SY))
    S₁=S*Y
    traceX=zeros(T,Threads.nthreads())
    Threads.@threads for i=1:2n
        tid=Threads.threadid()
        getindex!(a[tid],S,:,i)
        getindex!(b[tid],S₁,:,i)
        traceX[tid]+=BLAS.dotc(a[tid],b[tid])
    end
    return sum(traceX)
end

ext: jl


text:
function quote_string(str::AbstractString)::AbstractString
    repr(str)[2:end-1]
end

function runcmd(args::Cmd; wrap=identity)::Vector{String}
    julia = Base.julia_cmd()
    dbftp = joinpath("..", "bin", "dbftp.jl")
    lines = String[]
    cmd = `$julia $dbftp $args`
    open(wrap(cmd)) do io
        skipcount = 0
        for line in eachline(io)
            if skipcount > 0
                skipcount -= 1
                continue
            elseif startswith(line, "Julia Dropbox client")
                skipcount = 1
                continue
            elseif startswith(line, "Info: ") || startswith(line, "\rInfo: ")
                continue
            else
                push!(lines, line)
            end
        end
    end
    lines
end

@testset "Command version" begin
    lines = runcmd(`version`)
    @test length(lines) == 1
    m = match(r"^Version\s+(.*)", lines[1])
    @test m !== nothing
    version = VersionNumber(m.captures[1])
end

@testset "Option --verbose" begin
    lines = runcmd(`--verbose version`)
    @test length(lines) == 1
    m = match(r"^Version\s+(.*)", lines[1])
    @test m !== nothing
    version = VersionNumber(m.captures[1])
end

@testset "Option --debug" begin
    lines = runcmd(`--debug version`)
    @test length(lines) == 1
    m = match(r"^Version\s+(.*)", lines[1])
    @test m !== nothing
    version = VersionNumber(m.captures[1])
end

@testset "Command account" begin
    lines = runcmd(`account`)
    @test length(lines) == 1
    m = match(r"^Account: Name:\s+", lines[1])
    @test m !== nothing
end

@testset "Command du" begin
    lines = runcmd(`du`)
    @test length(lines) == 2
    m = match(r"^allocated:\s+(.*)", lines[1])
    @test m !== nothing
    m = match(r"^used:\s+(.*)", lines[2])
    @test m !== nothing
end

@testset "Command mkdir" begin
    lines = runcmd(`mkdir $folder`)
    @test length(lines) == 0
end

@testset "Command put" begin
    mktempdir() do dir
        filename = joinpath(dir, "hello")
        content = "Hello, World!\n"
        write(filename, content)
        lines = runcmd(`put $filename $folder`)
        @test length(lines) == 0

        lines = runcmd(`put $filename $folder/hello1`)
        @test length(lines) == 0

        # Upload same file again
        lines = runcmd(`put $filename $folder/hello1`)
        @test length(lines) == 0

        # Upload different file with same name
        filename1 = joinpath(dir, "hello1")
        content1 = "Hello, World 1!\n"
        write(filename1, content1)
        lines = runcmd(`put $filename1 $folder/hello1`)
        @test length(lines) == 0

        filename2 = joinpath(dir, "hello2")
        content2 = "Hello, World 2!\n"
        write(filename2, content2)
        lines = runcmd(`put $filename2 $folder/hello2`)
        @test length(lines) == 0

        dirname = joinpath(dir, "dir")
        mkdir(dirname)
        filename3 = joinpath(dirname, "hello3")
        content3 = "Hello, World 3!\n"
        write(filename3, content3)
        lines = runcmd(`put $dirname $folder`)
        @test length(lines) == 0
    end
end

@testset "Command cmp" begin
    mktempdir() do dir
        filename = joinpath(dir, "hello")
        content = "Hello, World!\n"
        write(filename, content)
        lines = runcmd(`cmp $filename $folder`)
        @test length(lines) == 0

        filename1 = joinpath(dir, "hello1")
        content1 = "Hello, World 1!\n"
        write(filename1, content1)
        lines = runcmd(`cmp $filename1 $folder/hello1`)
        @test length(lines) == 0

        lines = runcmd(`cmp $filename $folder/hello2`; wrap=ignorestatus)
        @test length(lines) == 1
        quoted_filename = quote_string(filename)
        @test lines[1] == "$quoted_filename: File size differs"

        filename2 = joinpath(dir, "hello2")
        content2 = "Hello, World 2!\n"
        write(filename2, content2)
        lines = runcmd(`cmp $filename2 $folder/hello2`)
        @test length(lines) == 0

        # TODO: compare recursively
        dirname = joinpath(dir, "dir")
        mkdir(dirname)
        filename3 = joinpath(dirname, "hello3")
        content3 = "Hello, World 3!\n"
        write(filename3, content3)
        lines = runcmd(`cmp $filename3 $folder/dir/hello3`)
        @test length(lines) == 0
    end
end

@testset "Command ls" begin
    lines = runcmd(`ls $folder`)
    @test lines == ["dir",
                    "hello",
                    "hello1",
                    "hello2",
                    ]

    lines = runcmd(`ls -R $folder`)
    @test lines == [".",
                    "dir",
                    "dir/hello3",
                    "hello",
                    "hello1",
                    "hello2",
                    ]

    lines = runcmd(`ls -l $folder`)
    @test length(lines) == 4
    @test startswith(lines[1], "d    ")
    @test startswith(lines[2], "- 14 ")
    @test startswith(lines[3], "- 16 ")
    @test startswith(lines[4], "- 16 ")
    @test endswith(lines[1], " dir")
    @test endswith(lines[2], " hello")
    @test endswith(lines[3], " hello1")
    @test endswith(lines[4], " hello2")
end

@testset "Command get" begin
    mktempdir() do dir
        filename = joinpath(dir, "hello")
        lines = runcmd(`get $folder/hello $dir`)
        @test length(lines) == 0
        content = read(filename, String)
        @test content == "Hello, World!\n"

        filename1 = joinpath(dir, "hello1")
        lines = runcmd(`get $folder/hello1 $filename1`)
        @test length(lines) == 0
        content = read(filename1, String)
        @test content == "Hello, World 1!\n"

        filename2 = joinpath(dir, "hello2")
        lines = runcmd(`get $folder/hello2 $filename2`)
        @test length(lines) == 0
        content2 = read(filename2, String)
        @test content2 == "Hello, World 2!\n"

        # TODO: download recursively
        dirname = joinpath(dir, "dir")
        mkdir(dirname)
        filename3 = joinpath(dirname, "hello3")
        lines = runcmd(`get $folder/dir/hello3 $filename3`)
        @test length(lines) == 0
        content3 = read(filename3, String)
        @test content3 == "Hello, World 3!\n"
    end
end

@testset "Command rm" begin
    lines = runcmd(`rm $folder`)
    @test length(lines) == 0
end


text:
# ---
# title: 1291. Sequential Digits
# id: problem1291
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/sequential-digits/description/>
# hidden: true
# ---
#
# An integer has _sequential digits_ if and only if each digit in the number is
# one more than the previous digit.
#
# Return a **sorted** list of all the integers in the range `[low, high]`
# inclusive that have sequential digits.
#
#
#
# **Example 1:**
#
#
#
#     Input: low = 100, high = 300
#     Output: [123,234]
#
#
# **Example 2:**
#
#
#
#     Input: low = 1000, high = 13000
#     Output: [1234,2345,3456,4567,5678,6789,12345]
#
#
#
#
# **Constraints:**
#
#   * `10 <= low <= high <= 10^9`
#
#
## @lc code=start
using LeetCode

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


text:
##############################################################################
##
## Reshaping
##
## Also, see issue # ??
##
##############################################################################

##############################################################################
##
## stack()
## melt()
##
##############################################################################

"""
Stacks a DataFrame; convert from a wide to long format

```julia
stack(df::AbstractDataFrame, [measure_vars], [id_vars];
      variable_name::Symbol=:variable, value_name::Symbol=:value)
melt(df::AbstractDataFrame, [id_vars], [measure_vars];
     variable_name::Symbol=:variable, value_name::Symbol=:value)
```

### Arguments

* `df` : the AbstractDataFrame to be stacked

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`. If neither `measure_vars`
  or `id_vars` are given, `measure_vars` defaults to all
  floating point columns.

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

* `variable_name` : the name of the new stacked column that shall hold the names
  of each of `measure_vars`

* `value_name` : the name of the new stacked column containing the values from
  each of `measure_vars`

### Result

* `::DataFrame` : the long-format DataFrame with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

See also `stackdf` and `meltdf` for stacking methods that return a
view into the original DataFrame. See `unstack` for converting from
long to wide format.

### Examples

```julia
d1 = DataFrame(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stack(d1, [:c, :d])
d1s2 = stack(d1, [:c, :d], [:a])
d1m = melt(d1, [:a, :b, :e])
d1s_name = melt(d1, [:a, :b, :e], variable_name=:somemeasure)
```

"""
function stack(df::AbstractDataFrame, measure_vars::AbstractVector{<:Integer},
               id_vars::AbstractVector{<:Integer}; variable_name::Symbol=:variable,
               value_name::Symbol=:value)
    N = length(measure_vars)
    cnames = names(df)[id_vars]
    insert!(cnames, 1, value_name)
    insert!(cnames, 1, variable_name)
    DataFrame(AbstractVector[repeat(_names(df)[measure_vars], inner=nrow(df)),   # variable
                  vcat([df[c] for c in measure_vars]...),             # value
                  [repeat(df[c], outer=N) for c in id_vars]...],      # id_var columns
              cnames)
end
function stack(df::AbstractDataFrame, measure_var::Int, id_var::Int;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(df, [measure_var], [id_var];
          variable_name=variable_name, value_name=value_name)
end
function stack(df::AbstractDataFrame, measure_vars::AbstractVector{<:Integer}, id_var::Int;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(df, measure_vars, [id_var];
          variable_name=variable_name, value_name=value_name)
end
function stack(df::AbstractDataFrame, measure_var::Int, id_vars::AbstractVector{<:Integer};
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(df, [measure_var], id_vars;
          variable_name=variable_name, value_name=value_name)
end
function stack(df::AbstractDataFrame, measure_vars, id_vars;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(df, index(df)[measure_vars], index(df)[id_vars];
          variable_name=variable_name, value_name=value_name)
end
# no vars specified, by default select only numeric columns
numeric_vars(df::AbstractDataFrame) =
    [T <: AbstractFloat || (T >: Missing && Missings.T(T) <: AbstractFloat)
     for T in eltypes(df)]

function stack(df::AbstractDataFrame, measure_vars = numeric_vars(df);
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    mv_inds = index(df)[measure_vars]
    stack(df, mv_inds, setdiff(1:ncol(df), mv_inds);
          variable_name=variable_name, value_name=value_name)
end

"""
Stacks a DataFrame; convert from a wide to long format; see
`stack`.
"""
function melt(df::AbstractDataFrame, id_vars::Union{Int,Symbol};
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    melt(df, [id_vars]; variable_name=variable_name, value_name=value_name)
end
function melt(df::AbstractDataFrame, id_vars;
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    id_inds = index(df)[id_vars]
    stack(df, setdiff(1:ncol(df), id_inds), id_inds;
          variable_name=variable_name, value_name=value_name)
end
function melt(df::AbstractDataFrame, id_vars, measure_vars;
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(df, measure_vars, id_vars; variable_name=variable_name,
          value_name=value_name)
end
melt(df::AbstractDataFrame; variable_name::Symbol=:variable, value_name::Symbol=:value) =
    stack(df; variable_name=variable_name, value_name=value_name)

##############################################################################
##
## unstack()
##
##############################################################################

"""
Unstacks a DataFrame; convert from a long to wide format

```julia
unstack(df::AbstractDataFrame, rowkeys::Union{Symbol, Integer},
        colkey::Union{Symbol, Integer}, value::Union{Symbol, Integer})
unstack(df::AbstractDataFrame, rowkeys::AbstractVector{<:Union{Symbol, Integer}},
        colkey::Union{Symbol, Integer}, value::Union{Symbol, Integer})
unstack(df::AbstractDataFrame, colkey::Union{Symbol, Integer},
        value::Union{Symbol, Integer})
unstack(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame to be unstacked

* `rowkeys` : the column(s) with a unique key for each row, if not given,
  find a key by grouping on anything not a `colkey` or `value`

* `colkey` : the column holding the column names in wide format,
  defaults to `:variable`

* `value` : the value column, defaults to `:value`

### Result

* `::DataFrame` : the wide-format DataFrame

If `colkey` contains `missing` values then they will be skipped and a warning will be printed.

If combination of `rowkeys` and `colkey` contains duplicate entries then last `value` will
be retained and a warning will be printed.

### Examples

```julia
wide = DataFrame(id = 1:12,
                 a  = repeat([1:3;], inner = [4]),
                 b  = repeat([1:4;], inner = [3]),
                 c  = randn(12),
                 d  = randn(12))

long = stack(wide)
wide0 = unstack(long)
wide1 = unstack(long, :variable, :value)
wide2 = unstack(long, :id, :variable, :value)
wide3 = unstack(long, [:id, :a], :variable, :value)
```
Note that there are some differences between the widened results above.
"""
function unstack(df::AbstractDataFrame, rowkey::Int, colkey::Int, value::Int)
    refkeycol = categorical(df[rowkey])
    droplevels!(refkeycol)
    keycol = categorical(df[colkey])
    droplevels!(keycol)
    valuecol = df[value]
    _unstack(df, rowkey, colkey, value, keycol, valuecol, refkeycol)
end

function _unstack(df::AbstractDataFrame, rowkey::Int,
                  colkey::Int, value::Int, keycol, valuecol, refkeycol)
    Nrow = length(refkeycol.pool)
    Ncol = length(keycol.pool)
    unstacked_val = [similar_missing(valuecol, Nrow) for i in 1:Ncol]
    hadmissing = false # have we encountered missing in refkeycol
    mask_filled = falses(Nrow+1, Ncol) # has a given [row,col] entry been filled?
    warned_dup = false # have we already printed duplicate entries warning?
    warned_missing = false # have we already printed missing in keycol warning?
    keycol_order = Vector{Int}(CategoricalArrays.order(keycol.pool))
    refkeycol_order = Vector{Int}(CategoricalArrays.order(refkeycol.pool))
    for k in 1:nrow(df)
        kref = keycol.refs[k]
        if kref <= 0 # we have found missing in colkey
            if !warned_missing
                @warn("Missing value in variable $(_names(df)[colkey]) at row $k. Skipping.")
                warned_missing = true
            end
            continue # skip processing it
        end
        j = keycol_order[kref]
        refkref = refkeycol.refs[k]
        if refkref <= 0 # we have found missing in rowkey
            if !hadmissing # if it is the first time we have to add a new row
                hadmissing = true
                # we use the fact that missing is greater than anything
                for i in eachindex(unstacked_val)
                    push!(unstacked_val[i], missing)
                end
            end
            i = length(unstacked_val[1])
        else
            i = refkeycol_order[refkref]
        end
        if !warned_dup && mask_filled[i, j]
            @warn("Duplicate entries in unstack at row $k for key "*
                  "$(refkeycol[k]) and variable $(keycol[k]).")
            warned_dup = true
        end
        unstacked_val[j][i] = valuecol[k]
        mask_filled[i, j] = true
    end
    levs = levels(refkeycol)
    # we have to handle a case with missings in refkeycol as levs will skip missing
    col = similar(df[rowkey], length(levs) + hadmissing)
    copyto!(col, levs)
    hadmissing && (col[end] = missing)
    df2 = DataFrame(unstacked_val, map(Symbol, levels(keycol)))
    insert!(df2, 1, col, _names(df)[rowkey])
end

unstack(df::AbstractDataFrame, rowkey::ColumnIndex,
        colkey::ColumnIndex, value::ColumnIndex) =
    unstack(df, index(df)[rowkey], index(df)[colkey], index(df)[value])

# Version of unstack with just the colkey and value columns provided
unstack(df::AbstractDataFrame, colkey::ColumnIndex, value::ColumnIndex) =
    unstack(df, index(df)[colkey], index(df)[value])

# group on anything not a key or value
unstack(df::AbstractDataFrame, colkey::Int, value::Int) =
    unstack(df, setdiff(_names(df), _names(df)[[colkey, value]]), colkey, value)

unstack(df::AbstractDataFrame, rowkeys, colkey::ColumnIndex, value::ColumnIndex) =
    unstack(df, rowkeys, index(df)[colkey], index(df)[value])

unstack(df::AbstractDataFrame, rowkeys::AbstractVector{<:Real}, colkey::Int, value::Int) =
    unstack(df, names(df)[rowkeys], colkey, value)

function unstack(df::AbstractDataFrame, rowkeys::AbstractVector{Symbol}, colkey::Int, value::Int)
    length(rowkeys) == 0 && throw(ArgumentError("No key column found"))
    length(rowkeys) == 1 && return unstack(df, rowkeys[1], colkey, value)
    g = groupby(df, rowkeys, sort=true)
    keycol = categorical(df[colkey])
    droplevels!(keycol)
    valuecol = df[value]
    _unstack(df, rowkeys, colkey, value, keycol, valuecol, g)
end

function _unstack(df::AbstractDataFrame, rowkeys::AbstractVector{Symbol},
                  colkey::Int, value::Int, keycol, valuecol, g)
    groupidxs = [g.idx[g.starts[i]:g.ends[i]] for i in 1:length(g.starts)]
    rowkey = zeros(Int, size(df, 1))
    for i in 1:length(groupidxs)
        rowkey[groupidxs[i]] .= i
    end
    df1 = df[g.idx[g.starts], g.cols]
    Nrow = length(g)
    Ncol = length(levels(keycol))
    unstacked_val = [similar_missing(valuecol, Nrow) for i in 1:Ncol]
    mask_filled = falses(Nrow, Ncol)
    warned_dup = false
    warned_missing = false
    keycol_order = Vector{Int}(CategoricalArrays.order(keycol.pool))
    for k in 1:nrow(df)
        kref = keycol.refs[k]
        if kref <= 0
            if !warned_missing
                @warn("Missing value in variable $(_names(df)[colkey]) at row $k. Skipping.")
                warned_missing = true
            end
            continue
        end
        j = keycol_order[kref]
        i = rowkey[k]
        if !warned_dup && mask_filled[i, j]
            @warn("Duplicate entries in unstack at row $k for key "*
                 "$(tuple((df[1,s] for s in rowkeys)...)) and variable $(keycol[k]).")
            warned_dup = true
        end
        unstacked_val[j][i] = valuecol[k]
        mask_filled[i, j] = true
    end
    df2 = DataFrame(unstacked_val, map(Symbol, levels(keycol)))
    hcat(df1, df2)
end

function unstack(df::AbstractDataFrame)
    Base.depwarn("In the future `unstack(df)` will call `unstack(df, :variable, :value)`." *
                 " use `unstack(df, :id, :variable, :value)` to treat `:id` as the only " *
                 "`rowkeys` column",
                 :unstack)
    unstack(df, :id, :variable, :value)
end

##############################################################################
##
## Reshaping using referencing (issue #145)
## New AbstractVector types (all read only):
##     StackedVector
##     RepeatedVector
##
##############################################################################

"""
An AbstractVector{Any} that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

### Constructor

```julia
StackedVector(d::AbstractVector...)
```

### Arguments

* `d...` : one or more AbstractVectors

### Examples

```julia
StackedVector(Any[[1,2], [9,10], [11,12]])  # [1,2,9,10,11,12]
```

"""
struct StackedVector <: AbstractVector{Any}
    components::Vector{Any}
end

function Base.getindex(v::StackedVector,i::Integer)
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    v.components[j][k]
end

Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.ndims(v::StackedVector) = 1
Base.eltype(v::StackedVector) = promote_type(map(eltype, v.components)...)
Base.similar(v::StackedVector, T::Type, dims::Union{Integer, AbstractUnitRange}...) =
    similar(v.components[1], T, dims...)

CategoricalArrays.CategoricalArray(v::StackedVector) = CategoricalArray(v[:]) # could be more efficient

"""
An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

### Constructor

```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

### Arguments

* `parent` : the AbstractVector that's repeated
* `inner` : the numer of times each element is repeated
* `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

### Examples

```julia
RepeatedVector([1,2], 3, 1)   # [1,1,1,2,2,2]
RepeatedVector([1,2], 1, 3)   # [1,2,1,2,1,2]
RepeatedVector([1,2], 2, 2)   # [1,2,1,2,1,2,1,2]
```

"""
struct RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

function Base.getindex(v::RepeatedVector, i::Integer)
    N = length(v.parent)
    idx = Base.fld1(mod1(i,v.inner*N),v.inner)
    v.parent[idx]
end

Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(v.parent)
Base.ndims(v::RepeatedVector) = 1
Base.eltype(v::RepeatedVector{T}) where {T} = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(v.parent), v.inner, v.outer)
Base.similar(v::RepeatedVector, T::Type, dims::Dims) = similar(v.parent, T, dims)
Base.unique(v::RepeatedVector) = unique(v.parent)

function CategoricalArrays.CategoricalArray(v::RepeatedVector)
    res = CategoricalArrays.CategoricalArray(v.parent)
    res.refs = repeat(res.refs, inner = [v.inner], outer = [v.outer])
    res
end

##############################################################################
##
## stackdf()
## meltdf()
## Reshaping using referencing (issue #145), using the above vector types
##
##############################################################################

"""
A stacked view of a DataFrame (long format)

Like `stack` and `melt`, but a view is returned rather than data
copies.

```julia
stackdf(df::AbstractDataFrame, [measure_vars], [id_vars];
        variable_name::Symbol=:variable, value_name::Symbol=:value)
meltdf(df::AbstractDataFrame, [id_vars], [measure_vars];
       variable_name::Symbol=:variable, value_name::Symbol=:value)
```

### Arguments

* `df` : the wide AbstractDataFrame

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

### Result

* `::DataFrame` : the long-format DataFrame with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

The result is a view because the columns are special AbstractVectors
that return indexed views into the original DataFrame.

### Examples

```julia
d1 = DataFrame(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stackdf(d1, [:c, :d])
d1s2 = stackdf(d1, [:c, :d], [:a])
d1m = meltdf(d1, [:a, :b, :e])
```

"""
function stackdf(df::AbstractDataFrame, measure_vars::AbstractVector{<:Integer},
                 id_vars::AbstractVector{<:Integer}; variable_name::Symbol=:variable,
                 value_name::Symbol=:value)
    N = length(measure_vars)
    cnames = names(df)[id_vars]
    insert!(cnames, 1, value_name)
    insert!(cnames, 1, variable_name)
    DataFrame(AbstractVector[RepeatedVector(_names(df)[measure_vars], nrow(df), 1), # variable
                             StackedVector(Any[df[c] for c in measure_vars]),       # value
                             [RepeatedVector(df[c], 1, N) for c in id_vars]...],    # id_var columns
              cnames)
end
function stackdf(df::AbstractDataFrame, measure_var::Int, id_var::Int;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdf(df, [measure_var], [id_var]; variable_name=variable_name,
            value_name=value_name)
end
function stackdf(df::AbstractDataFrame, measure_vars, id_var::Int;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdf(df, measure_vars, [id_var]; variable_name=variable_name,
            value_name=value_name)
end
function stackdf(df::AbstractDataFrame, measure_var::Int, id_vars;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdf(df, [measure_var], id_vars; variable_name=variable_name,
            value_name=value_name)
end
function stackdf(df::AbstractDataFrame, measure_vars, id_vars;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdf(df, index(df)[measure_vars], index(df)[id_vars];
            variable_name=variable_name, value_name=value_name)
end
function stackdf(df::AbstractDataFrame, measure_vars = numeric_vars(df);
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    m_inds = index(df)[measure_vars]
    stackdf(df, m_inds, setdiff(1:ncol(df), m_inds);
            variable_name=variable_name, value_name=value_name)
end

"""
A stacked view of a DataFrame (long format); see `stackdf`
"""
function meltdf(df::AbstractDataFrame, id_vars; variable_name::Symbol=:variable,
                value_name::Symbol=:value)
    id_inds = index(df)[id_vars]
    stackdf(df, setdiff(1:ncol(df), id_inds), id_inds;
            variable_name=variable_name, value_name=value_name)
end
function meltdf(df::AbstractDataFrame, id_vars, measure_vars;
                variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdf(df, measure_vars, id_vars; variable_name=variable_name,
            value_name=value_name)
end
meltdf(df::AbstractDataFrame; variable_name::Symbol=:variable, value_name::Symbol=:value) =
    stackdf(df; variable_name=variable_name, value_name=value_name)


text:
using LinearAlgebra: Symmetric, isposdef
using Distributions: MvNormal, logpdf
using Distributions

import ..kalman_filter_mv
import ..kalman_smoother_mv

struct LinearGaussianCached{T<:Any} <: AbstractLinearGaussian
    θ::T
end

function LinearGaussianCached(; kwargs...)
    LinearGaussianCached{typeof(kwargs.data)}(kwargs.data)
end

function Base.show(io::IO, ::LinearGaussianCached)
    print(io, "LinearGaussianCached()")
end

# Generative Model
function dpr_base(m::LinearGaussianCached, Zs)
    logpdf!(m.θ._dpr, MvNormal(m.θ.μ0, m.θ.Tau0), Zs)
end

function rpr_base(model::LinearGaussianCached, rng)
    rand(
        rng,
        MvNormal(model.θ.states[:, 1], Symmetric(model.θ.variances[:, :, 1])),
        model.θ.K,
    )
end

function dt_base(m::LinearGaussianCached, Zs, Zs_new, t)
    μ = m.θ.G(t) * Zs
    Σ = m.θ.Tau(t)
    if size(m.θ._dt, 1) != size(μ, 2)
        resize!(m.θ._dt, size(μ, 2))
    end
    d = logpdf!(m.θ._dt, MvNormal(Σ), Zs_new - μ)
    @assert !any(isnan.(d))
    return d
end

function dm_base(m::LinearGaussianCached, Zs, t)
    yt = view(m.θ.Y, t, :)
    if any(ismissing.(yt))
        return zeros(eltype(Zs), size(Zs, 2))
    else
        d = logpdf!(m.θ._dm, MvNormal(m.θ.Σ(t)), m.θ.F(t) * Zs .- yt)
        @assert !any(isnan.(d))
        return d
    end
end

# Inference Model
function rinit_base(m::LinearGaussianCached, rng)
    rand(rng, MvNormal(m.θ.states_sm[:, 1], Symmetric(m.θ.variances_sm[:, :, 1])), m.θ.K)
end

function dinit_base(m::LinearGaussianCached, Zs)
    logpdf!(
        m.θ._dinit,
        MvNormal(m.θ.states_sm[:, 1], Symmetric(m.θ.variances_sm[:, :, 1])),
        Zs,
    )
end

function calc_prop(m::LinearGaussianCached, Zs, t)
    A = @view(m.θ.Covs_sm[:, :, t-1]) / @view(m.θ.variances_sm[:, :, t-1])
    eps = m.θ._eps_prop .= (Zs .- @view(m.θ.states_sm[:, t-1]))
    μ = (A * eps) .+= @view(m.θ.states_sm[:, t])

    Σ = Symmetric(
        m.θ.variances_sm[:, :, t] -
        (m.θ.Covs_sm[:, :, t-1] / m.θ.variances_sm[:, :, t-1]) * m.θ.Covs_sm[:, :, t-1]',
    )

    @assert !any(isnan.(μ))
    @assert !any(isnan.(Σ))

    return μ, Σ
end

function rp_base(m::LinearGaussianCached, rng, Zs, t)
    μ, Σ = calc_prop(m, Zs, t)
    @assert isposdef(Σ)
    Zs_new = (Distributions.rand!(rng, MvNormal(Σ), m.θ._eps[1]) .+= μ)
    m.θ._eps[1] = m.θ._eps[2]
    m.θ._eps[2] = Zs_new
    return Zs_new
end

function dp_base(m::LinearGaussianCached, Zs, Zs_new, t)
    μ, Σ = calc_prop(m, Zs, t)
    @assert isposdef(Σ)
    d = logpdf!(m.θ._dp, MvNormal(Σ), Zs_new - μ)
    return d
end

function dpre_base(m::LinearGaussianCached, Zs, t)
    d_smooth = logpdf!(
        m.θ._dsmooth,
        MvNormal(m.θ.states_sm[:, t-1], Symmetric(m.θ.variances_sm[:, :, t-1])),
        Zs,
    )
    d_filt = logpdf!(
        m.θ._dfilt,
        MvNormal(m.θ.states[:, t-1], Symmetric(m.θ.variances[:, :, t-1])),
        Zs,
    )

    d = (d_smooth .-= d_filt)
    d .+= sum(@view(m.θ.lls[t:end]))
    @assert !any(isnan.(d))
    return d
end

# function make_kalman_proposal(F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)
#     ll, states, variances, lls = kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
#     states_sm, variances_sm, Covs_sm = kalman_smoother_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
#     return @NT(ll, states, variances, lls, states_sm, variances_sm, Covs_sm)
# end

function lgc_smc_model(K, F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)
    ll, states, variances, lls, states_sm, variances_sm, Covs_sm =
        make_kalman_proposal(F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)

    m = LinearGaussianCached(
        F = F,
        G = G,
        Tau = Tau,
        Σ = Σ,
        μ0 = μ0,
        Tau0 = Tau0,
        K = K,
        states = states,
        variances = variances,
        lls = lls,
        states_sm = states_sm,
        variances_sm = variances_sm,
        Covs_sm = Covs_sm,
        Y = Y,
        _eps = [
            Matrix{Float64}(undef, size(μ0, 1), K),
            Matrix{Float64}(undef, size(μ0, 1), K),
        ],
        _eps_prop = Matrix{Float64}(undef, size(μ0, 1), K),
        _dpr = Vector{Float64}(undef, K),
        _dt = Vector{Float64}(undef, K),
        _dm = Vector{Float64}(undef, K),
        _dinit = Vector{Float64}(undef, K),
        _dp = Vector{Float64}(undef, K),
        _dsmooth = Vector{Float64}(undef, K),
        _dfilt = Vector{Float64}(undef, K),
    )

    dpr(Zs) = dpr_base(m, Zs)
    dt(Zs, Zs_new, t) = dt_base(m, Zs, Zs_new, t)
    dm(Zs, t) = dm_base(m, Zs, t)

    rinit(rng) = rinit_base(m, rng)
    dinit(Zs) = dinit_base(m, Zs)
    rp(rng, Zs, t) = rp_base(m, rng, Zs, t)
    dp(Zs, Zs_new, t) = dp_base(m, Zs, Zs_new, t)

    dpre(Zs, t) = dpre_base(m, Zs, t)

    return @NT(rinit, rp, dinit, dpr, dp, dt, dm, dpre)
end


text:
#
# client.jl --
#
# Implement XPA client methods.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2020, Éric Thiébaut (https://github.com/JuliaAstro/XPA.jl).
#

"""
    XPA.Client()

yields a persistent XPA client connection which can be used for calls to
[`XPA.set`](@ref) and [`XPA.get`](@ref) methods.  Persistence means that a
connection to an XPA server is not closed when one of the above calls is
completed but will be re-used on successive calls.  Using `XPA.Client()`
therefore saves the time it takes to connect to a server, which could be
significant with slow connections or if there will be a large number of
exchanges with a given access point.

!!! note
    To avoid the delay for connecting to the XPA server, all XPA methods that
    perform XPA client requests now automatically use a connection that is kept
    open for the calling thread.  Directly calling `XPA.Client()` should be
    unnecessary, this method is kept for backward compatibility.

See also [`XPA.set`](@ref), [`XPA.get`](@ref), [`XPA.list`](@ref) and
[`XPA.find`](@ref).

"""
Client() = Client(_open())

"""
    XPA.connection()

yields a persistent XPA client connection that is kept open for the calling
thread (a different connection is memorized for each Julia thread).

Per-thread client connections are automatically open (or even re-open) as
needed.

""" connection

const CONNECTIONS = Client[]

function connection()
    id = Threads.threadid()
    while length(CONNECTIONS) < id
        push!(CONNECTIONS, Client(C_NULL))
    end
    if !isopen(CONNECTIONS[id])
        # Automatically (re)open the client connection.
        CONNECTIONS[id].ptr = _open()
    end
    return CONNECTIONS[id]
end

function _open()
    # The argument of XPAOpen is currently ignored (it is reserved for future
    # use).
    ptr = ccall((:XPAOpen, libxpa), Ptr{Cvoid}, (Ptr{Cvoid},), C_NULL)
    ptr != C_NULL || error("failed to create a persistent XPA connection")
    return ptr
end

function Base.close(conn::Client)
    if (ptr = conn.ptr) != C_NULL
        conn.ptr = C_NULL # avoid closing more than once!
        ccall((:XPAClose, libxpa), Cvoid, (Ptr{Cvoid},), ptr)
    end
    return nothing
end

"""
    XPA.list(conn=XPA.connection())

yields a list of available XPA access points.  The result is a vector of
[`XPA.AccessPoint`](@ref) instances.  Optional argument `conn` is a persistent
XPA client connection (created by [`XPA.Client`](@ref)); if omitted, a
per-thread connection is used (see [`XPA.connection`](@ref)).

See also [`XPA.Client`](@ref), [`XPA.connection`](@ref) and [`XPA.find`](@ref).

"""
function list(conn::Client = connection())
    lst = AccessPoint[]
    for str in split(chomp(get(String, conn, "xpans")), r"\n|\r\n?";
                     keepempty=false)
        arr = split(str; keepempty=false)
        if length(arr) != 5
            @warn "expecting 5 fields per access point (\"$str\")"
            continue
        end
        access = zero(GET)
        for c in arr[3]
            if c == 'g'
                access |= GET
            elseif c == 's'
                access |= SET
            elseif c == 'i'
                access |= INFO
            else
                @warn "unexpected access string (\"$(arr[3])\")"
                continue
            end
        end
        push!(lst, AccessPoint(arr[1], arr[2], arr[4], arr[5], access))
    end
    return lst
end

"""
    XPA.find([conn=XPA.connection(),] ident) -> apt

yields the accesspoint of the first XPA server matching `ident` or `nothing` if
none is found.  If a match is found, the result `apt` is an instance of
`XPA.AccessPoint` and has the following members:

    apt.class   # class of the access point (String)
    apt.name    # name of the access point
    apt.addr    # socket access method (host:port for inet,
    apt.user    # user name of access point owner
    apt.access  # allowed access

all members are `String`s but the last one, `access`, which is an `UInt`.

Argument `ident` may be a regular expression or a string of the
form `CLASS:NAME` where `CLASS` and `CLASS` are matched against the server
class and name respectively (they may be `"*"` to match any).

Optional argument `conn` is a persistent XPA client connection (created by
[`XPA.Client`](@ref)); if omitted, a per-thread connection is used (see
[`XPA.connection`](@ref)).

Keyword `user` may be used to specify the user name of the owner of the server
process, for instance `ENV["user"]` to match your servers.  The default
is `user="*"` which matches any user.

Keyword `throwerrors` may be set true (it is false by default) to automatically
throw an exception if no match is found (instead of returning `nothing`).

See also [`XPA.Client`](@ref), [`XPA.address`](@ref) and [`XPA.list`](@ref).

"""
find(ident::Union{AbstractString,Regex}; kwds...) =
    find(connection(), ident; kwds...)

function find(conn::Client,
              ident::AbstractString;
              user::AbstractString = "*",
              throwerrors::Bool = false)::Union{AccessPoint,Nothing}
    i = findfirst(isequal(':'), ident)
    if i === nothing
        # allow any class
        class = "*"
        name = ident
    else
        class = ident[1:i-1]
        name = ident[i+1:end]
    end
    anyuser = (user == "*")
    anyclass = (class == "*")
    anyname = (name == "*")
    lst = list(conn)
    for j in eachindex(lst)
        if ((anyuser || lst[j].user == user) &&
            (anyclass || lst[j].class == class) &&
            (anyname || lst[j].name == name))
            return lst[j]
        end
    end
    throwerrors && throw_no_servers_match(ident)
    return nothing
end

function find(conn::Client,
              ident::Regex;
              user::AbstractString = "*",
              throwerrors::Bool = false)::Union{AccessPoint,Nothing}
    anyuser = (user == "*")
    lst = list(conn)
    for j in eachindex(lst)
        if ((anyuser || lst[j].user == user) &&
            occursin(ident, lst[j].class*":"*lst[j].name))
            return lst[j]
        end
    end
    throwerrors && throw_no_servers_match(ident)
    return nothing
end

@noinline throw_no_servers_match(ident::AbstractString) =
    error("no XPA servers match pattern \"$(ident)\"")

@noinline throw_no_servers_match(ident::Regex) =
     error("no XPA servers match regular expression \"$(ident.pattern)\"")

"""
    XPA.address(apt) -> addr

yields the address of XPA accesspoint `apt` which can be: an instance of
`XPA.AccessPoint`, a string with a valid XPA server address or a server
`class:name` identifier.  In the latter case, [`XPA.find`](@ref) is called to
find a matching server which is much longer.

"""
address(apt::XPA.AccessPoint) =
    apt.addr

function address(apt::AbstractString)
    i = findfirst(isequal(':'), apt)
    if (i === nothing ||
        tryparse(UInt, apt[1:i-1],  base = 16) === nothing ||
        tryparse(UInt, apt[i+1:end],  base = 10) === nothing)
        return address(XPA.find(apt; throwerrors = true))
    end
    return apt
end

"""
    XPA.get([T, [dims,]] [conn,] apt, args...; kwds...)

retrieves data from one or more XPA access points identified by `apt` (a
template name, a `host:port` string or the name of a Unix socket file) with
arguments `args...` (automatically converted into a single string where the
arguments are separated by a single space).  Optional argument `conn` is a
persistent XPA client connection (created by [`XPA.Client`](@ref)); if omitted,
a per-thread connection is used (see [`XPA.connection`](@ref)).  The returned
value depends on the optional arguments `T` and `dims`.

If neither `T` nor `dims` are specified, an instance of [`XPA.Reply`](@ref) is
returned with all the answer(s) from the XPA server(s).  The following keywords
are available:

* Keyword `nmax` specifies the maximum number of answers, `nmax=1` by default.
  Specify `nmax=-1` to use the maximum number of XPA hosts.

* Keyword `throwerrors` specifies whether to check for errors.  If this keyword
  is set true, an exception is thrown for the first error message encountered
  in the list of answers.  By default, `throwerrors` is false.

* Keyword `mode` specifies options in the form `"key1=value1,key2=value2"`.

* Keyword `users` specifies the list of possible users owning the access-point.
  This (temporarily) overrides the settings in environment variable
  `XPA_NSUSERS`.  By default and if the environment variable `XPA_NSUSERS` is
  not set, the access-point must be owned the caller (see Section
  *Distinguishing Users* in XPA documentation).  The value is a string wich may
  be a list of comma separated user names or `"*"` to access all users on a
  given machine.

If `T` and, possibly, `dims` are specified, a single answer and no errors are
expected (as if `nmax=1` and `throwerrors=true`) and the data part of the
answer is converted according to `T` which must be a type and `dims` which is
an optional list of dimensions:

* If only `T` is specified, it can be `String` to return a string interpreting
  the data as ASCII characters or a type like `Vector{S}` to return the largest
  vector of elements of type `S` that can be extracted from the returned data.

* If both `T` and `dims` are specified, `T` can be a type like `Array{S}` or
  `Array{S,N}` and `dims` a list of `N` dimensions to retrieve the data as an
  array of type `Array{S,N}`.

See also [`XPA.Client`](@ref), [`XPA.get_data`](@ref), [`XPA.set`](@ref) and
[`XPA.verify`](@ref).

"""
function get(conn::Client,
             apt::AbstractString,
             cmd::AbstractString;
             mode::AbstractString = "",
             nmax::Integer = 1,
             throwerrors::Bool = false,
             users::Union{Nothing,AbstractString} = nothing)
    return _get(conn, apt, cmd, mode, _nmax(nmax), throwerrors, users)
end

get(apt::AccessPoint, args...; kwds...) =
    get(address(apt), args...; kwds...)

get(apt::AbstractString, args...; kwds...) =
    get(connection(), apt, args...; kwds...)

get(conn::Client, apt::AccessPoint, args...; kwds...) =
    get(conn, address(apt), args...; kwds...)

get(conn::Client, apt::AbstractString, args...; kwds...) =
    get(conn, apt, join_arguments(args); kwds...)

_get1(args...; kwds...) = get(args...; nmax = 1, throwerrors = true, kwds...)

function get(::Type{Vector{T}},
             args...; kwds...) :: Vector{T} where {T}
    get_data(Vector{T}, _get1(args...; kwds...))
end

function get(::Type{Vector{T}}, dim::Integer,
             args...; kwds...) :: Vector{T} where {T}
    get_data(Vector{T}, dim, _get1(args...; kwds...))
end

function get(::Type{Array{T}}, dim::Integer,
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Vector{T}, dim, _get1(args...; kwds...))
end

function get(::Type{Array{T}}, dims::NTuple{N,Integer},
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Array{T,N}, dims, _get1(args...; kwds...))
end

function get(::Type{Array{T,N}}, dims::NTuple{N,Integer},
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Array{T,N}, dims, _get1(args...; kwds...))
end

function get(::Type{String}, args...; kwds...)
    return get_data(String, _get1(args...; kwds...))
end

function _get(conn::Client, apt::AbstractString, params::AbstractString,
              mode::AbstractString, nmax::Int, throwerrors::Bool,
              users::Union{Nothing,AbstractString})
    lengths = fill!(Vector{Csize_t}(undef, nmax), 0)
    buffers = fill!(Vector{Ptr{Byte}}(undef, nmax*3), Ptr{Byte}(0))
    address = pointer(buffers)
    offset = nmax*sizeof(Ptr{Byte})
    prevusers = _override_nsusers(users)
    replies = ccall((:XPAGet, libxpa), Cint,
                    (Client, Cstring, Cstring, Cstring, Ptr{Ptr{Byte}},
                     Ptr{Csize_t}, Ptr{Ptr{Byte}}, Ptr{Ptr{Byte}}, Cint),
                    conn, apt, params, mode, address, lengths,
                    address + offset, address + 2*offset, nmax)
    _restore_nsusers(prevusers)
    0 ≤ replies ≤ nmax || error("unexpected number of replies from XPAGet")
    rep = finalizer(_free, Reply(replies, lengths, buffers))
    throwerrors && verify(rep; throwerrors=true)
    return rep
end

# Override environment variable XPA_NSUSERS.
_override_nsusers(::Nothing) = nothing
_override_nsusers(users::AbstractString) = begin
    prev = Base.get(ENV, "XPA_NSUSERS", "")
    ENV["XPA_NSUSERS"] = users
    return prev
end

# Restore environment variable XPA_NSUSERS.
_restore_nsusers(::Nothing) = nothing
_restore_nsusers(users::AbstractString) = begin
    if users == ""
        delete!(ENV, "XPA_NSUSERS")
    else
        ENV["XPA_NSUSERS"] = users
    end
    nothing
end

function _free(rep::Reply)
    nmax = _nmax(rep)
    fill!(rep.lengths, 0)
    for i in 0:2,
        j in 1:length(rep)
        k = i*nmax + j
        if (ptr = rep.buffers[k]) != NULL
            rep.buffers[k] = NULL
            _free(ptr)
        end
    end
end

# i-th server name is at index i + nmax.
_get_srv(rep::Reply, i::Integer) = _get_srv(rep, Int(i))
_get_srv(rep::Reply, i::Int) :: Ptr{Byte} =
    (1 ≤ i ≤ length(rep) ? rep.buffers[i + _nmax(rep)] : NULL)

# i-th message is at index i + 2*nmax.
_get_msg(rep::Reply, i::Integer) = _get_msg(rep, Int(i))
_get_msg(rep::Reply, i::Int) :: Ptr{Byte} =
    (1 ≤ i ≤ length(rep) ? rep.buffers[i + _nmax(rep)*2] : NULL)

# i-th data buffer is at index i.
_get_buf(rep::Reply, i::Integer, preserve::Bool) =
    _get_buf(rep, Int(i), preserve)

function _get_buf(rep::Reply, i::Int, preserve::Bool) :: Tuple{Ptr{Byte},Int}
    local ptr::Ptr{Byte}, len::Int
    if 1 ≤ i ≤ length(rep)
        ptr, len = rep.buffers[i], rep.lengths[i]
        (ptr == NULL ? len == 0 : len ≥ 0) || error("invalid buffer length")
        if ! preserve && ptr != NULL
            rep.lengths[i] = 0
            rep.buffers[i] = NULL
        end
    else
        ptr, len = NULL, 0
    end
    return (ptr, len)
end

"""
    XPA.join_arguments(args)

joins a tuple of arguments into a single string where arguments are separated
by a single space.  It is implemented so as to be faster than `join(args, " ")`
when `args` has less than 2 arguments.  It is intended to build XPA command
string from arguments.

"""
join_arguments(args::Tuple) = join(args, " ")
join_arguments(args::Tuple{AbstractString}) = args[1]
join_arguments(args::Tuple{Any}) = string(args[1])
join_arguments(::Tuple{}) = ""

"""

Private method `_nmax(n::Integer)` yields the maximum number of expected
answers to a get/set request.  The result is `n` if `n ≥ 1` or
`getconfig("XPA_MAXHOSTS")` otherwise.  The call `_nmax(rep::Reply)` yields
the maximum number of answer that can be stored in `rep`.

"""
_nmax(n::Integer) = (n == -1 ? Int(getconfig("XPA_MAXHOSTS")) : Int(n))
_nmax(rep::Reply) = length(rep.lengths)

"""
    XPA.get_server(rep, i=1)

yields the XPA identifier of the server which sent the `i`-th reply in XPA
answer `rep`.  An empty string is returned if there is no `i`-th reply.

See also [`XPA.get`](@ref), [`XPA.get_message`](@ref).

"""
get_server(rep::Reply, i::Integer=1) = _string(_get_srv(rep, i))

"""
    XPA.get_message(rep, i=1)

yields the message associated with the `i`-th reply in XPA answer `rep`.  An
empty string is returned if there is no `i`-th reply.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref), [`XPA.has_error`](@ref),
[`XPA.get_server`](@ref).

"""
get_message(rep::Reply, i::Integer=1) = _string(_get_msg(rep, i))

"""
    XPA.has_error(rep, i=1) -> boolean

yields whether `i`-th XPA answer `rep` contains an error message.  The error
message can be retrieved by calling `XPA.get_message(rep, i)`.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref),
[`XPA.get_message`](@ref).

"""
has_error(rep::Reply, i::Integer=1) =
    _startswith(_get_msg(rep, i), _XPA_ERROR)

const _XPA_ERROR_PREFIX = "XPA\$ERROR "
const _XPA_ERROR = Tuple(map(Byte, collect(_XPA_ERROR_PREFIX)))

function Base.show(io::IO, rep::Reply)
    print(io, "XPA.Reply")
    n = length(rep)
    if n == 0
        print(io, " (no replies)")
    else
        print(io, " (", n, " repl", (n > 1 ? "ies" : "y"), "):\n")
        for i in 1:n
            # Check whether all bytes in the data buffer are printable ASCII
            # characters.
            print(io, "  ", i, ": server = ",
                  repr(get_server(rep, i); context=io), ", message = ",
                  repr(get_message(rep, i); context=io), ", data = ")
            ptr, len = _get_buf(rep, i, true)
            if ptr == C_NULL
                print(io, "NULL")
            elseif len == 0
                print(io, repr("";  context=io))
            else
                cstring = true
                for j in 1:len
                    b = unsafe_load(ptr, j)
                    if (b & 0x80) != 0
                        # Not ASCII.
                        cstring = false
                        break
                    end
                    c = Char(b)
                    if !(isprint(c) || c == '\n' || c == '\r' || c == '\t')
                        # Not Printable.
                        cstring = false
                        break
                    end
                end
                if cstring
                    print(io, repr(unsafe_string(ptr, len);  context=io))
                else
                    print(io, len, (len > 1 ? " bytes" : " byte"))
                end
            end
            i < n && print(io, "\n")
        end
    end
end

"""
    XPA.has_errors(rep) -> boolean

yields whether answer `rep` contains any error messages.

See also [`XPA.get`](@ref), [`XPA.has_error`](@ref), [`XPA.get_message`](@ref).

"""
function has_errors(rep::Reply) :: Bool
    for i in 1:length(rep)
        if has_error(rep, i)
            return true
        end
    end
    return false
end

"""
    XPA.has_message(rep, i=1) -> boolean

yields whether `i`-th XPA answer `rep` contains an error message.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref).

"""
has_message(rep::Reply, i::Integer=1) =
    _startswith(_get_msg(rep, i), _XPA_MESSAGE)

const _XPA_MESSAGE_PREFIX = "XPA\$MESSAGE "
const _XPA_MESSAGE = Tuple(map(Byte, collect(_XPA_MESSAGE_PREFIX)))

function _startswith(ptr::Ptr{Byte}, tup::NTuple{N,Byte}) where {N}
    if ptr == NULL
        return false
    end
    for i in 1:N
        if unsafe_load(ptr, i) != tup[i]
            return false
        end
    end
    return true
end

"""
    XPA.verify(rep [, i]; throwerrors::Bool=false) -> boolean

verifies whether answer(s) in the result `rep` from an [`XPA.get`](@ref) or
[`XPA.set`](@ref) request has no errors.  If index `i` is specified only that
specific answer is considered; otherwise, all answers are verified.  If keyword
`throwerrors` is true, an exception is thrown for the first error found if any.

"""
function verify(rep::Reply; kwds...)
    for i in 1:length(rep)
        verify(rep, i; kwds...) || return false
    end
    return true
end

function verify(rep::Reply, i::Integer; throwerrors::Bool=false)
    if has_error(rep, i)
        if throwerrors
            # Strip error message prefix (which we know that it is present).
            msg = get_message(rep, i)
            j = length(_XPA_ERROR_PREFIX) + 1
            while j ≤ length(msg) && isspace(msg[j])
                j += 1
            end
            error(msg[j:end])
        else
            return false
        end
    end
    return true
end

"""
    XPA.get_data([T, [dims,]] rep, i=1; preserve=false)

yields the data associated with the `i`-th reply in XPA answer `rep`.  The
returned value depends on the optional leading arguments `T` and `dims`:

* If neither `T` nor `dims` are specified, a vector of bytes (`UInt8`) is
  returned.

* If only `T` is specified, it can be `String` to return a string interpreting
  the data as ASCII characters or a type like `Vector{S}` to return the largest
  vector of elements of type `S` that can be extracted from the data.

* If both `T` and `dims` are specified, `T` can be an array type like
  `Array{S}` or `Array{S,N}` and `dims` a list of `N` dimensions to retrieve
  the data as an array of type `Array{S,N}`.

Keyword `preserve` can be used to specifiy whether or not to preserve the
internal data buffer in `rep` for another call to `XPA.get_data`.  By default,
`preserve=true` when `T = String` is specified and `preserve=false` otherwise.

In any cases, the type of the result is predictible, so there should be no type
instability issue.

See also [`XPA.get`](@ref), [`XPA.get_message`](@ref),
[`XPA.get_server`](@ref).

"""
get_data(rep::Reply, args...; kwds...) =
    get_data(Vector{Byte}, rep, args...; kwds...)

function get_data(::Type{String}, rep::Reply, i::Integer=1;
                  preserve::Bool = true) :: String
    ptr, len = _get_buf(rep, i, preserve)
    ptr == NULL && return ""
    str = unsafe_string(ptr, len)
    preserve || _free(ptr)
    return str
end

function get_data(::Type{Vector{T}}, rep::Reply, i::Integer=1;
                  preserve::Bool = false) :: Vector{T} where {T}
    isbitstype(T) || error("invalid Array element type")
    ptr, len = _get_buf(rep, i, preserve)
    cnt = div(len, sizeof(T))
    if ptr == NULL || cnt ≤ 0
        # Empty vector.
        return T[]
    elseif preserve
        # Make a copy of the buffer in the form of a Julia vector.
        return _memcpy!(Vector{T}(undef, cnt), ptr, cnt*sizeof(T))
    else
        # Transfer ownership of buffer to Julia.
        return unsafe_wrap(Array, Ptr{T}(ptr), cnt, own=true)
    end
end

function get_data(::Type{Array{T}}, dim::Integer,
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Vector{T} where {T}
    return _get_buf(Vector{T}, _dimensions(dim), rep, i, preserve)
end

function get_data(::Type{Vector{T}}, dim::Integer,
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Vector{T} where {T}
    return _get_buf(Vector{T}, _dimensions(dim), rep, i, preserve)
end

function get_data(::Type{Array{T}}, dims::NTuple{N,Integer},
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Array{T,N} where {T,N}
    return _get_buf(Array{T,N}, _dimensions(dims), rep, i, preserve)
end

function get_data(::Type{Array{T,N}}, dims::NTuple{N,Integer},
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Array{T,N} where {T,N}
    return _get_buf(Array{T,N}, _dimensions(dims), rep, i, preserve)
end

"""

Private method `_get_buf(rep,i,preserve)` yields `(ptr,len)` the address and
length (in bytes) of internal buffer corresponding to the data for the `i`-th
reply in `rep`.  If `preserve` is false, then the internal buffer is set to
NULL and the caller is responsible to free it.  If `i` is out of range or if
there are no data associated with the `i`-th reply in `rep`, `(NULL,0)` is
returned.

The call:

```julia
_get_buf(::Type{Array{T,N}}, dims::NTuple{N,Int},
         rep::Reply, i::Integer, preserve::Bool)
```

yields the contents of the internal data buffer as a Julia array.

"""
function _get_buf(::Type{Array{T,N}}, dims::NTuple{N,Int},
                  rep::Reply, i::Int, preserve::Bool) :: Array{T,N} where {T,N}
    isbitstype(T) || error("invalid Array element type")
    minimum(dims) ≥ 0 || error("invalid Array dimensions")
    ptr, len = _get_buf(rep, i, preserve)
    cnt = prod(dims)
    cnt*sizeof(T) ≤ len || error("Array size too large for buffer")
    if ptr == NULL || cnt ≤ 0
        # Empty array.
        return Array{T,N}(undef, dims)
    elseif preserve
        # Make a copy of the buffer in the form of a Julia array.
        return _memcpy!(Array{T,N}(undef, dims), ptr, cnt*sizeof(T))
    else
        # Transfer ownership of buffer to Julia.
        return unsafe_wrap(Array, Ptr{T}(ptr), dims, own=true)
    end
end

_dimensions(dim::Integer) = (Int(dim),)
_dimensions(dim::Int) = (dim,)
_dimensions(dims::Tuple{}) = dims
_dimensions(dims::TupleOf{Integer}) = map(Int, dims)
_dimensions(dims::TupleOf{Int}) = dims

# _string(ptr) copies data at address ptr as a string, if ptr is NULL, an empty
# string is returned.
_string(ptr::Ptr{Byte}) = (ptr == NULL ? "" : unsafe_string(ptr))

"""
    XPA.set([conn,] apt, args...; data=nothing, kwds...) -> rep

sends `data` to one or more XPA access points identified by `apt` with
arguments `args...` (automatically converted into a single string where the
arguments are separated by a single space).  The result is an instance of
[`XPA.Reply`](@ref).  Optional argument `conn` is a persistent XPA client
connection (created by [`XPA.Client`](@ref)); if omitted, a per-thread
connection is used (see [`XPA.connection`](@ref)).

The following keywords are available:

* Keyword `data` specifies the data to send, may be `nothing`, an array or a
  string.  If it is an array, it must have contiguous elements (as a for a
  *dense* array) and must implement the `pointer` method.

* Keyword `nmax` specifies the maximum number of recipients, `nmax=1` by
  default.  Specify `nmax=-1` to use the maximum possible number of XPA hosts.

* Keyword `mode` specifies options in the form `"key1=value1,key2=value2"`.

* Keyword `throwerrors` specifies whether to check for errors.  If this keyword
  is set `true`, an exception is thrown for the first error message encountered
  in the list of answers.  By default, `throwerrors` is false.

* Keyword `users` specifies the list of possible users owning the access-point.
  This (temporarily) overrides the settings in environment variable
  `XPA_NSUSERS`.  By default and if the environment variable `XPA_NSUSERS` is
  not set, the access-point must be owned the caller (see Section
  *Distinguishing Users* in XPA documentation).  The value is a string wich may
  be a list of comma separated user names or `"*"` to access all users on a
  given machine.

See also [`XPA.Client`](@ref), [`XPA.get`](@ref) and [`XPA.verify`](@ref).

"""
function set(conn::Client,
             apt::AbstractString,
             cmd::AbstractString;
             data = nothing,
             mode::AbstractString = "",
             nmax::Integer = 1,
             throwerrors::Bool = false,
             users::Union{Nothing,AbstractString} = nothing)
    return _set(conn, apt, cmd, mode, buffer(data), _nmax(nmax),
                throwerrors, users)
end

function set(conn::Client,
             apt::AbstractString,
             args::Union{AbstractString,Real}...;
             kwds...)
    return _set(conn, apt, join_arguments(args); kwds...)
end

set(apt::AbstractString, args::Union{AbstractString,Real}...; kwds...) =
    set(connection(), apt, join_arguments(args); kwds...)

function _set(conn::Client, apt::AbstractString, params::AbstractString,
              mode::AbstractString, data::Union{NullBuffer,DenseArray},
              nmax::Int, throwerrors::Bool,
              users::Union{Nothing,AbstractString})
    lengths = fill!(Vector{Csize_t}(undef, nmax), 0)
    buffers = fill!(Vector{Ptr{Byte}}(undef, nmax*3), Ptr{Byte}(0))
    address = pointer(buffers)
    offset = nmax*sizeof(Ptr{Byte})
    prevusers = _override_nsusers(users)
    replies = ccall((:XPASet, libxpa), Cint,
              (Client, Cstring, Cstring, Cstring, Ptr{Cvoid},
               Csize_t, Ptr{Ptr{Byte}}, Ptr{Ptr{Byte}}, Cint),
              conn, apt, params, mode, data, sizeof(data),
                    address + offset, address + 2*offset, nmax)
    _restore_nsusers(prevusers)
    0 ≤ replies ≤ nmax || error("unexpected number of replies from XPASet")
    rep = finalizer(_free, Reply(replies, lengths, buffers))
    throwerrors && verify(rep; throwerrors=true)
    return rep
end

"""
    buf = XPA.buffer(data)

yields an object `buf` representing the contents of `data` and which can be
used as an argument to `ccall` without the risk of having the data garbage
collected.  Argument `data` can be `nothing`, a dense array or a string.  If
`data` is an array `buf` is just an alias for `data`.  If `data` is a string,
`buf` is a temporary byte buffer where the string has been copied.

Standard methods `pointer` and `sizeof` can be applied to `buf` to retieve the
address and the size (in bytes) of the data and `convert(Ptr{Cvoid},buf)` can
also be used.

See also [`XPA.set`](@ref).

"""
function buffer(arr::A) :: A where {T,N,A<:DenseArray{T,N}}
    @assert isbitstype(T)
    return arr
end

function buffer(str::AbstractString)
    @assert isascii(str)
    len = length(str)
    buf = Vector{Cchar}(undef, len)
    @inbounds for i in 1:len
        buf[i] = str[i]
    end
    return buf
end

buffer(::Nothing) = NullBuffer()

Base.unsafe_convert(::Type{Ptr{T}}, ::NullBuffer) where {T} = Ptr{T}(0)
Base.pointer(::NullBuffer) = C_NULL
Base.sizeof(::NullBuffer) = 0


text:
using Test
using Random
using Distributions
using MEstimation
using NLsolve

include("new_logistic_regression.jl")

Random.seed!(123);
n = 50;
p = 10;
x = rand(n, p);
x[:, 1] = fill(1.0, n);
true_betas = rand(p) * 0.02;
my_data = new_logistic_regression.simulate(true_betas, x, fill(1, n));

logistic_template = MEstimation.objective_function_template(
    new_logistic_regression.nobs,
    new_logistic_regression.observation,
    new_logistic_regression.loglik
)

logistic_obj_template = MEstimation.objective_function_template(new_logistic_regression.nobs,
                                                    new_logistic_regression.observation,
                                                    new_logistic_regression.loglik)

obj_br_implicit = MEstimation.fit(logistic_obj_template, my_data, true_betas, estimation_method="RBM")

obj_br_explicit = MEstimation.fit(logistic_obj_template, my_data, true_betas, estimation_method="RBM", br_method="explicit_trace")

@test isapprox(coef(obj_br_implicit), coef(obj_br_explicit))


text:
@recipe function f(ls::AbstractVector{<:AbstractLinePrimitive})
    linecolor --> :black

    @series begin
        label --> "Edges"
        ls[1]
    end
    if length(ls) > 1
        for l in ls[2:end]
            @series begin
                label := nothing
                l
            end
        end
    end
end

include("Edge.jl")
include("Ellipse.jl")


text:
# Autogenerated wrapper script for Libassuan_jll for x86_64-apple-darwin
export libassuan

using Libgpg_error_jll
JLLWrappers.@generate_wrapper_header("Libassuan")
JLLWrappers.@declare_library_product(libassuan, "@rpath/libassuan.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Libgpg_error_jll)
    JLLWrappers.@init_library_product(
        libassuan,
        "lib/libassuan.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Test of binary heaps

@testset "BinaryHeaps" begin

    @testset "make heap" begin
        vs = [4, 1, 3, 2, 16, 9, 10, 14, 8, 7]

        @testset "make min heap" begin
            h = BinaryMinHeap(vs)

            @test length(h) == 10
            @test !isempty(h)
            @test first(h) == 1
            @test isequal(h.valtree, [1, 2, 3, 4, 7, 9, 10, 14, 8, 16])
            @test sizehint!(h, 100) === h
        end

        @testset "make max heap" begin
            h = BinaryMaxHeap(vs)

            @test length(h) == 10
            @test !isempty(h)
            @test first(h) == 16
            @test isequal(h.valtree, [16, 14, 10, 8, 7, 3, 9, 1, 4, 2])
            @test sizehint!(h, 100) === h
        end

        @testset "push!" begin
            @testset "push! hmin" begin
                hmin = BinaryMinHeap{Int}()
                @test length(hmin) == 0
                @test isempty(hmin)

                ss = Any[
                    [4],
                    [1, 4],
                    [1, 4, 3],
                    [1, 2, 3, 4],
                    [1, 2, 3, 4, 16],
                    [1, 2, 3, 4, 16, 9],
                    [1, 2, 3, 4, 16, 9, 10],
                    [1, 2, 3, 4, 16, 9, 10, 14],
                    [1, 2, 3, 4, 16, 9, 10, 14, 8],
                    [1, 2, 3, 4, 7, 9, 10, 14, 8, 16]]

                for i = 1 : length(vs)
                    push!(hmin, vs[i])
                    @test length(hmin) == i
                    @test !isempty(hmin)
                    @test isequal(hmin.valtree, ss[i])
                end

                @testset "pop! hmin" begin
                    @test isequal(extract_all!(hmin), [1, 2, 3, 4, 7, 8, 9, 10, 14, 16])
                    @test isempty(hmin)
                end

            end

            @testset "push! hmax" begin
                hmax = BinaryMaxHeap{Int}()
                @test length(hmax) == 0
                @test isempty(hmax)

                ss = Any[
                    [4],
                    [4, 1],
                    [4, 1, 3],
                    [4, 2, 3, 1],
                    [16, 4, 3, 1, 2],
                    [16, 4, 9, 1, 2, 3],
                    [16, 4, 10, 1, 2, 3, 9],
                    [16, 14, 10, 4, 2, 3, 9, 1],
                    [16, 14, 10, 8, 2, 3, 9, 1, 4],
                    [16, 14, 10, 8, 7, 3, 9, 1, 4, 2]]

                for i = 1 : length(vs)
                    push!(hmax, vs[i])
                    @test length(hmax) == i
                    @test !isempty(hmax)
                    @test isequal(hmax.valtree, ss[i])
                end

                @testset "pop! hmax" begin
                    @test isequal(extract_all!(hmax), [16, 14, 10, 9, 8, 7, 4, 3, 2, 1])
                    @test isempty(hmax)
                end
            end
        end

    end

    @testset "hybrid push! and pop!" begin
        h = BinaryMinHeap{Int}()

        @testset "push1" begin
            push!(h, 5)
            push!(h, 10)
            @test isequal(h.valtree, [5, 10])
        end

        @testset "pop1" begin
            @test pop!(h) == 5
            @test isequal(h.valtree, [10])
        end

        @testset "push2" begin
            push!(h, 7)
            push!(h, 2)
            @test isequal(h.valtree, [2, 10, 7])
        end

        @testset "pop2" begin
            @test pop!(h) == 2
            @test isequal(h.valtree, [7, 10])
        end
    end

    @testset "nlargest and nsmallest" begin
        ss = [100,103,-12,-109,67,4,65,-52,-97,-32,-24,114,-128,
            102,-56,-17,-41,25,-30,-84,26,-84,48,49,-5,-38,28,
            114,-54,96,-55,67,74,127,-61,124,11,-7,93,-51,110,
            -106,-84,-90,-18,-12,-116,21,115,50]
        for n = -1:length(ss) + 1
            @test sort(ss, lt = >)[1:min(n, end)] == nlargest(n, ss)
            @test sort(ss, lt = <)[1:min(n, end)] == nsmallest(n, ss)
        end
    end

    @testset "push! type conversion" begin # issue 399
        h = BinaryMinHeap{Float64}()
        push!(h, 3.0)
        push!(h, 5)
        push!(h, Rational(4, 8))
        push!(h, Complex(10.1, 0.0))

        @test isequal(h.valtree, [0.5, 5.0, 3.0, 10.1])
    end
end # @testset BinaryHeap


text:
"""
    get_version(m::Module)

Get the version of a given module. It will try to find the version
of `Project.toml` if the given module is a project module. If fails,
it returns `v"0.0.0"`.
"""
function get_version(m::Module)
    # no version for scripts
    m === Main && return
    # project module
    path = pathof(m)
    if path !== nothing
        envpath = dirname(dirname(path))
        project = Pkg.Types.read_project(Pkg.Types.projectfile_path(envpath))
        if project.name == string(nameof(m)) && project.version !== nothing
            return project.version
        end
    end

    if isdefined(m, :COMMAND_VERSION)
        return m.COMMAND_VERSION
    end
    return
end

"""
    set_cmd!(cmds::Dict, cmd)

register `cmd` in the command registry `cmds`, which usually is
a constant `CASTED_COMMANDS` under given module.
"""
function set_cmd!(cmds::Dict, cmd, name = cmd.name)
    if haskey(cmds, name)
        @warn "replacing command $name in the registry"
    end

    return cmds[name] = cmd
end

"""
    default_name(x::String)

Return the lowercase of a given package name. It will
ignore the suffix if it ends with ".jl".
"""
function default_name(x::String)
    if endswith(x, ".jl")
        name = x[1:end-3]
    else
        name = x
    end
    return replace(lowercase(name), '_' => '-')
end

default_name(x::Symbol) = default_name(string(x))

function write_cmd(filename, cmd::Entry)
    write(filename, string(emit_expr(cmd)))
end


text:
m = @model x0 begin
    μ ~ HalfNormal(1)
    σ ~ HalfCauchy(1)
    x ~ Normal(x0 + μ, σ)
end

m_fwd = m(:μ,:σ)
m_inv = m(:x)

r = makeRand(m_fwd)

grt = (μ = 0.1, σ = 1.2735923758662662, x0=0.0)

function realize(n; pars...)
    pars = Dict(pars)
    xs = zeros(n)
    for j in 1:n
        newx = r(;pars...).x
        xs[j] = newx
        pars = merge(pars, Dict(:x0=>newx))

    end
    xs
end

obs = realize(1000;grt...)

using Plots
pyplot()
plot(obs)


text:
# Generate plots based on data from dev/demo_unsteady.jl

# Speed up loading
module ModifiedStokesSolver
include("../src/types.jl")
end

using FileIO
using HDF5
using PyPlot
using Interpolations

include("../src/plottools.jl")
include("../src/streakplot.jl")

function load_geo(datafile)
    info("Loading ", datafile)
    D = load(datafile)
    X = D["X"]
    Y = D["Y"]
    dcurve = D["dcurve"]
    interior = D["interior"]
    return X, Y, dcurve, interior
end

function make_vortex_plot(;Re=100, vorticity=true, stepnum=15000,
                          streaklines=[], streamlines=false,
                          plotcolorbar=false)
    basedir = "data/vortex_Re$(Re)"
    geoname = "$basedir/geo.jld"
    dataname = i -> @sprintf("%s/vstep_%06d.jld", basedir, i)
    datafile = dataname(stepnum)
    X, Y, dcurve, interior = load_geo(geoname)
    vorticity, U1, U2 = [], [], []
    t = 0.0;
    h5open(datafile, "r") do io
        U1 = read(io, "U1")
        U2 = read(io, "U2")
        dU2dx = read(io, "dU2dx")
        dU1dy = read(io, "dU1dy")
        t = read(io, "ti")
        vorticity = dU2dx - dU1dy
    end
    @. vorticity[interior==false] = NaN

    h5open(datafile, "r") do io
        dU2dx = read(io, "dU2dx")
        dU1dy = read(io, "dU1dy")
        t = read(io, "ti")
        vorticity = dU2dx - dU1dy
    end
    @. vorticity[interior==false] = NaN
    clf()
    s = 1
    plot_boundary(dcurve)
    zoom_to_boundary(dcurve)
    axis("off")
    pcolormesh(X, Y, vorticity, cmap="coolwarm", shading="gouraud")
    clim(-3, 3)
    if plotcolorbar
        colorbar(label="Vorticity", shrink=0.8, pad=0.01)
    end
    if streamlines
        # Start points at inflow and in wake
        ystart = 0.1:0.5:3
        ystart = vcat(-reverse(ystart), ystart)
        xstart = zeros(length(ystart)) .- 4.9
        push!(xstart, -2.5)
        push!(ystart, 0.03)
        push!(xstart, -2.5)
        push!(ystart, -0.03)

        streamplot(y=Y[1,:],
                   x=X[:,1],
                   u=U1',
                   v=U2',
                   density=40,
                   start_points=hcat( xstart, ystart-0.1),
                   #integration_direction="both",
                   maxlength=4,
                   color="black"

                   )
    end

    if !isempty(streaklines)
        starty = 0.1:0.5:3
        starty = vcat(reverse(-starty), starty) .- 0.1
        startx = ones(size(starty)) .* (-5)
        streakplot(X=X, Y=Y, interior=interior,
                   startx=startx,
                   starty=starty,
                   dataname=dataname,
                   ilist=streaklines)
    end

    zoom_to_boundary(dcurve)
    tight_layout(0)

    #pngname = "../docs/paper/figs/vortex_Re$(Re).png"
    #savefig(pngname, dpi=200)
    #run(`mogrify -trim $pngname`)
end

close("all")
figsize = (5.5, 3)

# figure(1, figsize=figsize)
# make_vortex_plot(Re=100, stepnum=6250)

# figure(2, figsize=figsize)
# make_vortex_plot(Re=200, stepnum=14500, plotcolorbar=true)

#figure(3, figsize=figsize)
#sn = 14600
#make_vortex_plot(Re=50, stepnum=sn, plotcolorbar=true, streaklines=sn-4000:1:sn)

figure(4, figsize=figsize)
make_vortex_plot(Re=25, stepnum=10, streamlines=true)


text:
using AOPE
using Test

@testset "AOPE.jl" begin
    # Write your tests here.
end


text:
abstract type RenderTask end

"""
    SamplerRendererTask{S<:Scene, C<:Camera, I<:SurfaceIntegrator}

A wrapper for the scene, camera and integrator as well as a task number and total count of
tasks. This information represents a piece of the full render job and is used to render
the correct pixels of the output.

"""
struct SamplerRendererTask{S<:Scene, C<:Camera, I<:SurfaceIntegrator, T<:Sampler} <: RenderTask
    scene::S
    camera::C
    integrator::I
    sampler::T
    number::Int64
    count::Int64
end

"""
    enqueue_and_run(Tasks::Array{S,1})

Add the given tasks to a queue and run them. Currently this is a dummy version which just
runs all the tasks one by one in a single thread.

The computation is multi-threaded on a single processor. The number of threads is one by
default, and can be changed by setting the `JULIA_NUM_THREADS` environmental variable. A
lock is used to prevent multiple data from accessing the camera film simultaneously. In
many films, a ray can affect multiple pixels or bins, including ones "belonging" to a
different thread.

"""
function enqueue_and_run(Tasks::Array{S,1}) where S<:RenderTask
    lock = Threads.SpinLock()
    Threads.@threads for task in Tasks
        run(task, lock)
    end
end

function enqueue_debug(Tasks::Array{S,1}) where S<:RenderTask
    lock = Threads.SpinLock()
    for (t, task) in enumerate(Tasks)
        run(task, lock)
    end
end

"""
    render(scene::Scene, camera::Camera, integrator::SurfaceIntegrator, sampler::Sampler)

Render a given scene, as seen by a given camera, using given SurfaceIntegrator and Sampler.
"""
function render(scene::Scene, camera::Camera, integrator::SurfaceIntegrator,
                sampler::Sampler ; debug=false)
    # initialize sample storage
    # i.e. make an array of Samples to be filled?

    # create and launch tasks for rendering
    nTasks = count_tasks(camera, nprocs())
    tasks = [SamplerRendererTask(scene, camera, integrator, sampler, n, nTasks) for n = 1:nTasks]

    if debug
        println("Running $(nTasks) render tasks in single-thread debug mode...")
        enqueue_debug(tasks)
    else
        println("Running $(nTasks) render tasks on $(Threads.nthreads()) threads...")
        enqueue_and_run(tasks)
    end
end

"""
    run(task::SamplerRendererTask, write_lock::Threads.AbstractLock)

Run a SamplerRendererTask. This is the core function of the raytracer: here we generate a
ray from the camera, trace its possible intersection with the scene, compute reflection of
all light sources at that point and add up the intensity to the camera film.

The thread lock given by the second argument, `write_lock` is used to lock the camera film
so it can only be accessed by one thread at a time. If a significant amount of time is
spent in the `add_sample!` function, multi-threading could lead to a performance drop. This
depends on the film.

"""
function run(task::SamplerRendererTask, write_lock::Threads.AbstractLock)
    subsampler = get_subsampler(task.sampler, task.number, task.count)
    subsampler == nothing && return nothing

    max_samples = 0 # TODO
    state = 0

    Nsamples = subsampler.xs * subsampler.ys
    samples = Array{CameraSample}(undef, Nsamples)

    # TODO: This would more elegantly be an iterator, by defining iterate(subsampler, state)
    # and iterate(subsampler). But in that case the recycling of the `samples` array may
    # not work, which would lead to a performance drop. Experimentation needed.
    while true
        finished(subsampler, state) && break
        state = get_samples!(subsampler, state, samples)

        # generate camera rays
        for i = 1:length(samples)
            # find camera ray for sample[i]
            weight, ray = generate_ray(task.camera, samples[i])
            # evaluate radiance along ray and add it to camera film
            if weight > 0.0
                maybe_isect = intersect(ray, task.scene)
                if maybe_isect != nothing
                    Ls = weight .* intensity(task.integrator, task.scene, maybe_isect,
                                            ray, samples[i])
                else
                    Ls = task.scene.background
                end
                lock(write_lock)
                add_sample!(task.camera.film, samples[i], Ls, maybe_isect)
                unlock(write_lock)

            end
        end
    end
    return nothing
end

"""
    PressureRendererTask{S<:Scene, T<:Sampler} <: RenderTask

A structure that contains the information needed to render one part of the scene.

"""
struct PressureRendererTask{S<:Scene, T<:Sampler} <: RenderTask
    scene::S
    integrator::PressureIntegrator
    sampler::T
    nlight::Int64
    number::Int64
    count::Int64
    sphere::BoundingSphere
end

"""
    function render(scene::Scene, integrator::PressureIntegrator, sampler::Sampler)

Compute the radiation pressure and torque on a given scene, using given PressureIntegrator
and Sampler.

"""
function render(scene::Scene, integrator::PressureIntegrator, sampler::Sampler ; debug=false)
    # create and launch tasks for rendering

    N_lights = length(scene.lights)
    nbins = sampler.xend * sampler.yend
    rays_per_bin = sampler.xs * sampler.ys
    tasks_per_light = round_pow2(max(div(nbins * rays_per_bin, 1024), 4*nprocs()))
    N_tasks = N_lights * tasks_per_light
    tasks = PressureRendererTask[]
    sphere = BoundingSphere(scene.bounds)
    for (l, light) in enumerate(scene.lights)
        for n = 1:tasks_per_light
            push!(tasks, PressureRendererTask(scene, integrator, sampler, l,
                    n, tasks_per_light, sphere))
        end
        println()
    end

    if debug
        println("Running $(N_tasks) render tasks in single-thread debug mode...")
        enqueue_debug(tasks)
    else
        println("Running $(N_tasks) render tasks on $(Threads.nthreads()) threads...")
        enqueue_and_run(tasks)
    end

end

function run(task::PressureRendererTask, write_lock::Threads.AbstractLock)
    subsampler = get_subsampler(task.sampler, task.number, task.count)
    subsampler == nothing && return nothing
    Nsamples = subsampler.xs * subsampler.ys
    samples = Array{CameraSample}(undef, Nsamples)

    light = task.scene.lights[task.nlight]

    state = 0

    while !finished(subsampler, state)
        state = get_samples!(subsampler, state, samples)

        force = Vector3(0)
        torque = Vector3(0)

        for i = 1:length(samples)
            sample = normalize(samples[i], subsampler.xnorm, subsampler.ynorm)
            ray = generate_ray(light, task.sphere, sample)
            isect = intersect(ray, task.scene)
            if isect != nothing
                f, t = compute_pressure(task.integrator, light, isect, ray)
                force += f
                torque += t
            end
        end
        lock(write_lock)
        task.integrator.force[task.nlight] += force * cross_section(task.sphere)
        task.integrator.torque[task.nlight] += torque * cross_section(task.sphere)
        task.integrator.counts[task.nlight] += length(samples)
        unlock(write_lock)
    end
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using MKL

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

#import LinearAlgebra, OpenBLAS32_jll
#LinearAlgebra.BLAS.lbt_forward(OpenBLAS32_jll.libopenblas_path)

using PatternFormation
using LinearSolve
using Plots
using OrdinaryDiffEq, LinearAlgebra
using FLoops
using BenchmarkTools
using Distributions

# Grid and initial conditions
const N = 900
myEquation = NTF_Periodic!
tspan = (0.0, 30000.0)

Mᵤ = 1e-5
Mᵥ = 1e-5
γᵤ = 1e-2
γᵥ = 1e-2
a = 0.5
c = 0.00125
N_threads = 1
BLAS.set_num_threads(N_threads)
ex = ThreadedEx(simd = true)

println("Hello")

#u01, x, y, dx, dy = init_cond(:ErMnO3_MnOnly, N, M = 21, α=3000, γ = 500, BC=:Neumann, r₀=0.01, n_vd = 100, n_id = 100)
#u02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 21, α=500, γ = 3000, BC=:Neumann, r₀=0.01, n_vd = 100, n_id = 100)
#u01, x, y, dx, dy = init_cond(:ErMnO3_MnOnly, N, M = 21, α=3000, BC=:Neumann, r₀=0.01)
#u02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 21, α=500, BC=:Neumann, r₀=0.01)
#u01, x, y, dx, dy = init_cond(:ErMnO3_MnOnly, N, M = 9, α=3000, BC=:Periodic)
#u02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 9, α=500, BC=:Periodic)
#u01 = 1 .- u01
u01, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = 0.01, dx = 1/3)
u02, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = 0.25, dx = 1/3)
#u01 = sin.(x .+ y')
#u02 = cos.(x .+ y')
u01 .= 0.01 .+ rand(Uniform(-0.01/100, +0.01/100), N, N)
u02 .= 0.25 .+ rand(Uniform(-0.25/100, +0.25/100), N, N)
#u01, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 13, α=3000, BC=:Neumann, c01 = 0.01)
#u02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 13, α=500, BC=:Neumann, c01 = 0.25)
u0 = cat(u01, u02, dims=3)
heatmap(x, y, u01')
heatmap(x, y, u02')
p = [Mᵤ, Mᵥ, γᵤ, γᵥ, a, c, dx, dy, N]

func(du,u,p,t) = myEquation(du, u, p, t, ex)
prob = ODEProblem(func, u0, tspan, p)

# Solve!
#println("Solving")
@time sol = solve(prob, ROCK2(), saveat=range(0, stop=tspan[2], length=101), progress=true, progress_steps=1)
#@time sol = solve(prob, Tsit5(), progress=true, progress_steps=1)
println("Done!")

# Plot!
println("Plotting")
anim = @animate for i in 1:length(sol.t)
  tit = "t = $(sol.t[i])"
  #u = kron(ones(3,3), sol.u[i][:,:,1])
  u = sol.u[i][:,:,2]
  #heatmap(u, c=:grays, aspect_ratio=dx/dy, axis=([], false), colorbar=false, clims=(minimum(u),maximum(u)), title=tit)
  heatmap(u, c=:grays, aspect_ratio=dx/dy, axis=([], false), colorbar=true, title=tit)
end
#
gif(anim, fps=10)