text:
@doc raw"""
    Hyperbolic{N} <: AbstractEmbeddedManifold{DefaultEmbeddingType}

The hyperbolic space $ℍ^n$ represented by $n+1$-Tuples, i.e. embedded in the
[`Lorentz`](@ref)ian manifold equipped with the [`MinkowskiMetric`](@ref)
$⟨\cdot,\cdot⟩_{\mathrm{M}}$. The space is defined as

```math
ℍ^n = \Bigl\{p ∈ ℝ^{n+1} : ⟨p,p⟩_{\mathrm{M}}= -p_{n+1}^2
  + \displaystyle\sum_{k=1}^n p_k^2 = -1, p_{n+1} > 0\Bigr\},.
```

The tangent space $T_p ℍ^n$ is given by

````math
T_p ℍ^n := \bigl\{
X ∈ ℝ^{n+1} : ⟨p,X⟩_{\mathrm{M}} = 0
\bigr\}.
````
Note that while the [`MinkowskiMetric`](@ref) renders the [`Lorentz`](@ref) manifold (only)
pseudo-Riemannian, on the tangent bundle of the Hyperbolic space it induces a Riemannian
metric. The corresponding sectional curvature is $-1$.

# Constructor

    Hyperbolic(n)

Generate the $ℍ^{n} ⊂ ℝ^{n+1}$
"""
struct Hyperbolic{N} <: AbstractEmbeddedManifold{DefaultIsometricEmbeddingType} end

Hyperbolic(n::Int) = Hyperbolic{n}()

"""
    check_manifold_point(M::Hyperbolic, p; kwargs...)

Check whether `p` is a valid point on the [`Hyperbolic`](@ref) `M`, i.e. is a `p` vector with
inner product in the embedding of -1, see [`MinkowskiMetric`](@ref).
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::Hyperbolic, p; kwargs...)
    mpv = invoke(check_manifold_point, Tuple{supertype(typeof(M)), typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    if !isapprox(minkowski_metric(p, p), -1.0; kwargs...)
        return DomainError(
            minkowski_metric(p, p),
            "The point $(p) does not lie on $(M) since its Minkowski inner product is not -1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::Hyperbolic, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to `p` on the [`Hyperbolic`](@ref) `M`, i.e.
after [`check_manifold_point`](@ref)`(M,p)`, `X` has to be of the same dimension as `p`
and orthogonal to `p` with respect to the inner product from the embedding, see
[`MinkowskiMetric`](@ref). The optional parameter `check_base_point` indicates whether to
call [`check_manifold_point`](@ref)  for `p`. The tolerance for the last test can be set
using the `kwargs...`.
"""
function check_tangent_vector(M::Hyperbolic, p, X; check_base_point = true, kwargs...)
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)), typeof(p), typeof(X)},
        M,
        p,
        X;
        check_base_point = false, # already checked above
        kwargs...
    )
    mpv === nothing || return mpv
    if !isapprox(minkowski_metric(p, X), 0.0; kwargs...)
        return DomainError(
            abs(minkowski_metric(p, X)),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal (with respect to the Minkowski inner product) in the embedding.",
        )
    end
    return nothing
end

decorated_manifold(M::Hyperbolic{N}) where {N} = Lorentz(N+1, MinkowskiMetric())

default_metric_dispatch(::Hyperbolic, ::MinkowskiMetric) = Val(true)

@doc raw"""
    distance(M::Hyperbolic, p, q)

Compute the distance on the [`Hyperbolic`](@ref) `M`, which reads

````math
d_{ℍ^n}(p,q) = \operatorname{acosh}( - ⟨p, q⟩_{\mathrm{M}}),
````

where $⟨\cdot,\cdot⟩_{\mathrm{M}}$ denotes the [`MinkowskiMetric`](@ref) on the embedding,
the [`Lorentz`](@ref)ian manifold.
"""
distance(M::Hyperbolic, p, q) = acosh(max(-minkowski_metric(p, q), 1.0))

embed!(::Hyperbolic, q, p) = (q .= p)

embed!(::Hyperbolic, Y, p, X) = (Y .= X)

@doc raw"""
    exp(M::Hyperbolic, p, X)

Compute the exponential map on the [`Hyperbolic`](@ref) space $ℍ^n$ emanating
from `p` towards `X`. The formula reads

````math
\exp_p X = \cosh(\sqrt{⟨X,X⟩_{\mathrm{M}}})p
+ \sinh(\sqrt{⟨X,X⟩_{\mathrm{M}}})\frac{X}{\sqrt{⟨X,X⟩_{\mathrm{M}}}},
````

where $⟨\cdot,\cdot⟩_{\mathrm{M}}$ denotes the [`MinkowskiMetric`](@ref) on the embedding,
the [`Lorentz`](@ref)ian manifold.
"""
exp(::Hyperbolic, ::Any...)

function exp!(M::Hyperbolic, q, p, X)
    vn = sqrt(max(inner(M, p, X, X), 0.0))
    vn < eps(eltype(p)) && return copyto!(q, p)
    return copyto!(q, cosh(vn) * p + sinh(vn) / vn * X)
end

@doc raw"""
    injectivity_radius(M::Hyperbolic)
    injectivity_radius(M::Hyperbolic, p)

Return the injectivity radius on the [`Hyperbolic`](@ref), which is $∞$.
"""
injectivity_radius(H::Hyperbolic) = Inf
injectivity_radius(H::Hyperbolic, ::ExponentialRetraction) = Inf
injectivity_radius(H::Hyperbolic, ::Any) = Inf
injectivity_radius(H::Hyperbolic, ::Any, ::ExponentialRetraction) = Inf
eval(quote
    @invoke_maker 1 Manifold injectivity_radius(M::Hyperbolic, rm::AbstractRetractionMethod)
end)

@doc raw"""
    log(M::Hyperbolic, p, q)

Compute the logarithmic map on the [`Hyperbolic`](@ref) space $ℍ^n$, the tangent
vector representing the [`geodesic`](@ref) starting from `p`
reaches `q` after time 1. The formula reads for $p ≠ q$

```math
\log_p q = d_{ℍ^n}(p,q)
\frac{q-⟨p,q⟩_{\mathrm{M}} p}{\lVert q-⟨p,q⟩_{\mathrm{M}} p \rVert_2},
```

where $⟨\cdot,\cdot⟩_{\mathrm{M}}$ denotes the [`MinkowskiMetric`](@ref) on the embedding,
the [`Lorentz`](@ref)ian manifold. For $p=q$ the logarihmic map is equal to the zero vector.
"""
log(::Hyperbolic, ::Any...)

function log!(M::Hyperbolic, X, p, q)
    scp = minkowski_metric(p, q)
    w = q + scp * p
    wn = sqrt(max(scp .^ 2 - 1, 0.0))
    wn < eps(eltype(p)) && return zero_tangent_vector!(M, X, p)
    X .= acosh(max(1.0, -scp)) / wn .* w
    return X
end

@doc raw"""
    manifold_dimension(H::Hyperbolic)

Return the dimension of the hyperbolic space manifold $ℍ^n$, i.e. $\dim(ℍ^n) = n$.
"""
manifold_dimension(::Hyperbolic{N}) where {N} = N

"""
    mean(
        M::Hyperbolic,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = CyclicProximalPointEstimationMethod();
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` on the
[`Hyperbolic`](@ref) space using [`CyclicProximalPointEstimation`](@ref).
"""
mean(::Hyperbolic, ::Any...)

function mean!(M::Hyperbolic, p, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(M, p, x, w, CyclicProximalPointEstimation(); kwargs...)
end

@doc raw"""
    project(M::Hyperbolic, p, X)

Perform an orthogonal projection with respect to the Minkowski inner product of `X` onto
the tangent space at `p` of the [`Hyperbolic`](@ref) space `M`.

The formula reads
````math
Y = X + ⟨p,X⟩_{\mathrm{M}} p,
````
where $⟨\cdot, \cdot⟩_{\mathrm{M}}$ denotes the [`MinkowskiMetric`](@ref) on the embedding,
the [`Lorentz`](@ref)ian manifold.
"""
project(::Hyperbolic, ::Any, ::Any)

project!(M::Hyperbolic, Y, p, X) = (Y .= X .+ minkowski_metric(p, X) .* p)

show(io::IO, ::Hyperbolic{N}) where {N} = print(io, "Hyperbolic($(N))")

@doc raw"""
    vector_transport_to(M::Hyperbolic, p, X, q, ::ParallelTransport)

Compute the paralllel transport of the `X` from the tangent space at `p` on the
[`Hyperbolic`](@ref) space $ℍ^n$ to the tangent at `q` along the [`geodesic`](@ref)
connecting `p` and `q`. The formula reads

````math
\mathcal P_{q←p}X = X - \frac{⟨\log_p q,X⟩_p}{d^2_{ℍ^n}(p,q)}
\bigl(\log_p q + \log_qp \bigr),
````
where $⟨\cdot,\cdot⟩_p$ denotes the inner product in the tangent space at `p`.
"""
vector_transport_to(::Hyperbolic, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(M::Hyperbolic, Y, p, X, q, ::ParallelTransport)
    w = log(M, p, q)
    wn = norm(M, p, w)
    wn < eps(eltype(p + q)) && return copyto!(Y, X)
    return copyto!(Y, X - (inner(M, p, w, X) * (w + log(M, q, p)) / wn^2))
end


text:
"""
`AbstractExecutor` is an abstract type for executors

Executors are structs which are responsible of running
`Action` attached to a given `Job`
"""
abstract type AbstractExecutor end

include("debug.jl")

ext: jl


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

export normalize_text, isemoji
using Base.Unicode

#, language!
# using Languages
# using SnowballStemmer

# A symbol s in this list will be expanded to BLANK*s if the predecesor of s is neither s nor BLANK
# On changes from s to BLANK or [^s] it will produce also produce an extra BLANK
# Note that enabled del_punc will delete all these symbols without any of the previous expansions

const BLANK = ' '
const EMOJIS = Set([l[1] for l in readlines(joinpath(@__DIR__, "emojis.txt"))])

const RE_USER = r"""@[^;:,.@#&\\\-\"'/:\*\(\)\[\]\¿\?\¡\!\{\}~\<\>\|\s]+"""
const RE_URL = r"(http|ftp|https)://\S+"
const RE_NUM = r"\d+"

function isemoji(c::Char)
    c in EMOJIS
end

function _preprocessing(config::TextConfig, text)
    if config.lc
        text = lowercase(text)
    end

    if config.group_url
        text = replace(text, RE_URL => "_url")
    end

    if config.group_usr
        text = replace(text, RE_USER => "_usr")
    end

    if config.group_num
        text = replace(text, RE_NUM => "_num")
    end

    text
end

"""
    normalize_text(config::TextConfig, text::AbstractString, output::Vector{Char})

Normalizes a given text using the specified transformations of `config`
"""
function normalize_text(config::TextConfig, text::AbstractString, output::Vector{Char})
    text = _preprocessing(config, text)
    push!(output, BLANK)
    rep = 0

    @inbounds for u in Unicode.normalize(text, casefold=config.lc, stripmark=config.del_diac, stripcc=true, compat=true)
        isspace(u) && (u = BLANK)
        config.del_punc && ispunct(u) && !(u in ('@', '#', '_')) && (u = BLANK)
        config.group_emo && isemoji(u) && (u = '👾')
        rep = u === output[end] ? rep + 1 : 0
        config.del_dup && rep > 1 && continue
        push!(output, u)
    end

    push!(output, BLANK)
    output
end


text:
# check for method ambiguities
@test check_method_ambiguity_binary(intersection)
@test check_method_ambiguity_binary(issubset)
@test check_method_ambiguity_binary(is_intersection_empty)
@test check_method_ambiguity_binary(convex_hull)


text:
immutable NoBackend <: AbstractBackend end

const _backendType = Dict{Symbol, DataType}(:none => NoBackend)
const _backendSymbol = Dict{DataType, Symbol}(NoBackend => :none)
const _backends = Symbol[]
const _initialized_backends = Set{Symbol}()

backends() = _backends
backend_name() = CURRENT_BACKEND.sym
_backend_instance(sym::Symbol) = haskey(_backendType, sym) ? _backendType[sym]() : error("Unsupported backend $sym")

macro init_backend(s)
    str = lowercase(string(s))
    sym = symbol(str)
    T = symbol(string(s) * "Backend")
    esc(quote
        immutable $T <: AbstractBackend end
        export $sym
        $sym(; kw...) = (default(; kw...); backend(symbol($str)))
        backend_name(::$T) = symbol($str)
        push!(_backends, symbol($str))
        _backendType[symbol($str)] = $T
        _backendSymbol[$T] = symbol($str)
        include("backends/" * $str * ".jl")
    end)
end

@init_backend Immerse
@init_backend Gadfly
@init_backend PyPlot
@init_backend Qwt
@init_backend UnicodePlots
@init_backend Winston
@init_backend Bokeh
@init_backend Plotly
@init_backend PlotlyJS
@init_backend GR
@init_backend GLVisualize
@init_backend PGFPlots

include("backends/web.jl")
# include("backends/supported.jl")

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

plot(pkg::AbstractBackend; kw...)                       = error("plot($pkg; kw...) is not implemented")
plot!(pkg::AbstractBackend, plt::Plot; kw...)           = error("plot!($pkg, plt; kw...) is not implemented")
_update_plot(pkg::AbstractBackend, plt::Plot, d::KW)  = error("_update_plot($pkg, plt, d) is not implemented")
_update_plot_pos_size{P<:AbstractBackend}(plt::AbstractPlot{P}, d::KW) = nothing
subplot(pkg::AbstractBackend; kw...)                    = error("subplot($pkg; kw...) is not implemented")
subplot!(pkg::AbstractBackend, subplt::Subplot; kw...)  = error("subplot!($pkg, subplt; kw...) is not implemented")

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

type CurrentBackend
  sym::Symbol
  pkg::AbstractBackend
end
CurrentBackend(sym::Symbol) = CurrentBackend(sym, _backend_instance(sym))

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

function pickDefaultBackend()
    env_default = get(ENV, "PLOTS_DEFAULT_BACKEND", "")
    if env_default != ""
        try
            Pkg.installed(env_default)  # this will error if not installed
            sym = symbol(lowercase(env_default))
            if haskey(_backendType, sym)
                return backend(sym)
            else
                warn("You have set PLOTS_DEFAULT_BACKEND=$env_default but it is not a valid backend package.  Choose from:\n\t",
                     join(sort(_backends), "\n\t"))
            end
        catch
            warn("You have set PLOTS_DEFAULT_BACKEND=$env_default but it is not installed.")
        end
    end

    # the ordering/inclusion of this package list is my semi-arbitrary guess at
    # which one someone will want to use if they have the package installed...accounting for
    # features, speed, and robustness
    for pkgstr in ("PyPlot", "GR", "PlotlyJS", "Immerse", "Gadfly", "UnicodePlots")
        if Pkg.installed(pkgstr) != nothing
            return backend(symbol(lowercase(pkgstr)))
        end
    end

    # the default if nothing else is installed
    backend(:plotly)
end

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

"""
Returns the current plotting package name.  Initializes package on first call.
"""
function backend()

  global CURRENT_BACKEND
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end

  sym = CURRENT_BACKEND.sym
  if !(sym in _initialized_backends)

    # initialize
    println("[Plots.jl] Initializing backend: ", sym)

    inst = _backend_instance(sym)
    try
      _initialize_backend(inst)
    catch err
      warn("Couldn't initialize $sym.  (might need to install it?)")
      rethrow(err)
    end

    push!(_initialized_backends, sym)

  end
  CURRENT_BACKEND.pkg
end

"""
Set the plot backend.
"""
function backend(pkg::AbstractBackend)
  CURRENT_BACKEND.sym = backend_name(pkg)
  CURRENT_BACKEND.pkg = pkg
end

function backend(modname::Symbol)
  CURRENT_BACKEND.sym = modname
  CURRENT_BACKEND.pkg = _backend_instance(modname)
end

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

supportedAxes(::AbstractBackend) = [:left]
supportedTypes(::AbstractBackend) = []
supportedStyles(::AbstractBackend) = [:solid]
supportedMarkers(::AbstractBackend) = [:none]
supportedScales(::AbstractBackend) = [:identity]
subplotSupported(::AbstractBackend) = false
stringsSupported(::AbstractBackend) = false
nativeImagesSupported(::AbstractBackend) = false

supportedAxes() = supportedAxes(backend())
supportedTypes() = supportedTypes(backend())
supportedStyles() = supportedStyles(backend())
supportedMarkers() = supportedMarkers(backend())
supportedScales() = supportedScales(backend())
subplotSupported() = subplotSupported(backend())
stringsSupported() = stringsSupported(backend())
nativeImagesSupported() = nativeImagesSupported(backend())

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


text:
module wilsonmod

include("general.jl")
using Unitful
using NLsolve

∑(x)=sum(x)
λ(T::Quantity,a,b,c,i::Int,j::Int) = a[i,j]u"cal/mol"+ b[i,j]u"cal/(mol*K)"*T + c[i,j]u"cal/(mol*K^2)"*T^2
λ(T::Float64,a,b,c,i::Int,j::Int) = a[i,j]+ b[i,j]*T + c[i,j]*T^2
Λ(T::Quantity,a,b,c,v,i::Int,j::Int) = (v[j]/v[i])*exp(-λ(T,a,b,c,i,j)/(1.98721u"cal/(mol*K)"*T))
Λ(T::Float64,a,b,c,v,i::Int,j::Int) = (v[j]/v[i])*exp(-λ(T,a,b,c,i,j)/(1.98721*T))
γ(T,x,a,b,c,v,i,n)=exp(-log(∑(x[j]*Λ(T,a,b,c,v,i,j) for j=1:n)) + 1 -∑((x[k]*Λ(T,a,b,c,v,k,i))/(∑(x[j]*Λ(T,a,b,c,v,k,j) for j=1:n)) for k=1:n))
P(T,x,CA,a,b,c,v,i,n) = x[i]*γ(T,x,a,b,c,v,i,n)*general.P_Antoine(CA,T,i)
Pₜₒₜ(T,x,CA,a,b,c,v,n) = ∑(P(T,x,CA,a,b,c,v,i,n) for i=1:n)
y_calc(T,x,CA,a,b,c,v,i,n) = P(T,x,CA,a,b,c,v,i,n)/Pₜₒₜ(T,x,CA,a,b,c,v,n)

"""
    pxy(n,T,CA,a,b,c,v,xx;uni=u"Torr")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `T :: Float` Es la temperatura del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `a :: n×n Array{Float,2}` Parametros de interacción a en ``\\frac{cal}{mol}``
- `b :: n×n Array{Float,2}` Parametros de interacción b en ``\\frac{cal}{mol K}``
- `c :: n×n Array{Float,2}` Parametros de interacción c en ``\\frac{cal}{mol K^2}``
- `v :: n-element Array{Float,1}` Volumen molar
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en Torr

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la presión `P` en Torr (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `P :: m-element Array{Float64,1}` presión
"""
function pxy(n::Int,T,CA::Array,a::Array,b::Array,c::Array,v::Vector,xx;uni=u"Torr")
    T = T|> u"K"
    xx = general.nc(xx...)
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    P = Vector{Float64}(undef,m)*uni
    for i=1:m
        x = xx[:,i]
        for j=1:n
            y[j,i] = y_calc(T,x,CA,a,b,c,v,j,n)
        end
        P[i] = Pₜₒₜ(T,x,CA,a,b,c,v,n) |> uni
    end

    return xx,y,P

end

"""
    txy(n,Pobj,CA,a,b,c,v,xx;uni=u"K")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `Pobj :: Float` Es la presión del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `a :: n×n Array{Float,2}` Parametros de interacción a en ``\\frac{cal}{mol}``
- `b :: n×n Array{Float,2}` Parametros de interacción b en ``\\frac{cal}{mol K}``
- `c :: n×n Array{Float,2}` Parametros de interacción c en ``\\frac{cal}{mol K^2}``
- `v :: n-element Array{Float,1}` Volumen molar
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en grados kelvin

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la temperatura `T` en grados kelvin (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `T :: m-element Array{Float64,1}` temperatura
"""
function txy(n::Int,Pobj,CA::Array,a::Array,b::Array,c::Array,v::Vector,xx;uni=u"K")
    Pobj = Pobj |> u"Torr"
    xx = general.nc(xx...)
    #=
    if n==3
        xx = general.nc_ord(xx...)
    else
        xx = general.nc(xx...)
    end=#
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    T = Vector{Float64}(undef,m)*u"K"
    Tᵢ = 300.0
    #temperatura inicial
    for i=1:n
        if xx[i,1]==1
            T[1] = general.T_Antoine(CA,Pobj,i)
            Tᵢ  = ustrip(u"K",T[1])
        end
    end
    for i = 1:m
        x = xx[:,i]
        P2ₜₒₜ(T) = Pₜₒₜ(T,x,CA,a,b,c,v,n) - ustrip(u"Torr",Pobj)
        #println(i)
        Pₛₒₗ = nlsolve(n_ary(P2ₜₒₜ),[Tᵢ])
        Tᵢ = Pₛₒₗ.zero[1]
        T[i] = Tᵢ*u"K"
        for j = 1:n
            y[j,i] = y_calc(T[i],x,CA,a,b,c,v,j,n)
        end
    end

    return xx,y,T .|> uni

end
end  # module

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
using Metalhead, Test
using Flux

@testset "MLPMixer" begin
	@testset for mode in [:small, :base, :large, :huge]
		@testset for drop_path_rate in [0.0, 0.5, 0.99]
			m = MLPMixer(mode; drop_path_rate)
			@test size(m(rand(Float32, 224, 224, 3, 2))) == (1000, 2)
			@test_skip gradtest(m, rand(Float32, 224, 224, 3, 2))
			GC.gc()
		end
	end
end

@testset "ResMLP" begin
	@testset for mode in [:small, :base, :large, :huge]
		@testset for drop_path_rate in [0.0, 0.5, 0.99]
			m = ResMLP(mode; drop_path_rate)
			@test size(m(rand(Float32, 224, 224, 3, 2))) == (1000, 2)
			@test_skip gradtest(m, rand(Float32, 224, 224, 3, 1))
			GC.gc()
		end
	end
end

@testset "gMLP" begin
	@testset for mode in [:small, :base, :large, :huge]
		@testset for drop_path_rate in [0.0, 0.5, 0.99]
			m = gMLP(mode; drop_path_rate)
			@test size(m(rand(Float32, 224, 224, 3, 2))) == (1000, 2)
			@test_skip gradtest(m, rand(Float32, 224, 224, 3, 2))
			GC.gc()
		end
	end
end


text:
#for measuring consistency of a periodic 1d timeseries
#returns a sample vector one cycle in duration where each sample is equal to the largest absolue difference measured between pairs of corresponding samples
#can choose to ignore the first ignore_ncycs cycles
function largest_cycle_diff(a::T, cycle_dur::Union{Int, HasTimeUnits}, delay_cyc1::Int=0; ignore_ncycs::Int=0) where {T<:AbstractVector}
    cyc_mat = get_cycles(a, cycle_dur, delay_cyc1)
    largest_cycle_diff(view(cyc_mat, :, ignore_ncycs+1:size(cyc_mat,2)))
end

largest_cycle_diff(a::T) where {T<:AbstractMatrix} = maximum(a, dims=2) .- minimum(a, dims=2)
largest_cycle_diff(sig::ImagineSignal, cycle_dur, delay_cyc1::Int; ignore_ncycs::Int=0) = largest_cycle_diff(get_samples(sig), cycle_dur, delay_cyc1; ignore_ncycs=ignore_ncycs)
largest_cycle_diff(sig::ImagineSignal, cycle_dur, delay_cyc1::HasTimeUnits; ignore_ncycs::Int=0) = largest_cycle_diff(get_samples(sig), cycle_dur, ImagineInterface.calc_num_samps(delay_cyc1, samprate(sig)); ignore_ncycs=ignore_ncycs)
largest_cycle_diff(sig::ImagineSignal, cycle_dur; ignore_ncycs::Int=0) = largest_cycle_diff(get_samples(sig), cycle_dur, 0; ignore_ncycs=ignore_ncycs)

#same as largest_cycle_diff but returns per-sample std
function largest_cycle_std(a::T, cycle_dur::Union{Int, HasTimeUnits}, delay_cyc1::Int; ignore_ncycs::Int=0) where {T<:AbstractVector}
    cyc_mat = get_cycles(a, cycle_dur, delay_cyc1)
    largest_cycle_std(view(cyc_mat, :, ignore_ncycs+1:size(cyc_mat,2)))
end

largest_cycle_std(a::T) where {T<:AbstractMatrix} = std(a, 2)
largest_cycle_std(sig::ImagineSignal, cycle_dur::HasTimeUnits, delay_cyc1::Int; ignore_ncycs::Int=0) = largest_cycle_std(get_samples(sig), cycle_dur, delay_cyc1; ignore_ncycs=ignore_ncycs)
largest_cycle_std(sig::ImagineSignal, cycle_dur::HasTimeUnits, delay_cyc1::HasTimeUnits; ignore_ncycs::Int=0) = largest_cycle_std(get_samples(sig), cycle_dur, ImagineInterface.calc_num_samps(delay_cyc1, samprate(sig)); ignore_ncycs=ignore_ncycs)
largest_cycle_std(sig::ImagineSignal, cycle_dur::HasTimeUnits; ignore_ncycs::Int=0) = largest_cycle_std(get_samples(sig), cycle_dur, 0; ignore_ncycs=ignore_ncycs)

#Reshape a periodic 1D signal into a matrix with sample index on the first axis and cycle index on the second
get_cycles(sig::ImagineSignal, cycle_dur, delay::HasTimeUnits) = get_cycles(get_samples(sig), cycle_dur, delay)
get_cycles(sig::ImagineSignal, cycle_dur, delay_nsamps::Int) = get_cycles(get_samples(sig), cycle_dur, delay_nsamps)
get_cycles(sig::ImagineSignal, cycle_dur) = get_cycles(get_samples(sig), cycle_dur, 0)

isapprox_int(x) = isapprox(round(Int, x), x)

get_cycles(a::AxisArray{T,1,D,Ax}, cycle_dur::HasTimeUnits; delay::HasTimeUnits = 0.0s) where {T,D,Ax} = get_cycles(a, cycle_dur, round(Int, delay/(step(axisvalues(a)[1]))))

function get_cycles(a::AxisArray{T,1,D,Ax}, cycle_dur::HasTimeUnits, delay_nsamps::Int) where {T,D,Ax}
    samp_dur = step(axisvalues(a)[1])
    nsamps = cycle_dur/samp_dur
    if isapprox_int(nsamps)
        nsamps = round(Int,nsamps) #special case
    end
    get_cycles(a.data, nsamps, delay_nsamps)
end

#general case (allocates new array)
#interpolates new cycle vector to honor the cycle time
#if input vector doesn't contain a whole number of cycles then omit a partial cycle at the end
function get_cycles(a::T, cycle_nsamps::Float64, delay_nsamps::Int = 0) where {T<:AbstractVector}
    la = length(a) - delay_nsamps
    if la < cycle_nsamps
        error("Insufficient sample count for the specified cycle time")
    end
    ncycs = la/cycle_nsamps
    if isapprox_int(ncycs)
        ncycs = round(Int, ncycs)
    else
        ncycs = floor(Int, ncycs)
    end
    nsamps_interp = ceil(Int, cycle_nsamps)
    output = zeros(eltype(a), nsamps_interp, ncycs)
    a_itp = extrapolate(interpolate(a, BSpline(Linear())), Flat())
    step_size = cycle_nsamps/nsamps_interp
    for i = 1:ncycs
        #output[:,i] = a_itp[linspace(1+(i-1)*cycle_nsamps, i*cycle_nsamps, nsamps_interp)] #not sure why this fails for certain float ranges
        cyc_inds = range(delay_nsamps+1+1+(i-1)*cycle_nsamps, stop=i*cycle_nsamps, length=nsamps_interp)
        for (ii,ival) in enumerate(cyc_inds)
            output[ii,i] = a_itp(ival)
        end
    end
    return output
end

#special case (in-place reshape)
function get_cycles(a::T, cycle_nsamps::Int, delay_nsamps::Int = 0) where {T<:AbstractVector}
    la = length(a) - delay_nsamps
    d, r = divrem(la, cycle_nsamps)
    if la < cycle_nsamps
        error("Insufficient sample count for the specified cycle time")
    else
        a = view(a, (delay_nsamps+1):(delay_nsamps+d*cycle_nsamps))
    end
    return reshape(a, cycle_nsamps, d)
end


text:
using NLSolvers
NLE_PROBS = Dict()

# Rosenbrock
# Source: MINPACK
NLE_PROBS["rosenbrock"] = Dict()
NLE_PROBS["rosenbrock"]["array"] = Dict()
function F_rosenbrock!(Fx, x)
    Fx[1] = 1 - x[1]
    Fx[2] = 10(x[2] - x[1]^2)
    return Fx
end
function J_rosenbrock!(Jx, x)
    Jx[1, 1] = -1
    Jx[1, 2] = 0
    Jx[2, 1] = -20 * x[1]
    Jx[2, 2] = 10
    return Jx
end
function FJ_rosenbrock!(Fx, Jx, x)
    F_rosenbrock!(Fx, x)
    J_rosenbrock!(Jx, x)
    Fx, Jx
end
function Jvop_rosenbrock!(x)
    function JacV(Fv, v)
        Fv[1] = -1 * v[1]
        Fv[2,] = -20 * x[1] * v[1] + 10 * v[2]
    end
    LinearMap(JacV, length(x))
end

NLE_PROBS["rosenbrock"]["array"]["x0"] = [-1.2, 1.0]
NLE_PROBS["rosenbrock"]["array"]["mutating"] = NLSolvers.VectorObjective(
    F_rosenbrock!,
    J_rosenbrock!,
    FJ_rosenbrock!,
    Jvop_rosenbrock!,
)

function F_powell_singular!(x::Vector, Fx::Vector, Jx::Union{Nothing,Matrix} = nothing)
    if !(Fx isa Nothing)
        Fx[1] = x[1] + 10x[2]
        Fx[2] = sqrt(5) * (x[3] - x[4])
        Fx[3] = (x[2] - 2x[3])^2
        Fx[4] = sqrt(10) * (x[1] - x[4])^2
    end
    if !(Jx isa Nothing)
        fill!(Jx, 0)
        Jx[1, 1] = 1
        Jx[1, 2] = 10
        Jx[2, 3] = sqrt(5)
        Jx[2, 4] = -Jx[2, 3]
        Jx[3, 2] = 2(x[2] - 2x[3])
        Jx[3, 3] = -2Jx[3, 2]
        Jx[4, 1] = 2sqrt(10) * (x[1] - x[4])
        Jx[4, 4] = -Jx[4, 1]
    end
    objective_return(Fx, Jx)
end

function Jvop_powell_singular(x)
    function JacV(Fv, v)
        Fv[1] = v[1] + 10 * v[2]
        Fv[2] = (v[3] - v[4]) * sqrt(5)
        xx23 = 2 * x[2] - 4 * x[3]
        Fv[3] = v[2] * xx23 - v[3] * xx23 * 2
        xx41 = 2 * sqrt(10) * (x[1] - x[4])
        Fv[4] = v[1] * xx41 - v[4] * xx41
    end
    LinearMap(JacV, length(x))
end

#=
NLE_PROBS["quantile"] = Dict()
NLE_PROBS["quantile"]["number"] = Dict()
@inline quantile_f(Fx, x) = log(max(x, 0.000001))
@inline quantile_j(Jx, x) = 1.0/x
@inline function quantile_fj(Fx, Jx, x)
    Fx = quantile_f(Fx, x)
    Jx = quantile_j(Jx, x)
    Fx, Jx
end
NLE_PROBS["quantile"]["number"]["x0"] = 0.5
NLE_PROBS["quantile"]["number"]["mutating"] = quantobj
function f()
    quantile_f(Fx, x) = log(max(x, 0.000001))
    quantile_j(Jx, x) = 1.0/x
    function quantile_fj(Fx, Jx, x)
       Fx = quantile_f(Fx, x)
       Jx = quantile_j(Jx, x)
       Fx, Jx
    end
    quantobj = NLSolvers.VectorObjective(quantile_f, quantile_j, quantile_fj, nothing)

    quantproblem = NEqProblem(quantobj, nothing, NLSolvers.Euclidean(0), NLSolvers.OutOfPlace())

    method = LineSearch(Newton())
    options = NEqOptions()
    state = NLSolvers.init(quantproblem, method, 3.0)
    @time res = solve(quantproblem, 0.4, method, options, state)
    @time res = solve(quantproblem, 0.4, method, options, state)
end

quantile_f(Fx, x) = log(max(x, 0.000001))
quantile_j(Jx, x) = 1.0/x
function quantile_fj(Fx, Jx, x)
   Fx = quantile_f(Fx, x)
   Jx = quantile_j(Jx, x)
   Fx, Jx
end
const quantobj = NLSolvers.VectorObjective(quantile_f, quantile_j, quantile_fj, nothing)

const quantproblem = NEqProblem(quantobj, nothing, NLSolvers.Euclidean(0), NLSolvers.OutOfPlace())

method = LineSearch(Newton())
options = NEqOptions()
state = NLSolvers.init(quantproblem, method, 3.0)
@time res = solve(quantproblem, 0.4, method, options, state)
@time res = solve(quantproblem, 0.4, method, options, state)
=#

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
struct LabeledGraph
    n :: Int32 # |V|
    m :: Int32 # |E|
    V :: Array{Int32, 1} # V[i] = Real Index of node i
    L :: Array{Int32, 1} # L[i] = Label of node i
    E :: Array{Tuple{Int32, Int32}, 1} # (u, v) in Edge Set
end

function readLabeledGraph(graphName; dataDir = "../data/") # Read a labeled graph
    # Initialized
    n = 0
    V = Vector{Int32}()
    id = Dict{Int32, Int32}()
    edge = Set{Tuple{Int32, Int32}}()

    getid(x) = begin
        if !haskey(id, x)
            n += 1
            push!(V, x)
            id[x] = n
            return n
        end
        return id[x]
    end

    # read graph from edge list file
    open(dataDir*graphName*".edgelist") do f1
        for line in eachline(f1)
            buf = split(line)
            if size(buf, 1) != 2
                continue
            end
            u = parse(Int32, buf[1])
            v = parse(Int32, buf[2])
            u1 = getid(u)
            v1 = getid(v)
            if u1 == v1
                continue
            end
            if u1 > v1
                u1, v1 = v1, u1
            end
            push!(edge, (u1, v1))
        end
    end

    # read labels
    L = zeros(Int32, n)
    open(dataDir*graphName*".group") do f2
        for line in eachline(f2)
            buf = split(line)
            if size(buf, 1) != 2
                continue
            end
            x = id[parse(Int32, buf[1])]
            y = parse(Int32, buf[2])
            L[x] = y
        end
    end

    m = length(edge)
    E = Array{Tuple{Int32, Int32}, 1}(undef, m)

    i = 0
    for (u, v) in edge
        i = i + 1
        E[i] = (u, v)
    end

    return LabeledGraph(n, m, V, L, E)
end

function readUnLabeledGraph(graphName; dataDir = "../data/") # Read a graph without labels
    # Initialized
    n = 0
    V = Vector{Int32}()
    id = Dict{Int32, Int32}()
    edge = Set{Tuple{Int32, Int32}}()

    getid(x) = begin
        if !haskey(id, x)
            n += 1
            push!(V, x)
            id[x] = n
            return n
        end
        return id[x]
    end

    # read graph from edge list file
    open(dataDir*graphName*".txt") do f1
        for line in eachline(f1)
            buf = split(line)
            if size(buf, 1) != 2
                continue
            end
            u = parse(Int32, buf[1])
            v = parse(Int32, buf[2])
            u1 = getid(u)
            v1 = getid(v)
            if u1 == v1
                continue
            end
            if u1 > v1
                u1, v1 = v1, u1
            end
            push!(edge, (u1, v1))
        end
    end

    L = zeros(Int32, n)

    m = length(edge)
    E = Array{Tuple{Int32, Int32}, 1}(undef, m)

    i = 0
    for (u, v) in edge
        i = i + 1
        E[i] = (u, v)
    end

    return LabeledGraph(n, m, V, L, E)
end

function adjacentList(G)
    g = Array{Array{Int32, 1}, 1}(undef, G.n)
    foreach(i -> g[i] = [], 1 : G.n)
    for (u, v) in G.E
        push!(g[u], v)
        push!(g[v], u)
    end
    return g
end

function degrees(G)
    d = zeros(Int, G.n)
    for (u, v) in G.E
        d[u] += 1
        d[v] += 1
    end
    return d
end


text:
using CodecBzip2
using CodecZstd
using CodecXz
using CodecZlib
using TranscodingStreams

decompression_codecs = Dict(
    "bz2" => Bzip2DecompressorStream,
    "zst" => ZstdDecompressorStream,
    "xz" => XzDecompressorStream,
    "gz" => GzipDecompressorStream,
)

get_decompressor_stream(filename::String) =
    get(decompression_codecs, split(filename, ".") |> last, NoopStream)

function getfile(path::String)
    if isdir(path)
        files = readdir(path, join = true)
        if length(files) > 1
            @warn "More than one file in directory ($path). Using first file."
        end
        return files |> first
    else
        return path
    end
end


text:
export eigs,
       entropy,
       polar,
       random_orthog,
       random_unitary,
       Spectrum,
       svd,
       truncerror

#
# Linear Algebra of order 2 NDTensors
#
# Even though DenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray

function Base.:*(T1::Tensor{ElT1,2,StoreT1},
                 T2::Tensor{ElT2,2,StoreT2}) where
                                       {ElT1,StoreT1<:Dense,IndsT1,
                                        ElT2,StoreT2<:Dense,IndsT2}
  RM = matrix(T1)*matrix(T2)
  indsR = (ind(T1,1), ind(T2,2))
  return tensor(Dense(vec(RM)), indsR)
end

function LinearAlgebra.exp(T::DenseTensor{ElT,2}) where {ElT<:Union{Real,Complex}}
  expTM = exp(matrix(T))
  return tensor(Dense(vec(expTM)),inds(T))
end

function LinearAlgebra.exp(T::Hermitian{ElT,
                                        <:DenseTensor{ElT,2}}) where {ElT<:Union{Real,Complex}}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(matrix(T)))
  return tensor(Dense(vec(expTM)),inds(T))
end

"""
  Spectrum
contains the (truncated) density matrix eigenvalue spectrum which is computed during a
decomposition done by `svd` or `eigen`. In addition stores the truncation error.
"""
struct Spectrum{VecT <: Union{AbstractVector, Nothing}}
  eigs::VecT
  truncerr::Float64
end

eigs(s::Spectrum) = s.eigs
truncerror(s::Spectrum) = s.truncerr

function entropy(s::Spectrum)
  S = 0.0
  eigs_s = eigs(s)
  isnothing(eigs_s) && error("Spectrum does not contain any eigenvalues, cannot compute the entropy")
  for p in eigs_s
    p > 1e-13 && (S -= p*log(p))
  end
  return S
end

"""
    svd(T::DenseTensor{<:Number,2}; kwargs...)

svd of an order-2 DenseTensor
"""
function LinearAlgebra.svd(T::DenseTensor{ElT,2,IndsT};
                           kwargs...) where {ElT,IndsT}
  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  #
  # Keyword argument deprecations
  #
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs,
                              :absoluteCutoff,
                              use_absolute_cutoff)
  end

  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs,
                              :doRelCutoff,
                              use_relative_cutoff)
  end

  if haskey(kwargs, :fastsvd) || haskey(kwargs, :fastSVD)
    error("In svd, fastsvd/fastSVD keyword arguments are removed in favor of alg, see documentation for more details.")
  end

  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  use_absolute_cutoff::Bool = get(kwargs,
                                  :use_absolute_cutoff,
                                  use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs,
                                  :use_relative_cutoff,
                                  use_relative_cutoff)
  alg::String = get(kwargs, :alg, "recursive")

  if alg == "divide_and_conquer"
    MU,MS,MV = svd(matrix(T); alg = LinearAlgebra.DivideAndConquer())
  elseif alg == "qr_iteration"
    MU,MS,MV = svd(matrix(T); alg = LinearAlgebra.QRIteration())
  elseif alg == "recursive"
    MU,MS,MV = svd_recursive(matrix(T))
  else
    error("svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.")
  end
  conj!(MV)

  P = MS .^ 2
  if truncate
    truncerr, _ = truncate!(P; mindim = mindim,
                               maxdim = maxdim,
                               cutoff = cutoff,
                               use_absolute_cutoff = use_absolute_cutoff,
                               use_relative_cutoff = use_relative_cutoff)
  else
    truncerr = 0.0
  end
  spec = Spectrum(P, truncerr)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    resize!(MS,dS)
    MV = MV[:,1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T,1),u))
  Sinds = IndsT((u,v))
  Vinds = IndsT((ind(T,2),v))
  U = tensor(Dense(vec(MU)),Uinds)
  S = tensor(Diag(MS),Sinds)
  V = tensor(Dense(vec(MV)),Vinds)
  return U,S,V,spec
end

function LinearAlgebra.eigen(T::Hermitian{ElT,<:DenseTensor{ElT,2,IndsT}};
                             kwargs...) where {ElT<:Union{Real,Complex},IndsT}
  # Keyword argument deprecations
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs,
                              :absoluteCutoff,
                              use_absolute_cutoff)
  end
  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs,
                              :doRelCutoff,
                              use_relative_cutoff)
  end

  truncate = haskey(kwargs,:maxdim) || haskey(kwargs,:cutoff)
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  use_absolute_cutoff::Bool = get(kwargs,
                                  :use_absolute_cutoff,
                                  use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs,
                                  :use_relative_cutoff,
                                  use_relative_cutoff)

  DM, VM = eigen(matrix(T))

  # Sort by largest to smallest eigenvalues
  p = sortperm(DM; rev = true)
  DM = DM[p]
  VM = VM[:,p]

  if truncate
    truncerr,_ = truncate!(DM;maxdim=maxdim,
                              cutoff=cutoff,
                              use_absolute_cutoff=use_absolute_cutoff,
                              use_relative_cutoff=use_relative_cutoff)
    dD = length(DM)
    if dD < size(VM, 2)
      VM = VM[:, 1:dD]
    end
  else
    dD = length(DM)
    truncerr = 0.0
  end
  spec = Spectrum(DM,truncerr)

  # Make the new indices to go onto V
  l = eltype(IndsT)(dD)
  r = eltype(IndsT)(dD)
  Vinds = IndsT((dag(ind(T, 2)), dag(r)))
  Dinds = IndsT((l, dag(r)))
  V = tensor(Dense(vec(VM)), Vinds)
  D = tensor(Diag(DM), Dinds)
  return D, V, spec
end

"""
    random_orthog(n::Int,m::Int)

Return a random matrix O of dimensions (n,m)
such that if n >= m, transpose(O)*O is the
identity, or if m > n O*transpose(O) is the
identity.
"""
function random_orthog(n::Int,m::Int)::Matrix
  #TODO: optimize by using recursive Householder algorithm?
  if n < m
    return transpose(random_orthog(m,n))
  end
  F = qr(randn(n,m))
  Q = Matrix(F.Q)
  for c=1:size(Q,2)
    if real(F.R[c,c]) < 0.0
      Q[:,c] *= -1
    end
  end
  return Q
end

"""
    random_unitary(n::Int,m::Int)::Matrix

Return a random complex matrix U of dimensions (n,m)
such that if n >= m, U'*U is the
identity, or if m > n U*U' is the
identity.

Sampling is based on https://arxiv.org/abs/math-ph/0609050
such that in the case `n==m`, the unitary matrix will be sampled
according to the Haar measure.
"""
function random_unitary(n::Int,m::Int)
  if n < m
    return Matrix(random_unitary(m,n)')
  end
  F = qr(randn(ComplexF64,n,m))
  Q = Matrix(F.Q)
  A = diag(Matrix(F.R))
  A ./= abs.(A)
  for c in 1:size(Q,2)
      Q[:,c] *= A[c]
  end
  return Q
end

function qr_positive(M::AbstractMatrix)
  sparseQ,R = qr(M)
  Q = convert(Matrix,sparseQ)
  nc = size(Q,2)
  for c=1:nc
    if real(R[c,c]) < 0.0
      R[c,c:end] *= -1
      Q[:,c] *= -1
    end
  end
  return (Q,R)
end

function LinearAlgebra.eigen(T::DenseTensor{ElT,2,IndsT};
                             kwargs...) where {ElT<:Union{Real,Complex},IndsT}
  # Keyword argument deprecations
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs,
                              :absoluteCutoff,
                              use_absolute_cutoff)
  end
  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs,
                              :doRelCutoff,
                              use_relative_cutoff)
  end

  truncate = haskey(kwargs,:maxdim) || haskey(kwargs,:cutoff)
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  use_absolute_cutoff::Bool = get(kwargs,
                                  :use_absolute_cutoff,
                                  use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs,
                                  :use_relative_cutoff,
                                  use_relative_cutoff)

  DM, VM = eigen(matrix(T))

  # Sort by largest to smallest eigenvalues
  #p = sortperm(DM; rev = true)
  #DM = DM[p]
  #VM = VM[:,p]

  if truncate
    truncerr,_ = truncate!(DM;maxdim=maxdim,
                              cutoff=cutoff,
                              use_absolute_cutoff=use_absolute_cutoff,
                              use_relative_cutoff=use_relative_cutoff)
    dD = length(DM)
    if dD < size(VM, 2)
      VM = VM[:, 1:dD]
    end
  else
    dD = length(DM)
    truncerr = 0.0
  end
  spec = Spectrum(abs.(DM), truncerr)

  i1, i2 = inds(T)

  # Make the new indices to go onto D and V
  l = typeof(i1)(dD)
  r = dag(sim(l))
  Dinds = (l, r)
  Vinds = (dag(i2), r)
  D = complex(tensor(Diag(DM), Dinds))
  V = complex(tensor(Dense(vec(VM)), Vinds))
  return D, V, spec
end

function LinearAlgebra.qr(T::DenseTensor{ElT,2,IndsT}
                          ;kwargs...) where {ElT,IndsT}
  positive = get(kwargs,:positive,false)
  # TODO: just call qr on T directly (make sure
  # that is fast)
  if positive
    QM,RM = qr_positive(matrix(T))
  else
    QM,RM = qr(matrix(T))
  end
  # Make the new indices to go onto Q and R
  q,r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = tensor(Dense(vec(Matrix(QM))),Qinds)
  R = tensor(Dense(vec(RM)),Rinds)
  return Q,R
end

# TODO: support alg keyword argument to choose the svd algorithm
function polar(T::DenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM,RM = polar(matrix(T))
  dim = size(QM,2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = tensor(Dense(vec(QM)),Qinds)
  R = tensor(Dense(vec(RM)),Rinds)
  return Q,R
end


text:
"""
ドットプロットでは，データ点が重ならないようにランダムな数値（Jitter）を加えているので，その並びに何らかの意味がありそうに見えてしまう。Python の seaborn にある swarmplot() を Julia から使えるようになっているが，これも同じ特質を持っている。

そのようなことを避けるために
https://blog.goo.ne.jp/r-de-r/e/237f8bdfda925c5d724604aa9eabc8e3
にプログラムを呈示している。
"""

using FreqTables, Plots
function simple_dotplot(x, y; accu=0, stp=0, logflag=false, simple=false,
            symmetrical=true, sizew=600, sizeh=400)
    function logaxis(z)
        pos = zeros(100)
        minval = minimum(z)
        maxval = maximum(z)
        logmaxval = log10(maxval)
        logminval = log10(minval)
        width = logmaxval - logminval
        factor = logmaxval < 3.5 ? 15 : 8
        w = width / factor
        st = floor(logminval)
        delta = 10^st
        ic = 0
        while true
            for i = 1:9
                v = i * delta
                if minval - delta <= v <= maxval + delta
                    if (i == 1 || ic == 0) || (ic > 0 && log10(v) - pos[ic] > w)
                        if i == 1 && ic > 0 && log10(v) - pos[ic] <= w
                            ic -= 1
                        end
                        ic += 1
                        pos[ic] = log10(v)
                        prev = v
                    end
                end
            end
            delta *= 10
            if delta > maxval
                break
            end
        end
        return ic > 3 ? pos[2:ic-1] : [logminval, median([logminval, logmaxval]), logmaxval]
    end
    if logflag == true
        y0 = y
        y = log10.(y);
    end
    minx = minimum(y)
    maxx = maximum(y)
    rng = maxx - minx
    minx = minx - rng * 20 / sizeh
    maxx = maxx + rng * 20 / sizeh
    if accu == 0
        accu = (maxx - minx) / 50
    end
    if stp == 0
        stp = 30 / sizew
    end
    y = round.(y ./ accu, digits=0) * accu
    table = freqtable(y, x)
    yname, xname = names(table)
    nvalue, ngroup = size(table)
    # pyplot()
    p = plot(tick_direction=:out, xticks=false, grid=false,
                size=(sizew, sizeh), xlims = (0, length(xname)+1),
                ylims=(minx, maxx), label="")
    xticks!(1:ngroup, xname)
    for i = 1:ngroup
        for j = 1:nvalue
            if table[j, i] != 0
                n = table[j, i]
                if n == 1
                    height, position = [yname[j]], [i]
                else
                    height = repeat([yname[j]], n)
                    position = collect(range(i-stp*(n-1)/2, i+stp*(n-1)/2, length=n))
                end
                p = scatter!(position, height, labels="", alpha=0.7,
                                markercolor=i, markerstrokecolor=i)
            end
        end
    end
    if logflag
        pos = logaxis(y0)
        yticks!(pos, string.(round.(10 .^ pos, digits = 1)))
    end
    println("accu = $accu, stp = $stp")
    # display(p)
    p
end;

"""
使用法

julia> using RDatasets

julia> airquality = RDatasets.dataset("datasets", "airquality");

julia> x = string.(airquality.Month); # 型変換が必要

julia> y = airquality.Wind;

julia> simple_dotplot(x, y, accu=0.105, stp=0.05)
accu = 0.105, stp = 0.05

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
__precompile__()

module ExpressionPatterns
export @metadispatch, @metadestruct

path(name) = joinpath(name, "$name.jl")

include(path("Helper"))
include(path("PatternStructure"))
include(path("Analyzer"))
include(path("Matching"))
include(path("Destructuring"))
include(path("Dispatch"))
include(path("Docs"))

import .Dispatch: @metadestruct, @metadispatch
end


text:
const UNINITIALIZED_DATETIME = Dates.DateTime(Dates.Minute(0))
const UNINITIALIZED_LENGTH = 0
const UNINITIALIZED_PERIOD = Dates.Period(Dates.Minute(0))

mutable struct ForecastParameters <: InfrastructureSystemsType
    horizon::Int
    initial_timestamp::Dates.DateTime
    interval::Dates.Period
    count::Int
end

function ForecastParameters(;
    horizon=UNINITIALIZED_LENGTH,
    initial_timestamp=UNINITIALIZED_DATETIME,
    interval=UNINITIALIZED_PERIOD,
    count=UNINITIALIZED_LENGTH,
)
    return ForecastParameters(horizon, initial_timestamp, interval, count)
end

function check_params_compatibility(params::ForecastParameters, other::ForecastParameters)
    _is_uninitialized(params) && return true

    if other.count != params.count
        throw(
            ConflictingInputsError(
                "forecast count $(other.count) does not match system count $(params.count)",
            ),
        )
    end

    if other.horizon != params.horizon
        throw(
            ConflictingInputsError(
                "forecast horizon $(other.horizon) does not match system horizon $(params.horizon)",
            ),
        )
    end

    if other.initial_timestamp != params.initial_timestamp
        throw(
            ConflictingInputsError(
                "forecast initial_timestamp $(other.initial_timestamp) does not match system " *
                "initial_timestamp $(params.initial_timestamp)",
            ),
        )
    end

    return
end

function _is_uninitialized(params::ForecastParameters)
    return params.horizon == UNINITIALIZED_LENGTH &&
           params.initial_timestamp == UNINITIALIZED_DATETIME &&
           params.interval == UNINITIALIZED_PERIOD &&
           params.count == UNINITIALIZED_LENGTH
end

function reset_info!(params::ForecastParameters)
    params.horizon = UNINITIALIZED_LENGTH
    params.initial_timestamp = UNINITIALIZED_DATETIME
    params.interval = UNINITIALIZED_PERIOD
    params.count = UNINITIALIZED_LENGTH
    return
end

function get_forecast_initial_times(params::ForecastParameters)
    return get_initial_times(params.initial_timestamp, params.count, params.interval)
end

function set_parameters!(params::ForecastParameters, other::ForecastParameters)
    params.horizon = other.horizon
    params.initial_timestamp = other.initial_timestamp
    params.interval = other.interval
    params.count = other.count
    return
end

function set_parameters!(params::ForecastParameters, forecast::Forecast)
    set_parameters!(params, TimeSeriesParameters(forecast).forecast_params)
    return
end

mutable struct TimeSeriesParameters <: InfrastructureSystemsType
    resolution::Dates.Period
    forecast_params::ForecastParameters
end

function TimeSeriesParameters(;
    resolution=UNINITIALIZED_PERIOD,
    forecast_params=ForecastParameters(),
)
    return TimeSeriesParameters(resolution, forecast_params)
end

function TimeSeriesParameters(ts::StaticTimeSeries)
    return TimeSeriesParameters(resolution=get_resolution(ts))
end

function TimeSeriesParameters(ts::Forecast)
    forecast_params = ForecastParameters(
        count=get_count(ts),
        horizon=get_horizon(ts),
        initial_timestamp=get_initial_timestamp(ts),
        interval=get_interval(ts),
    )
    return TimeSeriesParameters(get_resolution(ts), forecast_params)
end

function TimeSeriesParameters(
    initial_timestamp::Dates.DateTime,
    resolution::Dates.Period,
    len::Int,
    horizon::Int,
    interval::Dates.Period,
)
    if interval == Dates.Second(0)
        count = 1
    else
        last_timestamp = initial_timestamp + resolution * (len - 1)
        last_initial_time = last_timestamp - resolution * (horizon - 1)

        # Reduce last_initial_time to the nearest interval if necessary.
        diff =
            Dates.Millisecond(last_initial_time - initial_timestamp) %
            Dates.Millisecond(interval)
        if diff != Dates.Millisecond(0)
            last_initial_time -= diff
        end
        count =
            Dates.Millisecond(last_initial_time - initial_timestamp) /
            Dates.Millisecond(interval) + 1
    end
    fparams = ForecastParameters(
        horizon=horizon,
        initial_timestamp=initial_timestamp,
        interval=interval,
        count=count,
    )
    return TimeSeriesParameters(resolution, fparams)
end

function reset_info!(params::TimeSeriesParameters)
    params.resolution = UNINITIALIZED_PERIOD
    reset_info!(params.forecast_params)
    @info "Reset system time series parameters."
end

function _is_uninitialized(params::TimeSeriesParameters)
    return params.resolution == UNINITIALIZED_PERIOD
end

"""
Return true if `params` match `other` or if one of them is uninitialized.
"""
function check_params_compatibility(
    params::TimeSeriesParameters,
    other::TimeSeriesParameters,
)
    _is_uninitialized(params) && return true

    if other.resolution != params.resolution
        throw(
            ConflictingInputsError(
                "time series resolution $(other.resolution) does not match system " *
                "resolution $(params.resolution)",
            ),
        )
    end

    # `other` might be for a static time series and not have forecast params
    if !_is_uninitialized(other.forecast_params)
        check_params_compatibility(params.forecast_params, other.forecast_params)
    end
end

function check_add_time_series(params::TimeSeriesParameters, ts::TimeSeriesData)
    check_params_compatibility(params, TimeSeriesParameters(ts))
    _check_time_series_lengths(ts)
    return
end

function set_parameters!(params::TimeSeriesParameters, ts::StaticTimeSeries)
    params.resolution = get_resolution(ts)
    return
end

function set_parameters!(params::TimeSeriesParameters, forecast::Forecast)
    params.resolution = get_resolution(forecast)
    set_parameters!(params.forecast_params, forecast)
    return
end

function set_parameters!(params::TimeSeriesParameters, other::TimeSeriesParameters)
    if _is_uninitialized(params)
        # This is the first time series added.
        params.resolution = other.resolution
    end

    if _is_uninitialized(params.forecast_params)
        set_parameters!(params.forecast_params, other.forecast_params)
    end

    return
end

function _check_time_series_lengths(ts::StaticTimeSeries)
    data = get_data(ts)
    if length(data) < 2
        throw(ArgumentError("data array length must be at least 2: $(length(data))"))
    end
    if length(data) != length(ts)
        throw(ConflictingInputsError("length mismatch: $(length(data)) $(length(ts))"))
    end

    timestamps = TimeSeries.timestamp(data)
    difft = timestamps[2] - timestamps[1]
    if difft != get_resolution(ts)
        throw(ConflictingInputsError("resolution mismatch: $difft $(get_resolution(ts))"))
    end
    return
end

function _check_time_series_lengths(ts::Forecast)
    horizon = get_horizon(ts)
    if horizon < 2
        throw(ArgumentError("horizon must be at least 2: $horizon"))
    end
    for window in iterate_windows(ts)
        if size(window)[1] != horizon
            throw(ConflictingInputsError("length mismatch: $(size(window)[1]) $horizon"))
        end
    end
end

get_forecast_window_count(params::TimeSeriesParameters) = params.forecast_params.count
get_forecast_initial_times(params::TimeSeriesParameters) =
    get_forecast_initial_times(params.forecast_params)
get_forecast_horizon(params::TimeSeriesParameters) = params.forecast_params.horizon
get_forecast_initial_timestamp(params::TimeSeriesParameters) =
    params.forecast_params.initial_timestamp
get_forecast_interval(params::TimeSeriesParameters) = params.forecast_params.interval
get_time_series_resolution(params::TimeSeriesParameters) = params.resolution

function get_forecast_total_period(p::TimeSeriesParameters)
    f = p.forecast_params
    _is_uninitialized(f) && return Dates.Second(0)
    return get_total_period(
        f.initial_timestamp,
        f.count,
        f.interval,
        f.horizon,
        p.resolution,
    )
end


text:
##BASIC UTILITY FUNCTIONS
#SOURCE PERMUTATION
function permute_source_weights(source::Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}, shift_freq::AbstractFloat, PWM_shift_dist::Distribution{Univariate,Continuous})
    dirty=false; source_length=size(source[1],1)
    new_source=deepcopy(source)

    for source_pos in 1:source_length
        if rand() <= shift_freq
            pos_WM = exp.(source[1][source_pos,:]) #leaving logspace, get the wm at that position
            new_source[1][source_pos,:] = log.(wm_shift(pos_WM, PWM_shift_dist)) #accumulate probabilty at a randomly selected base, reassign in logspace and carry on
            !dirty && (dirty=true)
        end
    end

    if !dirty #if no positions were shifted, pick one and shift
        rand_pos=rand(1:source_length)
        pos_WM = exp.(source[1][rand_pos,:])
        new_source[1][rand_pos,:]=log.(wm_shift(pos_WM, PWM_shift_dist))
    end

    return new_source
end

                function wm_shift(pos_WM::AbstractVector{<:AbstractFloat}, PWM_shift_dist::Distribution{Univariate,Continuous})
                    base_to_shift = rand(1:4) #pick a base to accumulate probability
                    permute_sign = rand(-1:2:1)
                    shift_size = rand(PWM_shift_dist)
                    new_wm=zeros(4)

                    for base in 1:4 #ACGT
                        if base == base_to_shift
                            new_wm[base] =
                            clamp(0, #no lower than 0 prob
                            (pos_WM[base]          #selected PWM posn
                            + permute_sign * shift_size), #randomly permuted by size param
                            1) #no higher than prob 1
                        else
                            size_frac = shift_size / 3 #other bases shifted in the opposite direction by 1/3 the shift accumulated at the base to permute
                            new_wm[base] =
                            clamp(0,
                            (pos_WM[base]
                            - permute_sign * size_frac),
                            1)
                        end
                    end
                    new_wm = new_wm ./ sum(new_wm) #renormalise to sum 1 - necessary in case of clamping at 0 or 1
                    !isprobvec(new_wm) && throw(DomainError(new_wm, "Bad weight vector generated in wm_shift!")) #throw assertion exception if the position WM is invalid
                    return new_wm
                end

function permute_source_length(source::Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}, prior::Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},<:Bool}, length_limits::UnitRange{<:Integer}, permute_range::UnitRange{<:Integer}=LENGTHPERM_RANGE, uninformative::Dirichlet=Dirichlet([.25,.25,.25,.25]))
    source_PWM, prior_idx = source
    source_length = size(source_PWM,1)

    permute_sign, permute_length = get_length_params(source_length, length_limits, permute_range)

    permute_sign==1 ? permute_pos = rand(1:source_length+1) :
        permute_pos=rand(1:source_length-permute_length)

    if permute_sign == 1 #if we're to add positions to the PWM
        ins_WM=zeros(permute_length,4)
        if prior==false
            for pos in 1:permute_length
                ins_WM[pos,:] = log.(transpose(rand(uninformative)))
            end
        else
            for pos in 1:permute_length
                prior_position=permute_pos+prior_idx-1
                prior_position<1 || prior_position>length(prior) ?
                    ins_WM[pos,:] = log.(transpose(rand(uninformative))) :
                    ins_WM[pos,:] = log.(transpose(rand(prior[prior_position])))
            end
        end
        upstream_source=source_PWM[1:permute_pos-1,:]
        downstream_source=source_PWM[permute_pos:end,:]
        source_PWM=vcat(upstream_source,ins_WM,downstream_source)
        permute_pos==1 && (prior_idx-=permute_length)
    else #if we're to remove positions
        upstream_source=source_PWM[1:permute_pos-1,:]
        downstream_source=source_PWM[permute_pos+permute_length:end,:]
        source_PWM=vcat(upstream_source,downstream_source)
        permute_pos==1 && (prior_idx+=permute_length)
    end

    return (source_PWM, prior_idx) #return a new source
end

                function get_length_params(source_length::Integer, length_limits::UnitRange{<:Integer}, permute_range::UnitRange{<:Integer})
                    extendable = length_limits[end]-source_length
                    contractable =  source_length-length_limits[1]

                    if extendable == 0 && contractable > 0
                        permute_sign=-1
                    elseif contractable == 0 && extendable > 0
                        permute_sign=1
                    else
                        permute_sign = rand(-1:2:1)
                    end

                    permute_sign==1 && extendable<permute_range[end] && (permute_range=permute_range[1]:extendable)
                    permute_sign==-1 && contractable<permute_range[end] && (permute_range=permute_range[1]:contractable)
                    permute_length = rand(permute_range)

                    return permute_sign, permute_length
                end

function erode_source(source::Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer},length_limits::UnitRange{<:Integer},info_thresh)
    pwm,prior_idx=source
    infovec=get_pwm_info(pwm)
    start_idx,end_idx=get_erosion_idxs(infovec, info_thresh, length_limits)

    return new_source=(pwm[start_idx:end_idx,:], prior_idx+start_idx-1)
end

    function get_pwm_info(pwm::AbstractMatrix{<:AbstractFloat}; logsw::Bool=true)
        wml=size(pwm,1)
        infovec=zeros(wml)
        for pos in 1:wml
            logsw ? wvec=deepcopy(exp.(pwm[pos,:])) : wvec=deepcopy(pwm[pos,:])
            !isprobvec(wvec) && throw(DomainError(wvec, "Bad wvec in get_pwm_info -Original sources must be in logspace!!"))
            wvec.+=10^-99
            infscore = (2.0 + sum([x*log(2,x) for x in wvec]))
            infovec[pos]=infscore
        end
        return infovec
    end

    function get_erosion_idxs(infovec::AbstractVector{<:AbstractFloat}, info_thresh::AbstractFloat, length_limits::UnitRange{<:Integer})
        srcl=length(infovec)
        contractable =  srcl-length_limits[1]
        contractable <=0 && throw(DomainError(contractable, "erode_source passed a source at its lower length limit!"))
        centeridx=findmax(infovec)[2]

        start_idx=findprev(info-><(info,info_thresh),infovec,centeridx)
        start_idx===nothing ? (start_idx=1) : (start_idx+=1)
        end_idx=findnext(info-><(info, info_thresh),infovec,centeridx)
        end_idx===nothing ? (end_idx=srcl) : (end_idx-=1)

        pos_to_erode=srcl-(end_idx-start_idx)
        if pos_to_erode > contractable
            pos_to_restore = pos_to_erode-contractable
            while pos_to_restore>0
                end_die=rand()
                if end_die <= .5
                    start_idx>1 && (pos_to_restore-=1; start_idx-=1)
                else
                    end_idx<srcl && (pos_to_restore-=1; end_idx+=1)
                end
            end
        end

        return start_idx, end_idx
    end

#MIX MATRIX FUNCTIONS
function mixvec_decorrelate(mix::BitVector, moves::Integer)
    new_mix=deepcopy(mix)
    idxs_to_flip=rand(1:length(mix), moves)
    new_mix[idxs_to_flip] .= .!mix[idxs_to_flip]
    return new_mix
end

function mix_matrix_decorrelate(mix::BitMatrix, moves::Integer)
    clean=Vector{Bool}(trues(size(mix,1)))
    new_mix=deepcopy(mix)
    indices_to_flip = rand(CartesianIndices(mix), moves)
    new_mix[indices_to_flip] .= .!mix[indices_to_flip]
    clean[unique([idx[1] for idx in indices_to_flip])] .= false #mark all obs that had flipped indices dirty
    return new_mix, clean
end

# function most_dissimilar(mix1, mix2)
#     S1=size(mix1,2);S2=size(mix2,2)
#     dist_mat=zeros(S1,S2)
#     for s1 in 1:S1, s2 in 1:S2
#         dist_mat[s1,s2]=sum(mix1[:,s1].==mix2[:,s2])
#     end
#     scores=vec(sum(dist_mat,dims=1))
#     return findmin(scores)[2]
# end

function most_dissimilar(src_mixvec, target_mixmat)
    src_sim = [sum(src_mixvec.==target_mixmat[:,s]) for s in 1:size(target_mixmat,2)] #compose array of elementwise equality comparisons between mixvectors and sum to score
    merge_s=findmin(src_sim)[2] #source from merger model will be the one with the highest equality comparison score
end

function most_similar(src_mixvec, target_mixmat)
    src_sim = [sum(src_mixvec.==target_mixmat[:,s]) for s in 1:size(target_mixmat,2)] #compose array of elementwise equality comparisons between mixvectors and sum to score
    merge_s=findmax(src_sim)[2] #source from merger model will be the one with the highest equality comparison score
end


text:
using Revise

using JDF, CSV, DataFrames, Blosc, Base.GC, BenchmarkTools, Serialization

# use 12 threads
Blosc.set_num_threads(6)
# GC.gc()
@time a = CSV.read(
    "C:/data/Performance_All/Performance_2010Q3.txt",
    delim = '|',
    header = false
);
GC.gc()

@time metadatas = savejdf("c:/data/large.jdf", a);

using Revise, JDF, DataFrames, Serialization
@time a2 = loadjdf("c:/data/large.jdf"; verbose);
GC.gc()

@time type_compress!(a2);

io = iow()
compress_then_write(a2.Column2, io)
close(io)

@time psavejdf("c:/data/largec.jdf", a2);

using Revise, JDF, DataFrames, Serialization
@time a2 = ploadjdf("c:/data/largec.jdf");
GC.gc()

@time type_compress!(a2, verbose = true)

@time psavejdf("c:/data/largec2.jdf", a2);
@time a2 = loadjdf("c:/data/largec.jdf");

using Revise, JDF, DataFrames, Serialization
@time a2 = ploadjdf("c:/data/largec.jdf");
GC.gc()

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))


text:
abstract type MultivariateFun{T,N} <: Function end
const BivariateFun{T} = MultivariateFun{T,2}

export grad, lap, curl

#implements coefficients/values/evaluate
space(f::MultivariateFun{T,N}) where {T,N}=mapreduce(k->space(f,k),⊗,1:N)
domain(f::MultivariateFun{T,N}) where {T,N}=mapreduce(k->domain(f,k),×,1:N)

domain(f::MultivariateFun,k::Integer)=domain(space(f,k))

differentiate(u::BivariateFun,i::Integer,j::Integer) =
    j==0 ? u : differentiate(differentiate(u,i),i,j-1)
grad(u::BivariateFun) = [differentiate(u,1),differentiate(u,2)]
lap(u::BivariateFun) = differentiate(u,1,2)+differentiate(u,2,2)
Base.div(u::AbstractVector{B}) where {B<:BivariateFun} =
    differentiate(u[1],1)+differentiate(u[2],2)
curl(u::AbstractVector{B}) where {B<:BivariateFun} = differentiate(u[2],1)-differentiate(u[1],2)

Base.chop(f::MultivariateFun) = chop(f,10eps())
cfstype(::MultivariateFun{T}) where {T} = T
cfstype(::Type{MultivariateFun{T,N}}) where {T,N} = T
cfstype(::Type{MF}) where {MF<:MultivariateFun} = cfstype(supertype(MF))

include("VectorFun.jl")
include("TensorSpace.jl")
include("LowRankFun.jl")
include("ProductFun.jl")

arglength(f)=length(Base.uncompressed_ast(f.code.def).args[1])

## Convert between Fun and MultivariateFun
# need to chop extra zeros
Fun(f::ProductFun) =
    Fun(space(f),chop!(fromtensor(space(f),coefficients(f)),0))
Fun(f::ProductFun,sp::TensorSpace) = Fun(ProductFun(f,sp))
Fun(f::LowRankFun) = Fun(ProductFun(f))

Fun(f::MultivariateFun,sp::Space) = Fun(Fun(f),sp)

Fun(f,d1::Domain,d2::Domain) = Fun(f,d1*d2)

coefficients(f::BivariateFun,sp::TensorSpace)=coefficients(f,sp[1],sp[2])

points(f::BivariateFun,k...)=points(space(f),size(f,1),size(f,2),k...)

function *(vx::LowRankFun,u0::ProductFun)
    ret=zero(space(u0))
    for k=1:length(vx.A)
        a,b=vx.A[k],vx.B[k]
        ret+=transpose(b*(transpose(a*u0)))
    end
    ret
end

*(a::ProductFun,b::LowRankFun)=b*a
*(a::MultivariateFun,b::MultivariateFun)=LowRankFun(a)*ProductFun(b)

for OP in (:+,:-,:*,:/)
    @eval begin
        $OP(f::Fun,g::MultivariateFun)=$OP(ProductFun(f),g)
        $OP(f::MultivariateFun,g::Fun)=$OP(f,ProductFun(g))
    end
end

Base.sum(f::Fun{TS},k::Integer) where {TS<:TensorSpace}=sum(ProductFun(f),k)
Base.sum(f::Fun{TS}) where {TS<:TensorSpace}=sum(ProductFun(f))

## kron
# TODO: generalize
function Base.kron(f::Fun,g::Fun)
    sp=space(f)⊗space(g)
    it=tensorizer(sp)
    N=ncoefficients(f);M=ncoefficients(g)
    cfs=Array{promote_type(cfstype(f),cfstype(g))}(undef,0)
    for (k,j) in it
        # Tensor product is N x M, so if we are outside
        # the (N+M)th diagonal we have no more entries
        if k+j > N+M
            break
        elseif k ≤ N && j ≤ M
            push!(cfs,f.coefficients[k]*g.coefficients[j])
        else
            push!(cfs,0)
        end
    end
    Fun(sp,cfs)
end
Base.kron(f::Fun,g::Number) = kron(f,Fun(g))
Base.kron(f::Number,g::Fun) = kron(Fun(f),g)


text:
#=
An sorted array of integers was rotated an unknown number of times.

Given such an array, find the index of the element in the array in faster than linear time. If the element doesn't exist in the array, return null.

For example, given the array [13, 18, 25, 2, 8, 10] and the element 8, return 4 (the index of 8 in the array).

You can assume all the integers in the array are unique.
=#

using Test
include("Solutions/problem58_findindex_sortedarray.jl")

@test findindex([13, 18, 25, 2, 8, 10], 8) == 5 #in Julia indices start at 1

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module FastSigmoid

  abstract Posit{N,ES} <: AbstractFloat
  abstract Valid{N,ES} <: AbstractFloat

  #categorize the different types of posit.
  bitstype 8 Posit8{ES} <: Posit{8,ES}
  bitstype 16 Posit16{ES} <: Posit{16,ES}
  bitstype 32 Posit32{ES} <: Posit{32,ES}
  bitstype 64 Posit64{ES} <: Posit{64,ES}

  #categorize the different types of valid.
  #note that each of the valids is a pair of
  #valids start:stop
  bitstype 16 Valid16{ES} <: Valid{16,ES}
  bitstype 32 Valid32{ES} <: Valid{32,ES}
  bitstype 64 Valid64{ES} <: Valid{64,ES}
  bitstype 128 Valid128{ES} <: Valid{128,ES}

  #for now: make these abstract values
  type QFloat <: AbstractFloat; end

  #allow for a NaN error type.
  type NaNError <: Exception; end

  export Posit, Posit8, Posit16, Posit32, Posit64
  export Valid, Valid16, Valid32, Valid64, Valid128
  export NaNError

  #doing things with fastsigmoids
  include("convertposit.jl")
  include("positops.jl")
  include("positprops.jl")

end # module


text:
ipa = Dict()
ipb = Dict()

# Input membership functions
ipa["small"] = TriangularMF(1, 2, 3)
ipa["large"] = TriangularMF(5, 6, 7)
ipb["small"] = TriangularMF(1, 2, 3)
ips = [ipa, ipb]

# Rules
rule1 = Rule(["small", "small"], [1.0, 1.0, 1.0])
rule2 = Rule(["large", "small"], [0, 0, 17.0])
rules = [rule1, rule2]

fis = FISSugeno(ips, rules)

# Input values
ival1 = [2.3, 1.2] # Runs completely on rule 1
ival2 = [6.4, 2.1] # Runs completely on rule 2

@assert eval_fis(fis, ival1) == 4.5
@assert eval_fis(fis, ival2) == 17.0


text:
function solve()
    x, n = [parse(Int, x) for x in split(readline())]
    p = Set(parse(Int, x) for x in split(readline()))

    d = 0
    while true
        for s in [-1, 1]
            a = x + d*s
            if a ∉ p
                return a
            end
        end
        d += 1
    end

end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#=

# prices - k-dimensional vector of allowed price levels
# demands - n x k matrix of demand predictions
# c - required sum of prices
# where k is the number of price levels, n is number of products
def optimal_prices_category(prices, demands, c):

    n, k = np.shape(demands)

    # prepare inputs
    r = np.multiply(np.tile(prices, n), np.array(demands).reshape(1, k*n))
    A = np.array([[
        1 if j >= k*(i) and j < k*(i+1) else 0
        for j in range(k*n)
    ] for i in range(n)])
    A = np.append(A, np.tile(prices, n), axis=0)
    b = [np.append(np.ones(n), c)]

    # solve the linear program
    res = linprog(-r.flatten(),
              A_eq = A,
              b_eq = b,
              bounds=(0, 1))

    np.array(res.x).reshape(n, k)

# test run
optimal_prices_category(
  np.array([[1.00, 1.50, 2.00, 2.50]]),   # allowed price levels
  np.array([
    [ 28, 23, 20, 13],                    # demands for product 1
    [ 30, 22, 16, 12],                    # demands for product 2
    [ 32, 26, 19, 15]                     # demands for product 3
  ]), 5.50)                               # sum of prices for all three products

# output
[[0 0   1 0  ]                            # price vector for product 1 ($2.00)
 [0 1   0 0  ]                            # price vector for product 2 ($1.50)
 [0 0.5 0 0.5]]                           # price vector for product 3 (
1.50
a
n
d
2.50)

=#


text:
# --- Decl ---
decl_newint(id::String, value::Exp; location = nothing) = decl_new(t_int, id, value=value, location=location)

# --- Stmt ---
stmt_newblock(stmts = []; location = nothing) = stmt_new(k_stmt_block, stmts=stmts, location=location)
stmt_newchain(stmts = []; location = nothing) = stmt_new(k_stmt_chain, stmts=stmts, location=location)
stmt_newdecl(decl::Decl; location = nothing) = stmt_new(k_stmt_decl, decl=decl, location=location)
stmt_newexp(exp::Exp; location = nothing) = stmt_new(k_stmt_exp, exp=exp, location=location)
stmt_newreturn(exp::Exp; location = nothing) = stmt_new(k_stmt_return, exp=exp, location=location)
stmt_newifelse(condition::Exp, iftrue, iffalse = nothing; location = nothing) =
        stmt_new(k_stmt_ifelse, exp=condition, ifbody=iftrue, elsebody=iffalse, location=location)
stmt_newloop(init::Union{Stmt, Nothing}, condition::Exp,
             iter::Union{Stmt, Nothing}, body::Stmt; location = nothing) =
        stmt_new(k_stmt_loop, exp=condition, initbody=init, iterbody=iter, loopbody=body, location=location)

# IO
stmt_newprintint(exp; location = nothing) = stmt_new(k_stmt_printint, exp=exp, location=location)
stmt_newprintstr(s; location = nothing) = stmt_new(k_stmt_printstr, exp=s, location=location)
stmt_newprintline(; location = nothing) = stmt_new(k_stmt_printstr, exp="\\n", location=location)
stmt_newprintspace(; location = nothing) = stmt_new(k_stmt_printstr, exp=" ", location=location)
stmt_newscan(exp; location = nothing) = stmt_new(k_stmt_scan, exp=exp, location=location)

# Wrappers
# loop with <id> from <from> to <to>
function stmt_newloopwith(id, from, to, body; location = nothing)
    init = stmt_newdecl(decl_newint(id, from))
    condition = exp_newtest(exp_newid(id), "<=", to)
    iter = exp_newset(exp_newid(id), exp_newadd(exp_newid(id), exp_newint(1)))

    return stmt_newloop(init, condition, stmt_newexp(iter), body, location=location)
end

# Items is an array of Exp / String
function stmt_newprint(items; location = nothing)
    stmts = []
    for (i, item) in enumerate(items)
        push!(stmts, (isa(item, String) ? stmt_newprintstr : stmt_newprintint)(item))
        push!(stmts, (i == length(items) ? stmt_newprintline : stmt_newprintspace)())
    end

    return stmt_newchain(stmts, location=location)
end

function stmt_newletscan(id; location = nothing)
    stmtlet = stmt_newdecl(decl_newint(id, exp_newint(0), location=location), location=location)
    stmtscan = stmt_newscan(exp_newid(id, location=location), location=location)

    return stmt_newchain([stmtlet, stmtscan], location=location)
end

# --- Exp ---
exp_newint(value::Int; location = nothing) = exp_new(k_exp_int, type=t_int, value=value, location=location)
exp_newid(id::String; location = nothing) = exp_new(k_exp_id, id=id, location=location)
exp_newadd(left, right; location = nothing) = exp_new(k_exp_add, left=left, right=right, location=location)
exp_newdiv(left, right; location = nothing) = exp_new(k_exp_div, left=left, right=right, location=location)
exp_newmod(left, right; location = nothing) = exp_new(k_exp_mod, left=left, right=right, location=location)
exp_newneg(val; location = nothing) = exp_new(k_exp_neg, left=val, location=location)
exp_newmul(left, right; location = nothing) = exp_new(k_exp_mul, left=left, right=right, location=location)
exp_newset(left, right; location = nothing) = exp_new(k_exp_set, left=left, right=right, location=location)
exp_newtest(left, operator, right; location = nothing) = exp_new(k_exp_test, left=left, right=right, operator=operator, location=location)
exp_newbool(left, operator, right; location = nothing) = exp_new(k_exp_bool, left=left, right=right, operator=operator, location=location)
exp_newcall(fun_id, args::Array{Exp}; location = nothing) = exp_new(k_exp_call, id=fun_id, args=args, location=location)

# --- DeclType ---
# proc doesn't return
type_newproc(args) = DeclType(k_proc_t, args)
# fn returns an int
type_newfn(args) = DeclType(k_fn_t, args)

t_int = decltype_new(k_int_t)
t_void = decltype_new(k_void_t)

# --- Arg ---
arg_new(id::String) = Arg(id)


text:
module GaussianTest

using Test
using ForneyLab
using ForneyLab: naturalParams, standardDistribution
using LinearAlgebra: Diagonal

@testset "sample" begin
    # Univariate
    @test isa(sample(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=1.0, v=2.0)), Float64)

    # Multivariate
    @test isa(sample(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.2, 2.7], v=[2.0 -0.5; -0.5 1.5])), Vector{Float64})
    @test isa(sample(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.2, 2.7], v=Diagonal([2.0, 1.5]))), Vector{Float64})
end

@testset "prod!" begin
    # Univariate
    @test ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) * ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=2.0)
    @test ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) * ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=2.0)
    @test ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) * ProbabilityDistribution(Univariate, PointMass, m=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)
    @test ProbabilityDistribution(Univariate, PointMass, m=1.0) * ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)

    # Multivariate
    @test ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=diageye(2)) * ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=diageye(2)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=2.0*diageye(2))
    @test ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) * ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=2.0*diageye(2))
    @test ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) * ProbabilityDistribution(Multivariate, PointMass, m=ones(2)) == ProbabilityDistribution(Multivariate, PointMass, m=ones(2))
    @test ProbabilityDistribution(Multivariate, PointMass, m=ones(2)) * ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=diageye(2)) == ProbabilityDistribution(Multivariate, PointMass, m=ones(2))
end

@testset "natural parameters" begin
    # Univariate
    d = ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=5.0)
    η = naturalParams(d)
    s = standardDistribution(Univariate, Gaussian, η=η)
    @test d.params[:xi] == s.params[:xi] # Test conversion consistency
    @test d.params[:w] == s.params[:w]

    x = [-0.1, 2.0, 15.0]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Univariate, Gaussian, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency

    # Multivariate
    d = ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0, 5.0], w=[2.0 0.1; 0.1 3.0])
    η = naturalParams(d)
    s = standardDistribution(Multivariate, Gaussian, η=η)
    @test d.params[:xi] == s.params[:xi] # Test conversion consistency
    @test d.params[:w] == s.params[:w]

    x = [[-0.1, 2.0], [-3.0, 5.0]]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Multivariate, Gaussian, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency
end

end #module


text:
"""
Abstract supertype of properties that can be checked.

Every concrete subtype should provide the following functions:
  - `inout_map_property(::Property, ::AbstractVector{<:AbstractVector{Int}},
                        ::AbstractVector{Int}, ::Int)::Property`
  - `check_property(::LazySet, ::Property)::Bool`
"""
abstract type Property end

"""
    projection_map(partition::AbstractVector{<:AbstractVector{Int}},
                   blocks::AbstractVector{Int}
                  )::Vector{Int}

Create a sorted list of all dimensions in `blocks`.

### Input

- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition

### Output

A sorted list containing all output dimensions.
"""
function projection_map(partition::AbstractVector{<:AbstractVector{Int}},
                        blocks::AbstractVector{Int}
                       )::Vector{Int}
    proj = Vector{Int}()
    for bi in blocks
        for i in partition[bi]
            push!(proj, i)
        end
    end
    return proj
end