text:
"""
    animate(filename::AbstractString,
            plots::Vector{Plot2D};
            rate=10)

    animate(plots::Vector{Plot2D};rate=10)

Make an `.mp4` video from `plots`, with frame rate `rate`

If "filename" is given, the movie file will be stored there.
Otherwise, the movie will be opened.

# Example
```julia
X = cumsum(randn(100000))
Y = cumsum(randn(100000))
plots = [Plot(Path(X[1:t],Y[1:t])) for t=10:10:10000]
animate(plots)
```
"""
function animate(filename::AbstractString,
                 plots::Vector{Plot2D};
                 rate=6,
                 density=400,
                 quality=100,
                 loop=1,
                 format=:gif)
    bb = boundingbox(map(boundingbox,plots))
    e = ceil(Integer,log10(length(plots)))+1
    N = 10^e
    tempdir = mktempdir()
    framenames = []
    ProgressMeter.@showprogress for (i,P) in enumerate(plots)
        ext = format == :gif ? "pdf" : "png"
        framename = "$tempdir/frame$(string(N+i-1)[2:end]).$ext"
        push!(framenames,framename)
        save(framename,P,bbox=bb)
    end
    if format == :mp4
        clipname = "mymovie.mp4"
        run(`ffmpeg -r $rate -i $tempdir/frame%0$(e)d.png -pix_fmt yuv420p $tempdir/$clipname`)
    elseif format == :gif
        clipname = "myanimated.gif"
        run(`convert -delay $rate -density $density -quality $quality -loop $loop $framenames $tempdir/$clipname`)
    else
        error("available formats are :mp4 and :gif")
    end
    if filename ≠ ""
        cp("$tempdir/$clipname",filename;force=true)
    else
        try
            if Sys.isapple()
                run(`open $tempdir/$clipname`)
            elseif Sys.islinux() || Sys.isbsd()
                run(`xdg-open $tempdir/$clipname`)
            elseif Sys.iswindows()
                run(`start $tempdir/$clipname`)
            end
        catch e
            error(string(
            "Failed to open the movie file.",
                "Error: "), sprint(Base.showerror, e))
        end
    end
end

animate(plots::Vector{Plot2D};kwargs...) =
                                animate("",plots;kwargs...)


text:
module QuantumNumbers

using Base.Iterators: Reverse, flatten, product, reverse
using Printf: @printf, @sprintf
using DataStructures: OrderedDict
using Random: MersenneTwister, seed!, shuffle!
using LinearAlgebra: norm
using ..Combinatorics: Combinations
using ...Prerequisites: Float
using ...Prerequisites.NamedVectors: HomoNamedVector

import ...Interfaces: ⊕, ⊗, dimension, expand, permute, decompose, regularize!, regularize

export AbelianNumber, @abeliannumber, SQN, PQN, SPQN, Z2QN, periods
export qncounts, qnindptr, qncompression, qnexpansion, qncontents, qnindices, qnbruteforce, qnmontecarlo
export AbelianNumbers, SQNS, PQNS, SzPQNS, SPQNS, Z2QNS
export ukron, toordereddict

"Abstract type for all concrete quantum numbers for a single basis."
abstract type AbelianNumber <: HomoNamedVector{Float} end
function periods end

"""
    dimension(::Type{<:AbelianNumber}) -> Int
    dimension(::AbelianNumber) -> Int

The dimension of the Hilbert space an `AbelianNumber` represents. Apparently, this is always 1.
"""
dimension(::Type{<:AbelianNumber}) = 1
dimension(::AbelianNumber) = 1

"""
    regularize!(::Type{QN}, array::AbstractVector{<:Real}) where QN<:AbelianNumber -> typeof(array)
    regularize!(::Type{QN}, array::AbstractMatrix{<:Real}) where QN<:AbelianNumber -> typeof(array)

Regularize the elements of an array in place so that it can represent quantum numbers.
"""
function regularize!(::Type{QN}, array::AbstractVector{<:Real}) where QN<:AbelianNumber
    @assert array|>length == QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i, period) in enumerate(QN|>periods)
        (period == Inf) || @inbounds begin
            remainder = array[i] % period
            array[i] = (remainder < 0) ? remainder+period : remainder
        end
    end
    return array
end
function regularize!(::Type{QN}, array::AbstractMatrix{<:Real}) where QN<:AbelianNumber
    @assert size(array, 1) == QN|>length "regularize! error: not consistent shape of input array and $QN."
    for (i, period) in enumerate(QN|>periods)
        (period == Inf) || @inbounds for j in 1:size(array, 2)
            remainder = array[i, j] % period
            array[i, j] = (remainder < 0) ? remainder+period : remainder
        end
    end
    return array
end

"""
    regularize(::Type{QN}, array::Union{AbstractVector{<:Real}, AbstractMatrix{<:Real}}) where QN<:AbelianNumber -> typeof(array)

Regularize the elements of an array and return a copy that can represent quantum numbers.
"""
function regularize(::Type{QN}, array::Union{AbstractVector{<:Real}, AbstractMatrix{<:Real}}) where QN<:AbelianNumber
    result = copy(array)
    regularize!(QN, result)
    return result
end

"""
    @abeliannumber typename fieldnames fieldperiods

Construct a concrete `AbelianNumber` with the type name being `typename`, fieldnames specified by `fieldnames` and periods specified by `fieldperiods`.
"""
macro abeliannumber(typename, fieldnames, fieldperiods)
    typename = Symbol(typename)
    fieldnames = tuple(eval(fieldnames)...)
    fieldperiods = tuple(eval(fieldperiods)...)
    arguments = ntuple(i->Symbol(:v, i), length(fieldnames))
    @assert length(fieldnames) == length(fieldperiods) "abeliannumber error: number of fieldnames($(length(fieldnames))) and fieldperiods($(length(fieldperiods))) not equal."
    @assert all(isa(name, Symbol) for name in fieldnames) "abeliannumber error: all field names should be Symbol."
    @assert all(fieldperiods .> 0) "abeliannumber error: all field periods should be positive."
    if all(fieldperiods .== Inf)
        title = Expr(:call, :($(esc(typename))), (:($arg::Float) for arg in arguments)...)
        body = Expr(:call, :new, arguments...)
    else
        title = Expr(:call, :($(esc(typename))), (:($arg::Float) for arg in arguments)..., Expr(:kw, :(regularize::Bool), :true))
        body = Expr(:call, :new, ((p == Inf) ? arg : :(regularize ? (r = $arg % $p; (r < 0) ? (r + $p) : r) : $arg) for (p, arg) in zip(fieldperiods, arguments))...)
    end
    newtype = Expr(:struct, false, :($(esc(typename))<:AbelianNumber), Expr(:block, (:($field::Float) for field in fieldnames)..., Expr(:(=), title, body)))
    functions = Expr(:block, :(Base.fieldnames(::Type{<:$(esc(typename))}) = $fieldnames), :(periods(::Type{<:$(esc(typename))}) = $fieldperiods))
    return Expr(:block, :(global periods), :(Base.@__doc__($newtype)), functions)
end

"""
    SQN(Sz::Real)

The concrete `AbelianNumber` of a quantum system with spin z-component `Sz` conserved.
"""
@abeliannumber "SQN" (:Sz,) (Inf,)

"""
    PQN(N::Real)

The concrete `AbelianNumber` of a quantum system with particle number `N` conserved.
"""
@abeliannumber "PQN" (:N,) (Inf,)

"""
    SPQN(N::Real, Sz::Real)

The concrete `AbelianNumber` of a quantum system with both particle number `N` and spin z-component `Sz` conserved.
"""
@abeliannumber "SPQN" (:N, :Sz) (Inf, Inf)

"""
    Z2QN(N::Real)

The concrete `AbelianNumber` of a quantum system with a Z₂-like conserved quantity.
"""
@abeliannumber "Z2QN" (:N,) (2,)

abstract type QNProtocol end
struct QNIndptr <: QNProtocol end
struct QNCounts <: QNProtocol end
"""
    qnindptr

Indicate that methods with `AbelianNumbers` use the index pointer of the compressed contents.
"""
const qnindptr = QNIndptr()
"""
    qncounts

Indicate that methods with `AbelianNumbers` use the count number of the compressed contents.
"""
const qncounts = QNCounts()

"""
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, counts::Vector{Int}, ::QNCounts)
    AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, ::QNIndptr)

The whole quantum numbers of the total bases of a Hilbert space.

The default constructors construct an `AbelianNumbers` from a vector of concrete quantum numbers and an vector containing their counts or indptr.
"""
struct AbelianNumbers{QN<:AbelianNumber}
    form::Char
    contents::Vector{QN}
    indptr::Vector{Int}
    function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, counts::Vector{Int}, ::QNCounts)
        @assert form|>uppercase ∈ ('G', 'U', 'C') "AbelianNumbers error: 'form'($form) is not 'G', 'U' or 'C'."
        @assert length(contents) == length(counts) "AbelianNumbers error: dismatch lengths of contents and counts ($(length(contents))!=$length(counts))."
        return new{contents|>eltype}(form|>uppercase, contents, [0, cumsum(counts)...])
    end
    function AbelianNumbers(form::Char, contents::Vector{<:AbelianNumber}, indptr::Vector{Int}, ::QNIndptr)
        @assert form|>uppercase ∈ ('G', 'U', 'C') "AbelianNumbers error: 'form'($form) is not 'G', 'U' or 'C'."
        @assert length(contents)+1 == length(indptr) "AbelianNumbers error: dismatch shapes of contents and indptr ($(length(contents))+1!=$length(indptr))."
        return new{contents|>eltype}(form|>uppercase, contents, indptr)
    end
end

"""
    AbelianNumbers(qn::AbelianNumber, count::Int=1)

Construct an `AbelianNumbers` with one unique quantum number which occurs `count` times.
"""
AbelianNumbers(qn::AbelianNumber, count::Int=1) = AbelianNumbers('C', [qn], [0, count], qnindptr)

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:Integer})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their counts.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:Integer})
    contents = Vector{od|>keytype}(undef, length(od))
    indptr = zeros(Int, length(od)+1)
    for (i, (qn, count)) in enumerate(od)
        @inbounds contents[i] = qn
        @inbounds indptr[i+1] = indptr[i] + count
    end
    AbelianNumbers('U', contents, indptr, qnindptr)
end

"""
    AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:UnitRange{<:Integer}})

Construct an `AbelianNumbers` from an ordered dict containing concrete quantum numbers and their slices.
"""
function AbelianNumbers(od::OrderedDict{<:AbelianNumber, <:UnitRange{<:Integer}})
    contents = Vector{od|>keytype}(undef, length(od))
    indptr = zeros(Int, length(od)+1)
    for (i, (qn, slice)) in enumerate(od)
        @inbounds contents[i] = qn
        @inbounds @assert indptr[i]+1 == slice.start "AbelianNumbers error: slice not consistent."
        @inbounds indptr[i+1] = slice.stop
    end
    AbelianNumbers('U', contents, indptr, qnindptr)
end

"""
    ==(qns1::AbelianNumbers, qns2::AbelianNumbers) -> Bool

Overloaded equivalent operator.
Two `AbelianNumbers`es are equal to each other if and only if both their `contents`es and `indptr`s are elementwise equal to each other.
!!! note
    It is not necessary for two `AbelianNumbers`es to have the same eltype nor the same form to be equal to each other.
"""
Base.:(==)(qns1::AbelianNumbers, qns2::AbelianNumbers) = (qns1.contents == qns2.contents) && (qns1.indptr == qns2.indptr)

"""
    isequal(qns1::AbelianNumbers, qns2::AbelianNumbers) -> Bool

Overloaded equivalent operator.
Two `AbelianNumbers`es are equal to each other if and only if both their `contents`es and `indptr`s are elementwise equal to each other.
!!! note
    It is not necessary for two `AbelianNumbers`es to have the same eltype nor the same form to be equal to each other.
"""
Base.isequal(qns1::AbelianNumbers, qns2::AbelianNumbers) = isequal(qns1.contents, qns2.contents) && isequal(qns1.indptr, qns2.indptr)

"""
    show(io::IO, qns::AbelianNumbers)

Show an `AbelianNumbers`.
"""
Base.show(io::IO, qns::AbelianNumbers) = @printf io "QNS(%s)" join((@sprintf("%s=>%s", qn, slice) for (qn, slice) in pairs(qns, qnindptr)), ", ")

"""
    string(qns::AbelianNumbers) -> String

Convert an `AbelianNumbers` to string.
"""
Base.string(qns::AbelianNumbers) = @sprintf "QNS(%s, %s)" qns|>length qns|>dimension

"""
    length(qns::AbelianNumbers) -> Int

Get the number of unduplicate qunatum numbers in the `AbelianNumbers`.
"""
Base.length(qns::AbelianNumbers) = length(qns.contents)

"""
    eltype(::Type{<:AbelianNumbers{QN}}) where {QN<:AbelianNumber}
    eltype(qns::AbelianNumbers)

Get the type of the concrete `AbelianNumber` contained in an `AbelianNumbers`.
"""
Base.eltype(::Type{<:AbelianNumbers{QN}}) where {QN<:AbelianNumber} = QN
Base.eltype(qns::AbelianNumbers) = qns |> typeof |> eltype

"""
    getindex(qns::AbelianNumbers, index::Int) -> eltype(qns)
    getindex(qns::AbelianNumbers, slice::UnitRange{Int}) -> AbelianNumbers
    getindex(qns::AbelianNumbers, indices::Vector{Int}) -> AbelianNumbers

Overloaded `[]` operator.
!!! note
    1. For an `AbelianNumbers`, all these `getindex` functions act on its `contents`, i.e. its compressed data, but not on its expansion, i.e. the uncompressed data. This definition is consistent with the [`length`](@ref) function.
    2. When the index is an integer, the result is an `AbelianNumber`, while when the index is a unit range or a vector of intgers, the result is an `AbelianNumbers`. The logic is quite reasonable because such behaviors are much alike to those of a vector container.
"""
Base.getindex(qns::AbelianNumbers, index::Int) = qns.contents[index]
function Base.getindex(qns::AbelianNumbers, slice::UnitRange{Int})
    contents = qns.contents[slice]
    indptr = qns.indptr[slice.start:slice.stop+1]
    indptr .= indptr .- qns.indptr[slice.start]
    return AbelianNumbers(qns.form, contents, indptr, qnindptr)
end
function Base.getindex(qns::AbelianNumbers, indices::Vector{Int})
    contents = Vector{qns|>eltype}(undef, length(indices))
    indptr = zeros(Int, length(indices)+1)
    for (i, index) in enumerate(indices)
        contents[i] = qns.contents[index]
        indptr[i+1] = indptr[i] + qns.indptr[index+1] - qns.indptr[index]
    end
    return AbelianNumbers('U', contents, indptr, qnindptr)
end

"""
    iterate(qns::AbelianNumbers, state::Int=1)
    iterate(rv::Iterators.Reverse{<:AbelianNumbers}, state::Int=length(rv.itr, false))

Iterate or reversely iterate over the concrete `AbelianNumber`s contained in an `AbelianNumbers`.
"""
Base.iterate(qns::AbelianNumbers, state::Int=1) = (state > length(qns)) ? nothing : (@inbounds(qns.contents[state]), state+1)
Base.iterate(rv::Reverse{<:AbelianNumbers}, state::Int=length(rv.itr)) = (state < 1) ? nothing : (@inbounds(rv.itr.contents[state]), state-1)

"""
    keys(qns::AbelianNumbers) -> Vector{qns|>eltype}

Iterate over the concrete `AbelianNumber`s contained in an `AbelianNumbers`.
"""
Base.keys(qns::AbelianNumbers) = qns.contents

"""
    values(qns::AbelianNumbers, ::QNIndptr)
    values(qns::AbelianNumbers, ::QNCounts)

Iterate over the slices/counts of the `AbelianNumbers`.
"""
@views Base.values(qns::AbelianNumbers, ::QNIndptr) = ((start+1):stop for (start, stop) in zip(qns.indptr[1:end-1], qns.indptr[2:end]))
@views Base.values(qns::AbelianNumbers, ::QNCounts) = (stop-start for (start, stop) in zip(qns.indptr[1:end-1], qns.indptr[2:end]))

"""
    pairs(qns::AbelianNumbers, choice::Union{QNIndptr, QNCounts})

Iterate over the `AbelianNumber=>slice` or `AbelianNumber=>count` pairs.
"""
Base.pairs(qns::AbelianNumbers, choice::Union{QNIndptr, QNCounts}) = Base.Generator(=>, keys(qns), values(qns, choice))

"""
    toordereddict(qns::AbelianNumbers, ::QNIndptr) -> OrderedDict{qns|>eltype, UnitRange{Int}}
    toordereddict(qns::AbelianNumbers, ::QNCounts) -> OrderedDict{qns|>eltype, Int}

Convert an `AbelianNumbers` to an ordered dict.
"""
function toordereddict(qns::AbelianNumbers, ::QNIndptr)
    @assert qns.form != 'G' "toordereddict error: input `AbelianNumbers` cannot be `G` formed."
    result = OrderedDict{qns|>eltype, UnitRange{Int}}()
    for i = 1:length(qns)
        @inbounds result[qns.contents[i]] = qns.indptr[i]+1:qns.indptr[i+1]
    end
    return result
end
function toordereddict(qns::AbelianNumbers, ::QNCounts)
    @assert qns.form != 'G' "toordereddict error: input `AbelianNumbers` cannot be `G` formed."
    result = OrderedDict{qns|>eltype, Int}()
    for i = 1:length(qns)
        @inbounds result[qns.contents[i]] = qns.indptr[i+1] - qns.indptr[i]
    end
    return result
end

"""
    dimension(qns::AbelianNumbers) -> Int

The dimension of the Hilbert space an `AbelianNumbers` represents.
"""
dimension(qns::AbelianNumbers) = @inbounds(qns.indptr[end])

"""
    sort(qns::AbelianNumbers) -> Tuple{AbelianNumbers, Vector{Int}}

Sort the quantum numbers of an `AbelianNumber`, return the sorted `AbelianNumber` and the permutation array that sorts the expansion of the original `AbelianNumbers`.
"""
function Base.sort(qns::AbelianNumbers)
    ctpts = sortperm(qns.contents, alg=Base.Sort.QuickSort)
    masks = Vector{Bool}(undef, length(qns.contents))
    masks[1] = true
    unduplicate = 1
    for i = 2:length(ctpts)
        @inbounds masks[i] = qns.contents[ctpts[i]] ≠ qns.contents[ctpts[i-1]]
        @inbounds masks[i] && (unduplicate += 1)
    end
    contents = Vector{qns|>eltype}(undef, unduplicate)
    indptr = zeros(Int, unduplicate+1)
    permutation = Vector{Int}(undef, dimension(qns))
    qncount, ptcount = 0, 0
    for (mask, index) in zip(masks, ctpts)
        @inbounds mask && (qncount += 1; contents[qncount] = qns.contents[index]; indptr[qncount+1] = indptr[qncount])
        @inbounds indptr[qncount+1] += qns.indptr[index+1] - qns.indptr[index]
        for (i, p) in enumerate(@inbounds(qns.indptr[index]+1:qns.indptr[index+1]))
            @inbounds permutation[ptcount+i] = p
        end
        @inbounds ptcount += qns.indptr[index+1] - qns.indptr[index]
    end
    return AbelianNumbers('C', contents, indptr, qnindptr), permutation
end

struct QNCompression <: QNProtocol end
struct QNExpansion <: QNProtocol end
"""
    qncompression

Indicate that [`findall`](@ref) and [`permute`](@ref) use the compressed contents.
"""
const qncompression = QNCompression()
"""
    qnexpansion

Indicate that [`findall`](@ref) and [`permute`](@ref) use the expanded contents.
"""
const qnexpansion = QNExpansion()

"""
    findall(target::QN, qns::AbelianNumbers{QN}, ::QNCompression) where QN<:AbelianNumber -> Vector{Int}
    findall(target::QN, qns::AbelianNumbers{QN}, ::QNExpansion) where QN<:AbelianNumber -> Vector{Int}
    findall(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}, ::QNCompression) where {N, QN<:AbelianNumber} -> Vector{Int}
    findall(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}, ::QNExpansion) where {N, QN<:AbelianNumber} -> Vector{Int}

Find all the indices of the target quantum numbers in the contents ([`qncompression`](@ref) case) or the expansion ([`qnexpansion`](@ref) case) of an `AbelianNumbers`.
"""
Base.findall(target::QN, qns::AbelianNumbers{QN}, ::QNCompression) where QN<:AbelianNumber = findall((target,), qns, qncompression)
Base.findall(target::QN, qns::AbelianNumbers{QN}, ::QNExpansion) where QN<:AbelianNumber = findall((target,), qns, qnexpansion)
function Base.findall(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}, ::QNCompression) where {N, QN<:AbelianNumber}
    result = Int[]
    if qns.form == 'C'
        for qn in targets
            range = searchsorted(qns.contents, qn)
            (range.start <= range.stop) && push!(result, range.start)
        end
    else
        for (i, qn) in enumerate(qns)
            (qn ∈ targets) && push!(result, i)
        end
    end
    result
end
function Base.findall(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}, ::QNExpansion) where {N, QN<:AbelianNumber}
    result = Int[]
    for index in findall(targets, qns, qncompression)
        @inbounds append!(result, qns.indptr[index]+1:qns.indptr[index+1])
    end
    result
end

"""
    filter(target::QN, qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    filter(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}) where {N, QN<:AbelianNumber} -> AbelianNumbers{QN}

Find a subset of an `AbelianNumbers` by picking out the quantum numbers in targets.
"""
Base.filter(target::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = qns[findall(target, qns, qncompression)]
Base.filter(targets::NTuple{N, QN}, qns::AbelianNumbers{QN}) where {N, QN<:AbelianNumber} = qns[findall(targets, qns, qncompression)]

"""
    permute(qns::AbelianNumbers, permutation::Vector{Int}, ::QNCompression) -> AbelianNumbers
    permute(qns::AbelianNumbers, permutation::Vector{Int}, ::QNExpansion) -> AbelianNumbers

Reorder the quantum numbers contained in an `AbelianNumbers` with a permutation and return the new one.

For [`qncompression`](@ref) case, the permutation is for the compressed contents of the original `AbelianNumbers` while for [`qnexpansion`](@ref) case, the permutation is for the expanded contents of the original `AbelianNumbers`.
"""
permute(qns::AbelianNumbers, permutation::Vector{Int}, ::QNCompression) = qns[permutation]
function permute(qns::AbelianNumbers, permutation::Vector{Int}, ::QNExpansion)
    contents = Vector{qns|>eltype}(undef, length(permutation))
    indptr = zeros(Int, length(permutation)+1)
    expansion = expand(qns, qnindices)
    for (i, p) in enumerate(permutation)
        contents[i] = qns.contents[expansion[p]]
        indptr[i+1] = indptr[i] + 1
    end
    return AbelianNumbers('G', contents, indptr, qnindptr)
end

struct QNContents <: QNProtocol end
struct QNIndices <: QNProtocol end
"""
    qncontents

Indicate that [`expand`](@ref) uses the compressed/expanded contents.
"""
const qncontents = QNContents()
"""
    qnindices

Indicate that [`expand`](@ref) uses the indices of the compressed/expanded contents.
"""
const qnindices = QNIndices()

"""
    expand(qns::AbelianNumbers, ::QNContents) -> Vector{qns|>eltype}
    expand(qns::AbelianNumbers, ::QNIndices) -> Vector{Int}

Expand the contents ([`qncontents`](@ref) case) or indices ([`qnindices`](@ref) case) of an `AbelianNumbers` to the uncompressed form.
"""
function expand(qns::AbelianNumbers, ::QNContents)
    result = Vector{qns|>eltype}(undef, dimension(qns))
    for i = 1:length(qns)
        for j = @inbounds(qns.indptr[i]+1:qns.indptr[i+1])
            @inbounds result[j] = qns.contents[i]
        end
    end
    return result
end
function expand(qns::AbelianNumbers, ::QNIndices)
    result = Vector{Int}(undef, dimension(qns))
    for i = 1:length(qns)
        for j = @inbounds(qns.indptr[i]+1:qns.indptr[i+1])
            @inbounds result[j] = i
        end
    end
    return result
end

"""
    +(qn::AbelianNumber) -> typeof(qn)
    +(qn::QN, qns::QN...) where QN<:AbelianNumber -> QN
    +(qns::AbelianNumbers) -> AbelianNumbers
    +(qn::QN, qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    +(qns::AbelianNumbers{QN}, qn::QN) where QN<:AbelianNumber -> AbelianNumbers{QN}

Overloaded `+` operator for `AbelianNumber` and `AbelianNumbers`.
!!! note
    1. The addition between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `+` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` for direct sum and `⊗` for direct product.
    3. To ensure type stability, two `AbelianNumber` can be added together if and only if they are of the same type.
    4. Similarly, an `AbelianNumber` and an `AbelianNumbers` can be added together if and only if the former's type is the same with the latter's eltype.
"""
Base.:+(qn::AbelianNumber) = qn
Base.:+(qn::QN, qns::QN...) where {QN<:AbelianNumber} = map(+, qn, qns...)
Base.:+(qns::AbelianNumbers) = qns
Base.:+(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = qns + qn
Base.:+(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber} = AbelianNumbers((qns.form == 'G') ? 'G' : 'U', [iqn+qn for iqn in qns], qns.indptr, qnindptr)

"""
    -(qn::AbelianNumber) -> typeof(qn)
    -(qn1::QN, qn2::QN) where QN<:AbelianNumber -> QN
    -(qns::AbelianNumbers) -> AbelianNumbers
    -(qn::QN, qns::AbelianNumbers{QN}) where QN<:AbelianNumber -> AbelianNumbers{QN}
    -(qns::AbelianNumbers{QN}, qn::QN) where QN<:AbelianNumber -> AbelianNumbers{QN}

Overloaded `-` operator for `AbelianNumber` and `AbelianNumbers`.
!!! note
    1. The subtraction between an `AbelianNumbers` and an `AbelianNumber` is just a global shift of the contents of the `AbelianNumbers` by the `AbelianNumber`, therefore, the result is an `AbelianNumbers`.
    2. `-` cannot be used between two `AbelianNumbers` because the result is ambiguous. Instead, use `⊕` with signs for direct sum and `⊗` with signs for direct product.
    3. To ensure type stability, an `AbelianNumber` can be subtracted by another `AbelianNumber` if and only if they are of the same type.
    4. Similarly, an `AbelianNumber` can be subtracted by an `AbelianNumbers` or vice versa if and only if the former's type is the same with the latter's eltype.
"""
Base.:-(qn::AbelianNumber) = map(-, qn)
Base.:-(qn1::QN, qn2::QN) where {QN<:AbelianNumber} = map(-, qn1, qn2)
Base.:-(qns::AbelianNumbers) = AbelianNumbers((qns.form == 'G') ? 'G' : 'U', -qns.contents, qns.indptr, qnindptr)
Base.:-(qn::QN, qns::AbelianNumbers{QN}) where {QN<:AbelianNumber} = AbelianNumbers((qns.form == 'G') ? 'G' : 'U', [qn-iqn for iqn in qns], qns.indptr, qnindptr)
Base.:-(qns::AbelianNumbers{QN}, qn::QN) where {QN<:AbelianNumber} = AbelianNumbers((qns.form == 'G') ? 'G' : 'U', [iqn-qn for iqn in qns], qns.indptr, qnindptr)

"""
    *(qn::AbelianNumber, factor::Integer) -> typeof(qn)
    *(factor::Integer, qn::AbelianNumber) -> typeof(qn)
    *(qns::AbelianNumbers, factor::Integer) -> AbelianNumbers
    *(factor::Integer, qns::AbelianNumbers) -> AbelianNumbers

Overloaded `*` operator for the multiplication between an integer and an `AbelianNumber` or an `AbelianNumbers`.
"""
@generated function Base.:*(qn::AbelianNumber, factor::Integer)
    exprs = Expr[:(getfield(qn, $i)*factor) for i = 1:(qn|>fieldnames|>length)]
    return :(typeof(qn)($(exprs...)))
end
Base.:*(factor::Integer, qn::AbelianNumber) = qn * factor
Base.:*(qns::AbelianNumbers, factor::Integer) = AbelianNumbers((qns.form == 'G') ? 'G' : 'U', [qn*factor for qn in qns.contents], qns.indptr, qnindptr)
Base.:*(factor::Integer, qns::AbelianNumbers) = qns * factor

"""
    ^(qn::AbelianNumber, factor::Integer) -> typeof(qn)
    ^(qns::AbelianNumbers, factor::Integer) -> AbelianNumbers

Overloaded `^` operator for `AbelianNumber` and `AbelianNumbers`. This operation translates into the direct product of `factor` copies of `qn` or `qns`.
"""
Base.:^(qn::AbelianNumber, factor::Integer) = kron(NTuple{factor, qn|>typeof}(qn for i = 1:factor)...)
Base.:^(qns::AbelianNumbers, factor::Integer) = kron(NTuple{factor, qns|>typeof}(qns for i = 1:factor)...)

"""
    ⊕(qns::AbelianNumber...) -> AbelianNumbers{qns|>eltype}
    ⊕(qnses::AbelianNumbers...) -> qnses|>eltype

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
"""
⊕(qns::AbelianNumber...) = union(qns...)
⊕(qnses::AbelianNumbers...) = union(qnses...)

"""
    ⊗(qns::AbelianNumber...) -> eltype(qns)
    ⊗(qnses::AbelianNumbers...) -> eltype(qnses)

Get the direct product of some `AbelianNumber`s or `AbelianNumbers`es.
"""
⊗(qns::AbelianNumber...) = kron(qns...)
⊗(qnses::AbelianNumbers...) = kron(qnses...)

"""
    union(qns::Vararg{<:AbelianNumber, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where N -> AbelianNumbers
    union(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N, QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct sum of some `AbelianNumber`s or `AbelianNumbers`es.
!!! note
    1. Physically, the direct sum of a couple of `AbelianNumber`s or `AbelianNumbers`es is defined by the direct sum of the bases of the Hilbert spaces they represent. Therefore, the input `AbelianNumber`s or `AbelianNumbers`es must be homogenous. Inhomogenous 'AbelianNumber's must be direct producted first to ensure homogenity before the direct sum.
    2. Apparently, the dimension of the result equals the summation of those of the inputs, which means, even for `AbelianNumber`s, the result will be naturally an `AbelianNumbers` because the dimension of the result is larger than 1.
    3. Signs of `AbelianNumber`s or `AbelianNumbers`es can be provided when getting their direct sums.
"""
function Base.union(qns::Vararg{<:AbelianNumber, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where N
    AbelianNumbers('G', [(sign == 1) ? qn : -qn for (sign, qn) in zip(signs, qns)], fill(1, qns|>length), qncounts)
end
function Base.union(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N, QN<:AbelianNumber}
    lengths = NTuple{N, Int}(length(qns) for qns in qnses)
    contents = Vector{QN}(undef, sum(lengths))
    indptr = zeros(Int, sum(lengths)+1)
    count = 0
    for (sign, qns, length) in zip(signs, qnses, lengths)
        for i = 1:length
            @inbounds contents[count+i] = (sign == 1) ? qns.contents[i] : -qns.contents[i]
            @inbounds indptr[count+i+1] = indptr[count+i] + qns.indptr[i+1] - qns.indptr[i]
        end
        count += length
    end
    AbelianNumbers('G', contents, indptr, qnindptr)
end

"""
    kron(::Type{QN}, qn1::AbelianNumber, qn2::AbelianNumber) where QN<:AbelianNumber -> QN
    kron(qns::Vararg{<:AbelianNumber, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where N -> eltype(qns)
    kron(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N, QN<:AbelianNumber} -> AbelianNumbers{QN}

Get the direct product of some `AbelianNumber`s or `AbelianNumbers`es.
!!! note
    1. Physically, the direct product of a couple of `AbelianNumber`s or `AbelianNumbers`es are defined by the direct product of the bases of the Hilbert spaces they represent. Therefore, `AbelianNumbers` with differenct types or `AbelianNumbers`es with differenct eltypes are allowed to be direct producted in principle. However, for simplicity, we only implement a method which handle the situation of two `AbelianNumber`s with differenct types. The type of the result should be provided as the first parameter. Note that in this situation, the `fieldnames` and `periods` of the result type must be exactly equal to the flattened fieldnames and periods of the two input `AbelianNumber`s, which means, even the order of the input `AbelianNumber`s matters.
    2. Apparently, the dimension of the result equals the product of those of the inputs. Therefore, the direct product of `AbelianNumber`s is also an `AbelianNumber` since its dimension is still one.
    3. For other situations except the one mentioned in Note.1, the input `AbelianNumber`s or `AbelianNumbers`es must be homogenous. Meanwhile, signs can also be provided for these situations. Note that each quantum number in the contents of the result is obtained by a summation of the corresponding quanum numbers out of the inputs with the correct signs. This is a direct observation of the Abelian nature of our quantum numbers.
"""
function Base.kron(::Type{QN}, qn1::AbelianNumber, qn2::AbelianNumber) where QN<:AbelianNumber
    qnnames, qn1names, qn2names = QN|>fieldnames, qn1|>typeof|>fieldnames, qn2|>typeof|>fieldnames
    qnperiods, qn1periods, qn2periods = QN|>periods, qn1|>typeof|>periods, qn2|>typeof|>periods
    @assert qnnames == NTuple{QN|>length, Symbol}(flatten((qn1names, qn2names))) "kron error: fieldnames not match ($(qnnames), $(qn1names), $(qn2names))."
    @assert qnperiods == NTuple{QN|>length, Float}(flatten((qn1periods, qn2periods))) "kron error: periods not match ($(qnperiods), $(qn1periods), $(qn2periods))."
    QN(qn1..., qn2...)
end
Base.kron(qns::Vararg{<:AbelianNumber, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N} = sum((sign == 1) ? qn : -qn for (sign, qn) in zip(signs, qns))
function Base.kron(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N, QN<:AbelianNumber}
    lengths = NTuple{N, Int}((i < N) ? dimension(qns) : length(qns) for (i, qns) in enumerate(qnses))
    contents = Vector{QN}(undef, prod(lengths))
    indptr = zeros(Int, prod(lengths)+1)
    cache = Vector{QN}(undef, N)
    @inbounds expansions = NTuple{N-1, Vector{Int}}(expand(qnses[i], qnindices) for i = 1:N-1)
    for (i, indices) in enumerate(product(NTuple{N, UnitRange{Int64}}(1:length for length in reverse(lengths))...))
        for (j, (sign, qns, index)) in enumerate(zip(signs, qnses, reverse(indices)))
            @inbounds pos = (j < N) ? expansions[j][index] : index
            @inbounds cache[j] = (sign == 1) ? qns.contents[pos] : -qns.contents[pos]
        end
        @inbounds indptr[i+1] = indptr[i] + qnses[end].indptr[indices[1]+1] - qnses[end].indptr[indices[1]]
        @inbounds contents[i] = sum(cache)
    end
    AbelianNumbers('G', contents, indptr, qnindptr)
end

"""
    ukron(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, N)) where {N, QN<:AbelianNumber} -> AbelianNumbers{QN}, Dict{QN, Dict{NTuple{N, QN}, UnitRange{Int}}}

Unitary Kronecker product of several `AbelianNumbers`es. The product result as well as the records of the product will be returned.
!!! note
    1. All input `AbelianNumbers` must be 'U' formed or 'C' formed.
    2. Since duplicate quantum number are not allowed in 'U' formed and 'C' formed `AbelianNumbers`es, in general, there exists a merge process of duplicate quantum numbers in the product result. Therefore, records are needed to keep track of this process, which will be returned along with the product result. The records are stored in a `Dict{QN, Dict{NTuple{N, QN}, UnitRange{Int}}}` typed dict, in which, for each unduplicate quantum number `qn` in the product result, there exist a record `Dict((qn₁, qn₂, ...)=>start:stop, ...)` telling what quantum numbers `(qn₁, qn₂, ...)` a mereged duplicate `qn` comes from and what slice `start:stop` this merged duplicate corresponds.
"""
function ukron(qnses::Vararg{AbelianNumbers{QN}, N}; signs::NTuple{N, Int}=ntuple(i->1, Val(N))) where {N, QN<:AbelianNumber}
    @assert all((qns.form == 'U') || (qns.form == 'C') for qns in qnses) "ukron error: all input qnses should be 'U' formed or 'C' formed."
    lengths = NTuple{N, Int}(length(qns) for qns in qnses)
    cache = Vector{QN}(undef, N)
    container = OrderedDict{QN, Int}()
    records = Dict{QN, Dict{NTuple{N, QN}, UnitRange{Int}}}()
    for indices in product(NTuple{N, UnitRange{Int64}}(1:length for length in reverse(lengths))...)
        qn, count = QN|>zero, 1
        for (j, (sign, qns, index)) in enumerate(zip(signs, qnses, reverse(indices)))
            @inbounds cache[j] = qns.contents[index]
            @inbounds qn = (sign == 1) ? (qn + cache[j]) : (qn - cache[j])
            @inbounds count = count * (qns.indptr[index+1]-qns.indptr[index])
        end
        container[qn] = get(container, qn, 0) + count
        !haskey(records, qn) && (records[qn] = Dict{NTuple{N, QN}, UnitRange{Int}}())
        records[qn][NTuple{N, QN}(cache)] = (container[qn]-count+1):container[qn]
    end
    contents = QN[qn for qn in keys(container)]
    sort!(contents)
    indptr = zeros(Int, length(container)+1)
    for (i, qn) in enumerate(contents)
        @inbounds indptr[i+1] = indptr[i] + container[qn]
    end
    AbelianNumbers('C', contents, indptr, qnindptr), records
end

struct QNBruteForce <: QNProtocol end
struct QNMonteCarlo <: QNProtocol end
"""
    qnbruteforce

Indicate that [`decompose`](@ref) uses the brute force method.
"""
const qnbruteforce = QNBruteForce()

"""
    qnmontecarlo

Indicate that [`decompose`](@ref) uses the Monte Carlo method.
"""
const qnmontecarlo = QNMonteCarlo()

"""
    decompose(qnses::NTuple{N, AbelianNumbers{QN}}, target::QN, signs::NTuple{N, Int}, ::QNBruteForce; nmax::Int=20) where {N, QN<:AbelianNumber} -> Vector{NTuple{N, Int}}
    decompose(qnses::NTuple{N, AbelianNumbers{QN}}, target::QN, signs::NTuple{N, Int}, ::QNMonteCarlo; nmax::Int=20) where {N, QN<:AbelianNumber} -> Vector{NTuple{N, Int}}

Find a couple of decompositions of `target` with respect to `qnses`.
!!! note
    A tuple of integers `(i₁, i₂, ...)` is called a decomposition of a given `target` with respect to the given `qnses` if and only if they satisfy the "decomposition rule":
    ```math
    \\sum_\\text{j} \\text{signs}[\\text{j}]\\times\\text{qnses}[\\text{j}][\\text{i}_{\\text{j}}]==\\text{target}
    ```
    This equation is in fact a kind of a set of restricted [linear Diophantine equations](https://en.wikipedia.org/wiki/Diophantine_equation#Linear_Diophantine_equations). Indeed, our quantum numbers are always discrete Abelian ones and all instances of a concrete `AbelianNumber` forms a [module](https://en.wikipedia.org/wiki/Module_(mathematics)) over the [ring](https://en.wikipedia.org/wiki/Ring_(mathematics)) of integers. Therefore, each quantum number can be represented as a integral multiple of the unit element of the Abelian module, which results in the final reduction of the above equation to a set of linear Diophantine equations. Then finding a decomposition is equivalent to find a solution of the reduced linear Diophantine equations, with the restriction that the quantum numbers constructed from the solution should be in the corresponding `qnses`. Here we provide two methods to find such decompositions, one is by brute force ([`qnbruteforce`](@ref) case), and the other is by Monte Carlo simultatioins ([`qnmontecarlo`](@ref) case).
"""
function decompose(qnses::NTuple{N, AbelianNumbers{QN}}, target::QN, signs::NTuple{N, Int}, ::QNBruteForce; nmax::Int=20) where {N, QN<:AbelianNumber}
    result = Set{NTuple{N, Int}}()
    cache = Vector{Int}(undef, N)
    dimensions = NTuple{N, Int}(dimension(qns) for qns in reverse(qnses))
    indices = findall(target, kron(qnses..., signs=signs), qnexpansion)
    (nmax < length(indices)) && (shuffle!(MersenneTwister(), indices); indices = @views indices[1:nmax])
    for index in indices
        for (i, dimension) in enumerate(dimensions)
            @inbounds cache[end+1-i] = (index-1)%dimension + 1
            @inbounds index = (index-1)÷dimension + 1
        end
        push!(result, NTuple{N, Int}(cache))
    end
    return collect(NTuple{N, Int}, result)
end
function decompose(qnses::NTuple{N, AbelianNumbers{QN}}, target::QN, signs::NTuple{N, Int}, ::QNMonteCarlo; nmax::Int=20) where {N, QN<:AbelianNumber}
    seed!()
    result = Set{NTuple{N, Int}}()
    expansions = Vector{QN}[expand(qns, qncontents) for qns in qnses]
    dimensions = NTuple{N, Int}(dimension(qns) for qns in qnses)
    diff = indices->norm(sum((sign == 1) ? expansion[index] : -expansion[index] for (sign, expansion, index) in zip(signs, expansions, indices))-target)
    count = 1
    while true
        oldindices, newindices = Int[rand(1:dimension) for dimension in dimensions], ones(Int, N)
        olddiff, newdiff = diff(oldindices), diff(newindices)
        while newdiff > 0
            pos = rand(1:N)
            index = rand(1:dimensions[pos]-1)
            @inbounds newindices[pos] = (index < newindices[pos]) ? index : index+1
            newdiff = diff(newindices)
            if (newdiff <= olddiff) || (exp(olddiff-newdiff) > rand())
                oldindices[:] = newindices
                olddiff = newdiff
            end
        newindices[:] = oldindices
        end
        count = count + 1
        push!(result, NTuple{N, Int}(newindices))
        ((length(result) >= nmax) || (count > nmax*5)) && break
    end
    return collect(NTuple{N, Int}, result)
end

"""
    SQNS(S::Real) -> AbelianNumbers{SQN}

Construct the `AbelianNumbers` of the Hilbert space of a signle spin `S`.
"""
SQNS(S::Real) = AbelianNumbers('C', [SQN(sz) for sz = -S:S], collect(0:Int(2*S+1)), qnindptr)

"""
    PQNS(N::Real) -> AbelianNumbers{PQN}

Construct the `AbelianNumbers` of the Hilbert space of a single-particle state with at most `N` identical particles.
"""
PQNS(N::Real) = AbelianNumbers('C', [PQN(np) for np = 0:N], collect(0:Int(N)+1), qnindptr)

"""
    SzPQNS(Sz::Real) -> AbelianNumbers{SPQN}

Construct the `AbelianNumbers` of the Hilbert space of a single-paritcle state with at most one particle whose spin-z component is `Sz`.
"""
SzPQNS(Sz::Real) = AbelianNumbers('C', [SPQN(0.0, 0.0), SPQN(1.0, Sz)], [0, 1, 2], qnindptr)

"""
    SPQNS(S::Real) -> AbelianNumbers{SPQN}

Construct the `AbelianNumbers` of the Hilbert space of a single site with internal degrees of freedom that can be ascribed to a spin `S`.
"""
function SPQNS(S::Real)
    sqns = [SQN(sz) for sz = -S:S]
    contents = [SPQN(0.0, 0.0)]
    for n = 1:length(sqns)
        pn = PQN(n*1.0)
        for spins in Combinations{n}(sqns)
            push!(contents, kron(SPQN, pn, sum(spins)))
        end
    end
    return sort(AbelianNumbers('G', contents, collect(0:length(contents)), qnindptr))[1]
end

"""
    Z2QNS() -> AbelianNumbers{Z2QN}

Construct the `AbelianNumbers` of a ``Z_2`` Hilbert space.
"""
Z2QNS() = AbelianNumbers('C', [Z2QN(0.0), Z2QN(1.0)], [0, 1, 2], qnindptr)

end #module

ext: jl


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

# D, Δx = 1.0, 0.35

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
# Γ = meshsphere(1.0, 0.08)
@show numcells(Γ)
X, Y = raviartthomas(Γ), buffachristiansen(Γ)

# Δt, Nt = 0.11, 200
Δt, Nt = 0.6, 200
# Δt, Nt = 0.05, 400
T = timebasisshiftedlagrange(Δt, Nt, 2)
δ = timebasisdelta(Δt, Nt)

V = X ⊗ T
W = Y ⊗ δ
# width, delay, scaling = 8.0, 12.0, 1.0
duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = derive(creategaussian(duration, delay, amplitude))

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
# E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)
H = direction × E

@hilbertspace j; @hilbertspace m′
# K, I, N = MWDoubleLayerTDIO(1.0, 1.0, 0), Identity(), NCross()
K = MWDoubleLayerTDIO(1.0, 1.0, 1)
N = BEAST.TemporalDifferentiation(NCross()⊗Identity())

M = 0.5*N + 1.0*K
Z_mfie = assemble(M, W, V, Val{:bandedstorage})
b_mfie = assemble(H, W)
dot_xmfie = marchonintime(inv(Z_mfie[:,:,1]), Z_mfie, b_mfie, Nt)

# Xmfie, Δω, ω0 = fouriertransform(xmfie, Δt, 0.0, 2)
# ω = collect(ω0 .+ (0:Nt-1)*Δω)
# _, i1 = findmin(abs.(ω .- 1.0))
#
# ω1 = ω[i1]
# um = Xmfie[:,i1] / fouriertransform(gaussian)(ω1)


text:
@POMDP_require solve(solver::DeepQLearningSolver, problem::Union{MDP,POMDP}) begin
    P = typeof(problem)
    S = state_type(problem)
    A = action_type(problem)
    @req n_actions(::P)
    @req initial_state(::P, ::AbstractRNG)
    @req convert_o(::Type{Vector{Float64}}, ::S, ::P)
    @req action_index(::P, ::A)

end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using InverseLaplace
using Random
using StatsFuns, SpecialFunctions
using ForwardDiff

struct LaplaceTransformDistribution{T,TAlg}
    f::Function # Laplace transform of the pdf
    c²::T
    alg::TAlg # Algorithm to compute the inverse Laplace transform
    function LaplaceTransformDistribution{T,TAlg}(f::Function,c²::T=0.0,alg::TAlg=BromwichInverseLaplace()) where {T<:Real,TAlg}
        @assert _check_f(f) "The function passed is not valid"# Do series of check on f
        @assert c²>0 "c² has to a be non-negative real"
        new{T,TAlg}(f,c²,alg)
    end
end

function _check_f(f)
    return true
end
_gradf(d::LaplaceTransformDistribution,x::Real) = ForwardDiff.gradient(dist.f,[x])[1]
_gradlogf(d::LaplaceTransformDistribution,x::Real) = ForwardDiff.gradient(log∘dist.f,[x])[1]
_hessianlogf(d::LaplaceTransformDistribution,x::Real) = ForwardDiff.hessian(log∘dist.f,[x])[1]

Distributions.pdf(dist::LaplaceTransformDistribution,x::Real) = apply_f(dist,x)
Distributions.mean(dist::LaplaceTransformDistribution) = _gradf(dist,dist.c²)/dist.f(dist.c²)
Distributions.var(dist::LaplaceTransformDistribution) = _hessianlogf(dist,dist.c²)/dist.f(dist.c²)-mean(dist)^2

struct BromwichInverseLaplace{T}
    A::Int
    l::Int
    m::Int
    n::Int
    expA2l::T
    ijπl::Vector{Complex{T}}
    expijπl::Vector{Complex{T}}
    coefs::Vector{T}
    altern::Vector{Int}
    b::Vector{T}
    function BromwichInverseLaplace(A::Int=19,l::Int=1,m::Int=11,n::Int=28)
        @assert A>0 "A should be a positive integer"
        @assert l>0 "l should be a positive integer"
        @assert m>0 "m should be a positive integer"
        @assert n>0 "n should be a positive integer"
        T = typeof(exp(A/(2*l)))
        ijπl = 1im.*(1:l).*π/l
        new{T}(A,l,m,n,exp(A/(2*l)),ijπl,exp.(ijπl),binomial.(m,0:m)./(2^m),(-1).^(0:(n+m)),zeros(n+m),zeros())
    end
end

function Random.rand(dist::LaplaceTransformDistribution,n::Int)
    nTries = 100
    i = 0
    while i < nTries
        try
            if i > 10
                @warn "This is taking more than 10 tries! Trie $i/$nTries!"
            end
            return laptrans(dist,n=n)
        catch e
            if e isa AssertionError
                @warn "$(e.msg), Trying again"
                i+=1
                continue
            else
                rethrow(e)
            end
        end
    end
    @error "Sampler failed to converge"
end

@inline function apply_f(dist::LaplaceTransformDistribution,t)
    if t==0
        return zero(t)
    else
        invlaplace(dist.f,t,dist.alg)
    end
end

@inline function apply_F(dist::LaplaceTransformDistribution,t)
    #Take care of the spe
    if t == 0
        return zero(t)
    else
        invlaplace(s -> dist.f(s) / s,t,dist.alg)
    end
end

function invlaplace(f::Function,t::Real,A::Int=19,l::Int=1,m::Int=11,n::Int=28)
    # n is the burnin, m is the number of components, A and l aare discretization parameters
    alg = BromwichInverseLaplace(A,l,m,n)
    invlaplace(f,t,alg)
end

function invlaplace(f::Function,t::Real,alg::BromwichInverseLaplace)
    scaled_t = 1/(2*alg.l*t)
    alg.b[1] = f(alg.A*scaled_t) + 2*sum(real.(f.(alg.A*scaled_t .+ alg.ijπl./t).*alg.expijπl))
    @inbounds for k in 1:(alg.n+alg.m)
        alg.b[k+1] = 2*sum(real.(f.(alg.A*scaled_t .+ alg.ijπl./t .+ 1im*π*k/t).*alg.expijπl))
    end
    acoeff = alg.eA2l*scaled_t
    cumsum!(alg.s,acoeff.*alg.b.*alg.altern)
    return dot(view(alg.s,(alg.n+1):end),alg.coeffs)
end

# Laplace transform implemented from "Computational Probability Grassmann 2000"
function laptrans(dist::LaplaceTransformSampler;n::Int=10,jmax::Int=500,kmax::Int=1000,b::Real=2.0,τ::Real=1e-7)
    # Step 1
    global u = sort!(rand(n))
    # Step 2
    global xmax = 0.1
    j = 0
    while apply_F(dist,xmax) < u[end] && j < jmax
        xmax = b*xmax
        j += 1
    end
    @assert j!=jmax "Xmax search failed to converge"
    x = zeros(n)
    # Step 3
    for i in 1:n
        x_L  = i == 1 ? 0.0 : x[i-1]; x_U = xmax; k = 0; t = x_L;
        while abs(apply_F(dist,t)-u[i])>τ && k <kmax
            k += 1
            t = t - (apply_F(dist,t)-u[i])/apply_f(dist,t)
            if t < x_L || t > x_U
                t = 0.5*(x_L+x_U)
            end
            if apply_F(dist,t) <= u[i]
                x_L = t
            else
                x_U = t
            end
        end
        @assert k != kmax "Modified Newton-Raphson procedure failed to converge for i=$i, x[$i]=$t"
        x[i] = t
    end
    # Step 4
    return shuffle!(x)
end

ν=3.0

testfunction(x) = (1+x/ν)^(-0.5*(ν+1))
testfunction(x) = sech(sqrt(x/2))
testfunction(x) = (1+sqrt(3*x))*exp(-sqrt(3*x))
testfunction(x) = (1+sqrt(5*x)+5/3*x)*exp(-sqrt(5*x))
testdist = LaplaceTransformSampler(testfunction,nothing)

laptrans(testdist,n=1)[1]
using Plots
using AugmentedGaussianProcesses
p = histogram(rand(testdist,10000),bins=range(0,3,length=100),normalize=true)|>display
using Distributions
histogram!(rand(Gamma((ν+1)/2,1/ν),10000),alpha=0.5,bins=range(0,10,length=100),normalize=true)
histogram!([AugmentedGaussianProcesses.draw(AugmentedGaussianProcesses.PolyaGammaDist(),1,0) for i in 1:10000],alpha=0.5,bins=range(0,10,length=100),normalize=true)


text:
using Caesar
using Test

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


text:
module AttrGroups

using JSON

_symbol_dict(x) = x
_symbol_dict(d::Associative) =
    Dict{Symbol,Any}([(Symbol(k), _symbol_dict(v)) for (k, v) in d])

function main()

    data = _symbol_dict(JSON.parsefile(joinpath(@__DIR__, "plotschema.json")))

    nms = Set{Symbol}()
    function add_to_names!(d::Associative)
        map(add_to_names!, keys(d))
        map(add_to_names!, values(d))
        nothing
    end
    add_to_names!(s::Symbol) = push!(nms, s)
    add_to_names!(x) = nothing

    add_to_names!(data[:schema][:layout][:layoutAttributes])
    for (_, v) in data[:schema][:traces]
        add_to_names!(v)
    end

    _UNDERSCORE_ATTRS = collect(
        filter(
            x-> contains(string(x), "_") && !startswith(string(x), "_"),
            nms
        )
    )

    _SRC_ATTRS = collect(filter(x -> endswith(string(x), "src"), nms))

    open(joinpath(@__DIR__, "src_attrs.csv"), "w") do f
        writecsv(f, map(string, _SRC_ATTRS))
    end

    open(joinpath(@__DIR__, "underscore_attrs.csv"), "w") do f
        writecsv(f, map(string, _UNDERSCORE_ATTRS))
    end

    _UNDERSCORE_ATTRS, _SRC_ATTRS
end

end  # module


text:
using Base: @propagate_inbounds
using CUDA
using Adapt
using OffsetArrays
using Statistics

using Oceananigans.Architectures
using Oceananigans.Utils
using Oceananigans.Grids: interior_indices, interior_parent_indices

import Base: minimum, maximum, extrema
import Statistics: mean
import Oceananigans: location, instantiated_location
import Oceananigans.Architectures: architecture
import Oceananigans.BoundaryConditions: fill_halo_regions!
import Oceananigans.Grids: interior_x_indices, interior_y_indices, interior_z_indices
import Oceananigans.Grids: total_size, topology, nodes, xnodes, ynodes, znodes, xnode, ynode, znode
import Oceananigans.Utils: datatuple

const ArchOrNothing = Union{AbstractArchitecture, Nothing}
const GridOrNothing = Union{AbstractGrid, Nothing}

"""
    AbstractField{X, Y, Z, A, G, T, N}

Abstract supertype for fields located at `(X, Y, Z)` on architecture `A`
and defined on a grid `G` with eltype `T` and `N` dimensions.
"""
abstract type AbstractField{X, Y, Z, A <: ArchOrNothing, G <: GridOrNothing, T, N} <: AbstractArray{T, N} end

"""
    AbstractDataField{X, Y, Z, A, G, T, N}

Abstract supertype for fields with concrete data in settable underlying arrays,
located at `(X, Y, Z)` on architecture `A` and defined on a grid `G` with eltype `T`
and `N` dimensions.
"""
abstract type AbstractDataField{X, Y, Z, A, G, T, N} <: AbstractField{X, Y, Z, A, G, T, N} end

Base.IndexStyle(::AbstractField) = IndexCartesian()

function validate_field_data(X, Y, Z, data, grid)
    Tx, Ty, Tz = total_size((X, Y, Z), grid)

    if size(data) != (Tx, Ty, Tz)
        e = "Cannot construct field at ($X, $Y, $Z) with size(data)=$(size(data)). " *
            "`data` must have size ($Tx, $Ty, $Tz)."
        throw(ArgumentError(e))
    end

    return nothing
end

# Endpoint for recursive `datatuple` function:
@inline datatuple(obj::AbstractField) = data(obj)

#####
##### Computing AbstractField
#####

# Note: overload compute! for custom fields to produce non-default behavior

"""
    compute!(field)

Computes `field.data`.
"""
compute!(field) = nothing

"""
    compute_at!(field, time)

Computes `field.data` at `time`. Falls back to compute!(field).
"""
compute_at!(field, time) = compute!(field)

mutable struct FieldStatus{T}
    time :: T
end

Adapt.adapt_structure(to, status::FieldStatus) = (time = status.time,)

"""
    conditional_compute!(field, time)

Computes `field.data` if `time != field.status.time`.
"""
function conditional_compute!(field, time)

    if time == zero(time) || time != field.status.time
        compute!(field, time)
        field.status.time = time
    end

    return nothing
end

# This edge case occurs if `fetch_output` is called with `model::Nothing`.
# We do the safe thing here and always compute.
conditional_compute!(field, ::Nothing) = compute!(field, nothing)

"""
    @compute(exprs...)

Call compute! on fields after defining them.
"""
macro compute(def)
    expr = Expr(:block)
    field = def.args[1]
    push!(expr.args, :($(esc(def))))
    push!(expr.args, :(compute!($(esc(field)))))
    return expr
end

#####
##### AbstractField functionality
#####

@inline location(a) = (Nothing, Nothing, Nothing)

"Returns the location `(X, Y, Z)` of an `AbstractField{X, Y, Z}`."
@inline location(::AbstractField{X, Y, Z}) where {X, Y, Z} = (X, Y, Z) # note no instantiation
@inline location(f, i) = location(f)[i]

@inline instantiated_location(::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = (LX(), LY(), LZ())

"Returns the architecture of on which `f` is defined."
architecture(f::AbstractField) = f.architecture

"Returns the topology of a fields' `grid`."
@inline topology(f, args...) = topology(f.grid, args...)

"""
    size(f::AbstractField)

Returns the size of an `AbstractField{X, Y, Z}` located at `X, Y, Z`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
Base.size(f::AbstractField) = size(location(f), f.grid)

"Returns the length of a field's `data`."
@inline Base.length(f::AbstractField) = prod(size(f))

"""
    total_size(field::AbstractField)

Returns a 3-tuple that gives the "total" size of a field including
both interior points and halo points.
"""
total_size(f::AbstractField) = total_size(location(f), f.grid)

Base.fill!(f::AbstractDataField, val) = fill!(parent(f), val)

#####
##### Accessing wrapped arrays
#####

"Returns `f.data` for `f::Field` or `f` for `f::AbstractArray."
@inline data(a) = nothing # fallback
@inline data(f::AbstractDataField) = f.data

@inline cpudata(a) = data(a)

@inline cpudata(f::AbstractField{X, Y, Z, <:AbstractGPUArchitecture}) where {X, Y, Z} =
    offset_data(Array(parent(f)), f.grid, location(f))

"Returns `f.data.parent` for `f::Field`."
@inline Base.parent(f::AbstractField) = parent(data(f))

"Returns a view of `f` that excludes halo points."
@inline interior(f::AbstractDataField{X, Y, Z}) where {X, Y, Z} =
    view(parent(f), interior_parent_indices(X, topology(f, 1), f.grid.Nx, f.grid.Hx),
                    interior_parent_indices(Y, topology(f, 2), f.grid.Ny, f.grid.Hy),
                    interior_parent_indices(Z, topology(f, 3), f.grid.Nz, f.grid.Hz))

@inline interior(f::AbstractField) = f

@inline interior_copy(f::AbstractField{X, Y, Z}) where {X, Y, Z} =
    parent(f)[interior_parent_indices(X, topology(f, 1), f.grid.Nx, f.grid.Hx),
              interior_parent_indices(Y, topology(f, 2), f.grid.Ny, f.grid.Hy),
              interior_parent_indices(Z, topology(f, 3), f.grid.Nz, f.grid.Hz)]

#####
##### getindex
#####

# Don't use axes(f) to checkbounds; use axes(f.data)
Base.checkbounds(f::AbstractDataField, I...) = Base.checkbounds(f.data, I...)

@propagate_inbounds Base.getindex(f::AbstractDataField, inds...) = getindex(f.data, inds...)

# Linear indexing
@propagate_inbounds Base.getindex(f::AbstractDataField, i::Int)  = parent(f)[i]

#####
##### setindex
#####

@propagate_inbounds function Base.setindex!(f::AbstractDataField, val, i, j, k)
    f.data[i, j, k] = val
    return f
end

# Linear indexing
@propagate_inbounds function Base.setindex!(f::AbstractDataField, val, i::Int)
    parent(f)[i] = val
    return f
end

#####
##### Coordinates of fields
#####

@propagate_inbounds xnode(i, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = xnode(LX(), i, ψ.grid)
@propagate_inbounds ynode(j, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = ynode(LY(), j, ψ.grid)
@propagate_inbounds znode(k, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = znode(LZ(), k, ψ.grid)

@propagate_inbounds Base.lastindex(f::AbstractDataField) = lastindex(f.data)
@propagate_inbounds Base.lastindex(f::AbstractDataField, dim) = lastindex(f.data, dim)

xnodes(ψ::AbstractField) = xnodes(location(ψ, 1), ψ.grid)
ynodes(ψ::AbstractField) = ynodes(location(ψ, 2), ψ.grid)
znodes(ψ::AbstractField) = znodes(location(ψ, 3), ψ.grid)

nodes(ψ::AbstractField; kwargs...) = nodes(location(ψ), ψ.grid; kwargs...)

#####
##### fill_halo_regions!
#####

fill_halo_regions!(field::AbstractField, arch, args...) =
    fill_halo_regions!(field.data, field.boundary_conditions, arch, field.grid, args...)

#####
##### Some conveniences
#####

for f in (:+, :-)
    @eval Base.$f(ϕ::AbstractArray, ψ::AbstractField) = $f(ϕ, interior(ψ))
    @eval Base.$f(ϕ::AbstractField, ψ::AbstractArray) = $f(interior(ϕ), ψ)
end

Base.isapprox(ϕ::AbstractDataField, ψ::AbstractDataField; kw...) = isapprox(interior(ϕ), interior(ψ); kw...)


text:
const GAEVF = JuMP.GenericAffExpr{Float64, VariableRef}
const GQEVF = JuMP.GenericQuadExpr{Float64, VariableRef}

function moi_tests(
    op_problem::OperationsProblem,
    params::Bool,
    vars::Int,
    interval::Int,
    lessthan::Int,
    greaterthan::Int,
    equalto::Int,
    binary::Bool,
)

    JuMPmodel = op_problem.psi_container.JuMPmodel
    @test (:params in keys(JuMPmodel.ext)) == params
    @test JuMP.num_variables(JuMPmodel) == vars
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.Interval{Float64}) == interval
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.LessThan{Float64}) == lessthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.GreaterThan{Float64}) == greaterthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.EqualTo{Float64}) == equalto
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(JuMPmodel)) ==
          binary

    return

end

function psi_constraint_test(
    op_problem::OperationsProblem,
    constraint_names::Vector{Symbol},
)
    for con in constraint_names
        @test !isnothing(get(op_problem.psi_container.constraints, con, nothing))
    end
    return
end

function psi_checkbinvar_test(
    op_problem::OperationsProblem,
    bin_variable_names::Vector{Symbol},
)
    for variable in bin_variable_names
        for v in PSI.get_variable(op_problem.psi_container, variable)
            @test JuMP.is_binary(v)
        end
    end
    return
end

function psi_checkobjfun_test(op_problem::OperationsProblem, exp_type)
    @test JuMP.objective_function_type(op_problem.psi_container.JuMPmodel) == exp_type
    return
end

function moi_lbvalue_test(op_problem::OperationsProblem, con_name::Symbol, value::Number)
    for con in op_problem.psi_container.constraints[con_name]
        @test JuMP.constraint_object(con).set.lower == value
    end
    return
end

function psi_checksolve_test(op_problem::OperationsProblem, status)
    JuMP.optimize!(op_problem.psi_container.JuMPmodel)
    @test termination_status(op_problem.psi_container.JuMPmodel) in status
end

function psi_checksolve_test(
    op_problem::OperationsProblem,
    status,
    expected_result,
    tol = 0.0,
)
    res = solve!(op_problem)
    @test termination_status(op_problem.psi_container.JuMPmodel) in status
    @test isapprox(get_total_cost(res)[:OBJECTIVE_FUNCTION], expected_result, atol = tol)
end


text:
mutable struct Violin <: Geometry
    label::String
    shape::Poly
    nsample::Int
end

function Violin( ; label, shape, nsample)
    return Violin(label, shape, nsample)
end


text:
# p21.jl - eigenvalues of Mathieu operator -u_xx + 2qcos(2x)u
#         (compare p8.jl and p. 724 of Abramowitz & Stegun)

N = 42; h = 2*pi/N; x = h*(1:N);
D2 = toeplitz([-pi^2/(3*h^2)-1/6; @. -.5*(-1)^(1:N-1)/sin(h*(1:N-1)/2)^2]);
qq = 0:.2:15; data = [];
for q = qq;
    e = sort(eigvals(-D2 + 2*q*diagm(cos.(2*x))));
    data = [data; e[1:11]'];
end
clf(); subplot(1,2,1);
plot(qq,data[:,1:2:end],"b-");
plot(qq,data[:,2:2:end],"b--");
xlabel("q"); ylabel("λ");
axis([0,15,-24,32]); yticks(-24:4:32);


text:
# Autogenerated wrapper script for cargo_license_jll for x86_64-linux-gnu
export cargo_license

JLLWrappers.@generate_wrapper_header("cargo_license")
JLLWrappers.@declare_executable_product(cargo_license)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cargo_license,
        "bin/cargo-license",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    MPS

A finite size matrix product state type.
Keeps track of the orthogonality center.
"""
mutable struct MPS <: AbstractMPS
  data::Vector{ITensor}
  llim::Int
  rlim::Int
  function MPS(A::Vector{<:ITensor},
               llim::Int = 0,
               rlim::Int = length(A) + 1)
    new(A, llim, rlim)
  end
end

@doc """
    MPS(v::Vector{<:ITensor})

Construct an MPS from a Vector of ITensors.
""" MPS(v::Vector{<:ITensor})

"""
    MPS()

Construct an empty MPS with zero sites.
"""
MPS() = MPS(ITensor[], 0, 0)

"""
    MPS(N::Int)

Construct an MPS with N sites with default constructed
ITensors.
"""
MPS(N::Int) = MPS(Vector{ITensor}(undef, N))

"""
    MPS([::Type{ElT} = Float64, ]sites)

Construct an MPS filled with Empty ITensors of type `ElT` from a collection of indices.
"""
function MPS(::Type{T}, sites::Vector{<:Index}) where {T <: Number}
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  if N==1
    v[1] = emptyITensor(T,sites[1])
    return MPS(v)
  end

  l = [Index(1, "Link,l=$ii") for ii=1:N-1]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = emptyITensor(T,l[ii], s)
    elseif ii == N
      v[ii] = emptyITensor(T,l[ii-1], s)
    else
      v[ii] = emptyITensor(T,l[ii-1],s,l[ii])
    end
  end
  return MPS(v)
end

MPS(sites::Vector{<:Index}) = MPS(Float64, sites)

function randomU(s1::Index, s2::Index)
  if !hasqns(s1) && !hasqns(s2)
    mdim = dim(s1)*dim(s2)
    RM = randn(mdim,mdim)
    Q,_ = NDTensors.qr_positive(RM)
    G = itensor(Q,dag(s1),dag(s2),s1',s2')
  else
    M = randomITensor(QN(),s1',s2',dag(s1),dag(s2))
    U, S, V = svd(M,(s1',s2'))
    u = commonind(U,S)
    v = commonind(S,V)
    replaceind!(U, u, v)
    G = U * V
  end
  return G
end

function randomizeMPS!(M::MPS, sites::Vector{<:Index}, linkdim=1)
  N = length(sites)
  c = div(N,2)
  max_pass = 100
  for pass=1:max_pass,half=1:2
    if half==1
      (db,brange) = (+1, 1:1:N-1)
    else
      (db,brange) = (-1, N:-1:2)
    end
    for b=brange
      s1 = sites[b]
      s2 = sites[b+db]
      G = randomU(s1,s2)
      T = noprime(G*M[b]*M[b+db])
      rinds = uniqueinds(M[b],M[b+db])
      U,S,V = svd(T,rinds;maxdim=linkdim)
      M[b] = U
      M[b+db] = S*V
      M[b+db] /= norm(M[b+db])
    end
    if half==2 && dim(commonind(M[c],M[c+1])) >= linkdim
      break
    end
  end
  setleftlim!(M, 0)
  setrightlim!(M, 2)
  if dim(commonind(M[c],M[c+1])) < linkdim
    error("MPS center bond dim less than requested")
  end
end

function randomCircuitMPS(::Type{Float64},
                          sites::Vector{<:Index},
                          linkdim::Int;
                          kwargs...)
  N = length(sites)
  M = MPS(N)

  if N==1
    M[1] = ITensor(randn(dim(sites[1])),sites[1])
    M[1] /= norm(M[1])
    return M
  end

  l = Vector{Index}(undef,N)

  d = dim(sites[N])
  chi = min(linkdim,d)
  l[N-1] = Index(chi,"Link,l=$(N-1)")
  O = NDTensors.random_orthog(chi,d)
  M[N] = itensor(O,l[N-1],sites[N])

  for j=N-1:-1:2
    chi *= dim(sites[j])
    chi = min(linkdim,chi)
    l[j-1] = Index(chi,"Link,l=$(j-1)")
    O = NDTensors.random_orthog(chi,dim(sites[j])*dim(l[j]))
    T = reshape(O,(chi,dim(sites[j]),dim(l[j])))
    M[j] = itensor(T,l[j-1],sites[j],l[j])
  end

  O = NDTensors.random_orthog(1,dim(sites[1])*dim(l[1]))
  l0 = Index(1,"Link,l=0")
  T = reshape(O,(1,dim(sites[1]),dim(l[1])))
  M[1] = itensor(T,l0,sites[1],l[1])
  M[1] *= setelt(l0=>1)

  M.llim = 0
  M.rlim = 2

  return M
end

function randomCircuitMPS(sites::Vector{<:Index},
                          linkdim::Int;
                          kwargs...)
  return randomCircuitMPS(Float64, sites, linkdim; kwargs...)
end

"""
    randomMPS(::Type{ElT<:Number}, sites::Vector{<:Index}; linkdim=1)

Construct a random MPS with link dimension `linkdim` of
type `ElT`.
"""
function randomMPS(::Type{ElT},
                   sites::Vector{<:Index},
                   linkdim::Int=1) where {ElT<:Number}
  if hasqns(sites[1])
    error("initial state required to use randomMPS with QNs")
  end

  # For non-QN-conserving MPS, instantiate
  # the random MPS directly as a circuit:
  return randomCircuitMPS(ElT, sites, linkdim)
end

"""
    randomMPS(sites::Vector{<:Index}; linkdim=1)

Construct a random MPS with link dimension `linkdim` of
type `Float64`.
"""
randomMPS(sites::Vector{<:Index},
          linkdim::Int=1) = randomMPS(Float64, sites, linkdim)

"""
    randomMPS(sites::Vector{<:Index}, state; linkdim=1)

Construct a real, random MPS with link dimension `linkdim`,
made by randomizing an initial product state specified by
`state`. This version of `randomMPS` is necessary when creating
QN-conserving random MPS (consisting of QNITensors). The initial
`state` array provided determines the total QN of the resulting
random MPS.
"""
function randomMPS(sites::Vector{<:Index},
                   state,
                   linkdim::Int=1)::MPS
  M = productMPS(sites, state)
  if linkdim > 1
    randomizeMPS!(M,sites,linkdim)
  end
  return M
end

"""
    productMPS(::Type{T<:Number}, ivals::Vector{<:IndexVal})

Construct a product state MPS with element type `T` and
nonzero values determined from the input IndexVals.
"""
function productMPS(::Type{T},
                    ivals::Vector{<:IndexVal}) where {T<:Number}
  N = length(ivals)
  M = MPS(N)

  if N==1
    M[1] = emptyITensor(T,ind(ivals[1]))
    M[1][ivals[1]] = one(T)
    return M
  end

  if hasqns(ind(ivals[1]))
    links = [Index(QN()=>1;tags="Link,l=$n") for n=1:N]
  else
    links = [Index(1,"Link,l=$n") for n=1:N]
  end
  M[1] = emptyITensor(T,ind(ivals[1]), links[1])
  M[1][ivals[1],links[1](1)] = one(T)
  for n=2:N-1
    s = ind(ivals[n])
    M[n] = emptyITensor(T,dag(links[n-1]),s,links[n])
    M[n][links[n-1](1),ivals[n],links[n](1)] = one(T)
  end
  M[N] = emptyITensor(T,dag(links[N-1]),ind(ivals[N]))
  M[N][links[N-1](1),ivals[N]] = one(T)
  return M
end

"""
    productMPS(ivals::Vector{<:IndexVal})

Construct a product state MPS with element type `Float64` and
nonzero values determined from the input IndexVals.
"""
productMPS(ivals::Vector{<:IndexVal}) =
  productMPS(Float64,
             ivals::Vector{<:IndexVal})

"""
    productMPS(::Type{T},
               sites::Vector{<:Index},
               states::Union{Vector{String},
                             Vector{Int},
                             String,
                             Int})

Construct a product state MPS of element type `T`, having
site indices `sites`, and which corresponds to the initial
state given by the array `states`. The input `states` may
be an array of strings or an array of ints recognized by the
`state` function defined for the relevant Index tag type.
In addition, a single string or int can be input to create
a uniform state.

# Examples
```julia
N = 10
sites = siteinds("S=1/2", N)
states = [isodd(n) ? "Up" : "Dn" for n=1:N]
psi = productMPS(ComplexF64, sites, states)
phi = productMPS(sites, "Up")
```
"""
function productMPS(::Type{T},
                    sites::Vector{ <: Index},
                    states) where {T <: Number}
  if length(sites) != length(states)
    throw(DimensionMismatch("Number of sites and and initial states don't match"))
  end
  ivals = [state(sites[n],states[n]) for n=1:length(sites)]
  return productMPS(T, ivals)
end

function productMPS(::Type{T},
                    sites::Vector{ <: Index},
                    states::Union{String, Int}) where {T <: Number}
  ivals = [state(sites[n], states) for n in 1:length(sites)]
  return productMPS(T, ivals)
end

"""
    productMPS(sites::Vector{<:Index},states)

Construct a product state MPS having
site indices `sites`, and which corresponds to the initial
state given by the array `states`. The `states` array may
consist of either an array of integers or strings, as
recognized by the `state` function defined for the relevant
Index tag type.

# Examples
```julia
N = 10
sites = siteinds("S=1/2",N)
states = [isodd(n) ? "Up" : "Dn" for n=1:N]
psi = productMPS(sites,states)
```
"""
productMPS(sites::Vector{ <: Index},
           states) =
  productMPS(Float64, sites, states)

"""
    siteind(M::MPS, j::Int; kwargs...)

Get the site Index of the MPS.
"""
siteind(M::MPS, j::Int; kwargs...) = firstsiteind(M, j; kwargs...)

"""
    siteinds(M::MPS)

Get a vector of the site indices of the MPS.
"""
siteinds(M::MPS) = [siteind(M, j) for j in 1:length(M)]

function replace_siteinds!(M::MPS, sites)
  for j in eachindex(M)
    sj = siteind(M, j)
    replaceind!(M[j], sj, sites[j])
  end
  return M
end

replace_siteinds(M::MPS, sites) = replace_siteinds!(copy(M), sites)

"""
    replacebond!(M::MPS, b::Int, phi::ITensor; kwargs...)

Factorize the ITensor `phi` and replace the ITensors
`b` and `b+1` of MPS `M` with the factors. Choose
the orthogonality with `ortho="left"/"right"`.
"""
function replacebond!(M::MPS, b::Int, phi::ITensor;
                      kwargs...)
  ortho::String = get(kwargs, :ortho, "left")
  swapsites::Bool = get(kwargs, :swapsites, false)
  which_decomp::Union{String, Nothing} = get(kwargs, :which_decomp, nothing)
  normalize::Bool = get(kwargs, :normalize, false)

  # Deprecated keywords
  if haskey(kwargs, :dir)
    error("""dir keyword in replacebond! has been replaced by ortho.
          Note that the options are now the same as factorize, so use `left` instead of `fromleft` and `right` instead of `fromright`.""")
  end

  indsMb = inds(M[b])
  if swapsites
    sb = siteind(M, b)
    sbp1 = siteind(M, b+1)
    indsMb = replaceind(indsMb, sb, sbp1)
  end

  L,R,spec = factorize(phi, indsMb;
                       which_decomp = which_decomp,
                       tags = tags(linkind(M, b)),
                       kwargs...)

  M[b]   = L
  M[b+1] = R
  if ortho == "left"
    leftlim(M) == b-1 && setleftlim!(M, leftlim(M)+1)
    rightlim(M) == b+1 && setrightlim!(M, rightlim(M)+1)
    normalize && (M[b+1] ./= norm(M[b+1]))
  elseif ortho == "right"
    leftlim(M) == b && setleftlim!(M, leftlim(M)-1)
    rightlim(M) == b+2 && setrightlim!(M, rightlim(M)-1)
    normalize && (M[b] ./= norm(M[b]))
  else
    error("In replacebond!, got ortho = $ortho, only currently supports `left` and `right`.")
  end
  return spec
end

"""
    replacebond(M::MPS, b::Int, phi::ITensor; kwargs...)

Like `replacebond!`, but returns the new MPS.
"""
function replacebond(M0::MPS, b::Int, phi::ITensor;
                     kwargs...)
  M = copy(M0)
  replacebond!(M, b, phi; kwargs...)
  return M
end

"""
    sample!(m::MPS)

Given a normalized MPS m, returns a `Vector{Int}`
of `length(m)` corresponding to one sample
of the probability distribution defined by
squaring the components of the tensor
that the MPS represents. If the MPS does
not have an orthogonality center,
orthogonalize!(m,1) will be called before
computing the sample.
"""
function sample!(m::MPS)
  orthogonalize!(m, 1)
  return sample(m)
end

"""
    sample(m::MPS)

Given a normalized MPS m with `orthocenter(m)==1`,
returns a `Vector{Int}` of `length(m)`
corresponding to one sample of the
probability distribution defined by
squaring the components of the tensor
that the MPS represents
"""
function sample(m::MPS)
  N = length(m)

  if orthocenter(m) != 1
    error("sample: MPS m must have orthocenter(m)==1")
  end
  if abs(1.0-norm(m[1])) > 1E-8
    error("sample: MPS is not normalized, norm=$(norm(m[1]))")
  end

  result = zeros(Int,N)
  A = m[1]

  for j=1:N
    s = siteind(m,j)
    d = dim(s)
    # Compute the probability of each state
    # one-by-one and stop when the random
    # number r is below the total prob so far
    pdisc = 0.0
    r = rand()
    # Will need n,An, and pn below
    n = 1
    An = ITensor()
    pn = 0.0
    while n <= d
      projn = ITensor(s)
      projn[s[n]] = 1.0
      An = A*projn
      pn = real(scalar(dag(An)*An))
      pdisc += pn
      (r < pdisc) && break
      n += 1
    end

    result[j] = n

    if j < N
      A = m[j+1]*An
      A *= (1.0/sqrt(pn))
    end
  end
  return result
end

function HDF5.write(parent::Union{HDF5File,HDF5Group},
                    name::AbstractString,
                    M::MPS)
  g = g_create(parent,name)
  attrs(g)["type"] = "MPS"
  attrs(g)["version"] = 1
  N = length(M)
  write(g, "length", N)
  write(g, "rlim", M.rlim)
  write(g, "llim", M.llim)
  for n=1:N
    write(g,"MPS[$(n)]", M[n])
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{MPS})
  g = g_open(parent,name)
  if read(attrs(g)["type"]) != "MPS"
    error("HDF5 group or file does not contain MPS data")
  end
  N = read(g, "length")
  rlim = read(g, "rlim")
  llim = read(g, "llim")
  v = [read(g,"MPS[$(i)]",ITensor) for i in 1:N]
  return MPS(v, llim, rlim)
end


text:
using Devectorize

function update!(net::Network)
    for i in 1:length(net.layers)
        update!(net.layers[i], net.optimizers[i], net.rnn_optimizers[i])
    end
end

function update!(layer::Layer, opt::Optimizer, dummy::Optimizer)
    update_weight!(opt, layer.W, layer.ΔW)
    update_bias!(opt, layer.b, layer.Δb)
end

function update!(layer::RecurrentLayer, opt::Optimizer, rnn_opt::Optimizer)
    update_weight!(opt, layer.W, layer.ΔW)
    update_weight!(rnn_opt, layer.Wh, layer.ΔWh)
    update_bias!(opt, layer.b, layer.Δb)
end

function update_bias!(sgd::Optimizer, b, Δb)
    db = @in1! Δb .* -0.3
    nb = @in1! b   .+ Δb
end

# --------------------------------------------------------------------------------------------
immutable SGD <: Optimizer
    λ::Float32            # weight decrese
    ε::Float32            # learning rate

    function SGD(;λ=0.01, ε=1.0)
        new(λ, ε)
    end
end

function update_weight!(sgd::SGD, W, ΔW)
    ε=sgd.ε; λ=sgd.λ

    @inbounds for col=1:size(W,2), row=1:size(W,1)
        ΔW[row, col] += W[row, col] * λ
        ΔW[row, col] *= -ε
        W[row, col]   += ΔW[row, col]
    end
end

# --------------------------------------------------------------------------------------------
immutable AdaGrad <: Optimizer
    λ::Float32  # weight decrese
    ε::Float32  # learning rate
    μ::Float32  # RMSprop
    G::Weight    # sum of gradients
    H::Weight    # sum of gradients^2

    function AdaGrad(in_size, out_size; ε=0.0, λ=0.01, μ=0.0)
        new(λ, ε, μ, zeros(out_size, in_size), zeros(out_size, in_size))
    end
end

function update_weight!(ada::AdaGrad, W, ΔW)
    ε=ada.ε; λ=ada.λ; μ=ada.μ

    new_g = @in1! ada.G .+ ΔW

    if μ > 0.0
        copy!(ada.H, (μ * (ΔW .^ 2) + (1.0 - μ) * ada.H)) # RMSprop
    else
        new_h = @in1! ada.H .+ (ΔW .^ 2)                    # AdaGradRDA
    end

    H = sqrt(ada.H)
    for col=1:size(W,2), row=1:size(W,1)
        W[row,col] += (-ε * ΔW[row,col]) / (1.0 + H[row,col])

        # FOBOS
        if λ > 0.0
            if abs(W[row,col]) < λ
                W[row,col] = 0.0
            elseif W[row,col] > 0.0
                W[row,col] -= λ
            else
                W[row,col] += λ
            end
        end
    end
end

# --------------------------------------------------------------------------------------------
type AdaGradRDA <: Optimizer
    λ::Float32  # weight decrese
    μ::Float32  # RMSprop
    G::Weight    # sum of gradients
    H::Weight    # sum of gradients^2
    num_of_grads::Int64 # number of gradients

    function AdaGradRDA(in_size, out_size; λ=0.01, μ=0.0)
        new(λ, μ, zeros(out_size, in_size), zeros(out_size, in_size), 0)
    end
end

function update_weight!(ada::AdaGradRDA, W, ΔW)
    λ=ada.λ; μ=ada.μ

    ada.num_of_grads += 1
    ada.G            += ΔW

    if μ > 0.0
        ada.H = μ * (ΔW .^ 2) + (1.0 - μ) * ada.H  # RMSprop
    else
        ada.H += ΔW .^ 2                             # AdaGradRDA
    end

    G = ada.G ./ ada.num_of_grads

    for col=1:size(W,2), row=1:size(W,1)
        if abs(G[row,col]) < λ
            W[row,col] = 0.0
        elseif G[row,col] > λ
            W[row,col] = (-1.0 * G[row,col] * ada.num_of_grads + λ) / (1.0 + ada.H[row,col])
        else
            W[row,col] = (-1.0 * G[row,col] * ada.num_of_grads - λ) / (1.0 + ada.H[row,col])
        end
    end
end


text:
# The precision of forward-differiation is accurate at Float64 but not at Float32/16.
# Anyway the backward-differentiation is accurate at Float64/32/16.
# Just check at Float64/32/16 while keeping the same random seed.
# For the below example the gradients for CPU/GPU on my device are
# Float64 is 100.724249419
# Float32 is 100.72425
# Float16 is 100.75
# ochannels = 1

@testset "check gradient for PlainConv1d block" begin
    clear()
    # using Random
    # Random.seed!(1234);
    TYPE = Array{Float64};

    # [0] prepare model
    ichannels = 1;
    ochannels = 2;
    c1 = PlainConv1d(ichannels,4,3; type=TYPE);
    c2 = PlainConv1d(4,ochannels,2; type=TYPE);

    # [1] prepare input data and its label
    timeSteps = 128;
    batchsize = 32;
    x = Variable(rand(ichannels,timeSteps,batchsize); type=TYPE);
    l = Variable(rand(ochannels,      125,batchsize); type=TYPE);

    # [2] forward and backward propagation
    o1 = forward(c1,  x);
    o2 = forward(c2, o1);
    COST1 = mseCost(o2, l);
    backward();

    # [3] with a samll change of a weight
    GRAD = c1.w.delta[1];
    DELTA = 1e-6;
    c1.w.value[1] += DELTA;

    # [4] forward and backward propagation
    o1 = forward(c1,  x);
    o2 = forward(c2, o1);
    COST2 = mseCost(o2, l);
    backward();

    # [5] check if the auto-grad is true or not
    dLdW = (COST2 - COST1)/DELTA;   # numerical gradient
    err  = abs((dLdW-GRAD)/(GRAD+eps(Float64)))*100;  # relative error in %
    err  = err < 1e-1 ? 0.0 : err;
    @test err < 0.1
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Zlib_jll

@testset "Zlib_jll" begin
    @test VersionNumber(unsafe_string(ccall((:zlibVersion, libz), Cstring, ()))) == v"1.2.11"
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
struct LJRefParam <: EoSParam
    epsilon::PairParam{Float64}
    sigma::PairParam{Float64}
    Mw::SingleParam{Float64}
end

struct LJRefConsts <: EoSParam
    n::Vector{Float64}
    t::Vector{Float64}
    d::Vector{Int}
    c::Vector{Int}
    beta::Vector{Float64}
    gamma::Vector{Float64}
    eta::Vector{Float64}
    epsilon::Vector{Float64}
    function LJRefConsts()
        n = [
        0.52080730e-2,
        0.21862520e+1,
        -0.21610160e+1,
        0.14527000e+1,
        -0.20417920e+1,
        0.18695286e0,
        -0.90988445e-1,
        -0.49745610e0,
        0.10901431e0,
        -0.80055922e0,
        -0.56883900e0,
        -0.62086250e0,
        -0.14667177e+1,
        0.18914690e+1,
        -0.13837010e0,
        -0.38696450e0,
        0.12657020e0,
        0.60578100e0,
        0.11791890e+1,
        -0.47732679e0,
        -0.99218575e+1,
        -0.57479320e0,
        0.37729230e-2,
        ]

        t = [
        1.000
        0.320
        0.505
        0.672
        0.843
        0.898
        1.294
        2.590
        1.786
        2.770
        1.786
        1.205
        2.830
        2.548
        4.650
        1.385
        1.460
        1.351
        0.660
        1.496
        1.830
        1.616
        4.970
        ]

        d = [
        4
        1
        1
        2
        2
        3
        5
        2
        2
        3
        1
        1
        1
        1
        2
        3
        3
        2
        1
        2
        3
        1
        1
        ]
        c =
        [
        1
        2
        1
        2
        2
        1

        ]

        g = [
        2.067 0.625 0.710 0.2053
        1.522 0.638 0.860 0.4090
        8.820 3.910 1.940 0.6000
        1.722 0.156 1.480 1.2030
        0.679 0.157 1.490 1.8290
        1.883 0.153 1.945 1.3970
        3.925 1.160 3.020 1.3900
        2.461 1.730 1.110 0.5390
        28.20 383.0 1.170 0.9340
        0.753 0.112 1.330 2.3690
        0.820 0.119 0.240 2.4300]

        eta = g[:,1]
        beta = g[:,2]
        gamma = g[:,3]
        epsilon = g[:,4]
        return new(n,t,d,c,beta,gamma,eta,epsilon)
    end

end

is_splittable(::LJRefConsts) = false

struct LJRef <: EmpiricHelmholtzModel
    components::Vector{String}
    icomponents::UnitRange{Int}
    params::LJRefParam
    consts::LJRefConsts
    references::Vector{String}
end
@registermodel LJRef

idealmodel(x::LJRef) = x
export LJRef

function LJRef(components;
    userlocations=String[],
    verbose=false)
    params,sites = getparams(components, ["SAFT/PCSAFT"]; userlocations=userlocations, verbose=verbose)
    Mw = params["Mw"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], params["k"])
    params = LJRefParam(epsilon,sigma,Mw)
    consts = LJRefConsts()
    icomponents = 1:length(components)
    references = ["10.1063/1.4945000"]
    return LJRef(components,icomponents,params,consts,references)
end

function _f0(model::LJRef,ρ,T,z=SA[1.0],∑z = sum(z))
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    lnΣz = log(∑z)
    res = zero(ρ+T+first(z))
    for i  ∈ @comps
        τᵢ = 1.32/(T/ϵ[i])
        δᵢ = (N_A*ρ*σ[i]^3)/0.31
        aᵢ = log(δᵢ) + 1.5*log(τᵢ) + 1.515151515*τᵢ + 6.262265814
        res += z[i]*(aᵢ + log(z[i]) - lnΣz)
    end
    return res
end

function _fr(model::LJRef,δ,τ)
    ai = zero(δ+τ)
    n = model.consts.n
    t = model.consts.t
    d = model.consts.d
    c = model.consts.c
    β = model.consts.beta
    γ = model.consts.gamma
    η = model.consts.eta
    ε = model.consts.epsilon

    @inbounds begin
        for k ∈ 1:6
            ai += n[k]*(δ^d[k])*(τ^t[k])
        end
        for (k,k_) ∈ zip(7:12,1:6)
            ai += n[k]*(δ^d[k])*(τ^t[k])*exp(-δ^c[k_])
        end

        for (k,k_) ∈ zip(13:23,1:11)
            ai += n[k]*(δ^(d[k]))*(τ^(t[k]))*
            exp(-η[k_]*(δ - ε[k_])^2 - β[k_]*(τ -γ[k_])^2)
        end
    end
    return ai
end

#TODO: better relations? EoSRef was done with one fluid in mind.
#this is technically an unsafe extension.
function _v_scale(model::LJRef,z=SA[1.0])
    σᵢᵢ = model.params.sigma.diagvalues
    val = N_A*sum(z[i]*σᵢᵢ[i]^3 for i in 1:length(z))
    return val
end

function T_scale(model::LJRef,z=SA[1.0])
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    σϵ_mix = zero(eltype(z))
    σ_mix = zero(eltype(z))
    comps = length(model)
    for i ∈ 1:comps
        zi = z[i]
        zii = zi*zi
        σ3 = σ[i,i]^3
        σϵ_mix += zii*σ3*ϵ[i,i]
        σ_mix += zii*σ3
        for j ∈ 1:(i-1)
            σ3ij = σ[i,j]^3
            zij = zi*z[j]
            σϵ_mix += 2*zij*σ3ij*ϵ[i,j]
            σ_mix += 2*zij*σ3ij*ϵ[i,j]
        end
    end
    return σϵ_mix/σ_mix
end

function lb_volume(model::LJRef, z = SA[1.0])
    val = π/6*_v_scale(model,z)
    return val
end

function eos(model::LJRef,V,T,z = SA[1.0])
    Σz = sum(z)
    ρ = Σz/V
    α0 = _f0(model,ρ,T,Σz)
    τ = 1.32/(T/T_scale(model,z))
    δ = (ρ*_v_scale(model,z)/Σz)/0.31
    αr =  _fr(model,δ,τ)
    x1 = R̄*T*Σz*αr
    x2 =  R̄*T*α0
    return x1+x2
end

function a_ideal(model::LJRef,V,T,z = SA[1.0])
    Σz = sum(z)
    ρ = Σz/V
    α0 = _f0(model,ρ,T,Σz)
    return α0/Σz
end

function a_res(model::LJRef,V,T,z = SA[1.0])
    Σz = sum(z)
    ρ = Σz/V
    α0 = _f0(model,ρ,T,Σz)
    τ = 1.32/(T/T_scale(model,z))
    δ = (ρ*_v_scale(model,z)/Σz)/0.31
    return  _fr(model,δ,τ)
end

function eos_res(model::LJRef,V,T,z = SA[1.0])
    Σz = sum(z)
    ρ = Σz/V
    α0 = _f0(model,ρ,T,Σz)
    τ = 1.32/(T/T_scale(model,z))
    δ = (ρ*_v_scale(model,z)/Σz)/0.31
    αr =  _fr(model,δ,τ)
    return R̄*T*Σz*αr
end
#=
function ljref_psat(Tr,pc)
    n = (0.54000e+1,0.44704e01,-0.18530e+1,0.19890e0,-0.11250e+1)
    t = (1.,1.5,4.7,2.5,21.4)
    tr1 = one(Tr) - Tr
    res = sum(ni*tr1^ti for (ni,ti) in zip(n,t))
    return exp(res/Tr)*pc
end
=#
function ljref_rholsat(Tr)
    n = (0.1362e+1,0.2093e+1,-0.2110e+1,0.3290e0,0.1410e+1)
    t = (0.313 ,0.940,1.630,17.,2.4)
    tr1 = one(Tr) - Tr
    res = sum(ni*tr1^ti for (ni,ti) in zip(n,t))
    rhoc = 0.31
    return (1+res)*rhoc
end

function ljref_rhovsat(Tr)
    n = (-0.69655e+1,-0.10331e+3,-0.20325e+1,-0.44481e+2,-0.18463e+2,-0.26070e+3)
    t = (1.320 ,19.24,0.360,8.780,4.040,41.60)
    tr1 = one(Tr) - Tr
    res = sum(ni*tr1^ti for (ni,ti) in zip(n,t))
    rhoc = 0.31
    return exp(res)*rhoc
end

function x0_sat_pure(model::LJRef,T)
    σ = model.params.sigma.values[1]
    Tc = 1.32*T_scale(model)
    ρl =  ljref_rholsat(T/Tc)/(N_A*σ^3)
    ρv =  ljref_rhovsat(T/Tc)/(N_A*σ^3)
    return [log10(1/ρl),log10(1/ρv)]
end

function p_scale(model::LJRef,z = SA[1.0])
    rhoc = 1/(_v_scale(model,z))
    Tc = T_scale(model,z)
    return R̄*Tc*rhoc
end


text:
module WaveEquationCubedSphereTests

using Test
using Gridap
using GridapGeosciences
using Plots
using LinearAlgebra
using WriteVTK
using DataFrames
using CSV

# Initial depth
function h₀(xyz)
  x,y,z = xyz
  xᵢ,yᵢ,zᵢ = (1,0,0)
  1.0 + 0.0001*exp(-5*((xᵢ-x)^2+(yᵢ-y)^2+(zᵢ-z)^2))
end

# Initial velocity
u₀(xyz) = zero(xyz)

function new_vtk_step(Ω,file,hn,un)
  createvtk(Ω,
            file,
            cellfields=["hn"=>hn, "un"=>un],
            nsubcells=4)
end

function _ssrk2_update!(res,MM,MMchol,DIV,α,a,b)
  mul!(res, MM, a)            # res <- MM*a
  mul!(res, DIV, b, α, 1.0)   # res <- 1.0*res + α*DIV*b
  ldiv!(MMchol,res)           # res <- inv(MM)*res
end

function generate_energy_plots(outdir,N,ke,pe,kin_to_pot,pot_to_kin)
  plot(1:N,
       [ke,pe],
       label=["Ek (kinetic)" "Ep (potential)"],
       xlabel="Step",ylabel="Energy")
  savefig(joinpath(outdir,"energy.png"))
  plot(1:N,
       [kin_to_pot,-pot_to_kin,kin_to_pot-pot_to_kin],
       label=["hⁿ∇⋅uⁿ" "uⁿ⋅∇(hⁿ)" "Balance"],
       xlabel="Step",ylabel="Kinetic to potential energy balance",legend = :outertopleft)
  savefig(joinpath(outdir,"kinetic_to_potential_energy_balance.png"))
end

"""
  Solves the wave equation using a 2nd order
  Strong-Stability-Preserving Runge-Kutta
  explicit time integration scheme (SSRK2)

  g : acceleration due to gravity
  H : (constant) reference layer depth
  T : [0,T] simulation interval
  N : number of time subintervals
"""
function solve_wave_equation_ssrk2(
      model,order,degree,g,H,T,N;
      write_results=false,
      out_dir="wave_eq_ncells_$(num_cells(model))_order_$(order)_ssrk2",
      out_period=N/10)

  RT=ReferenceFE(raviart_thomas,Float64,order)
  DG=ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,RT; conformity=:Hdiv)
  Q = FESpace(model,DG; conformity=:L2)

  U = TrialFESpace(V)
  P = TrialFESpace(Q)

  Ω  = Triangulation(model)
  dΩ = Measure(Ω,degree)
  dω = Measure(Ω,degree,ReferenceDomain())

  # Build mass matrix in RT and DG spaces
  # and their sparse Cholesky factors
  amm(u,v) = ∫(v⋅u)dΩ
  RTMM=assemble_matrix(amm,U,V)
  L2MM=assemble_matrix(amm,P,Q)
  RTMMchol=lu(RTMM)
  L2MMchol=lu(L2MM)

  # Build g*DIV(v)*h and H*q*DIV(u)
  ad(h,v) = ∫(DIV(v)*h)dω
  divvh=assemble_matrix(ad,P,V)
  adt(u,q) = ∫(q*DIV(u))dω
  qdivu=assemble_matrix(adt,U,Q)

  # Interpolate initial condition into FE spaces
  hn=interpolate_everywhere(h₀,P); hnv=get_free_dof_values(hn)
  un=interpolate_everywhere(u₀,U); unv=get_free_dof_values(un)
  function run_simulation(pvd=nothing)
    # Allocate work space vectors
    h1v = similar(get_free_dof_values(hn))
    h2v = similar(h1v)
    u1v = similar(get_free_dof_values(un))
    u2v = similar(u1v)
    if (write_results)
      ke  = Vector{Float64}(undef,N)
      pe  = Vector{Float64}(undef,N)
      kin_to_pot = Vector{Float64}(undef,N)
      pot_to_kin = Vector{Float64}(undef,N)
      mass = Vector{Float64}(undef,N)
      initialize_csv(joinpath(out_dir,"wave_eq_geosciences_data.csv"),
                  "time", "hn_dot_div_un", "un_dot_grad_hn",
                  "mass", "kinetic", "potential")
    end
    dt  = T/N
    dtg = dt*g
    dtH = dt*H

    for step=1:N
       # 1st step
       # inv(L2MM)*(L2MM*hnv - dt*Hqdivu*unv)
       # inv(RTMM)*(RTMM*unv + dt*gdivvh*hnv)
       _ssrk2_update!(h1v, L2MM, L2MMchol, qdivu,-dtH, hnv, unv)
       _ssrk2_update!(u1v, RTMM, RTMMchol, divvh, dtg, unv, hnv)

       # 2nd step
       # inv(L2MM)*(L2MM*h1v - dt*Hqdivu*u1v)
       # inv(RTMM)*(RTMM*u1v + dt*gdivvh*h1v)
       _ssrk2_update!(h2v, L2MM, L2MMchol, qdivu, -dtH, h1v, u1v)
       _ssrk2_update!(u2v, RTMM, RTMMchol, divvh,  dtg, u1v, h1v)

       # Averaging steps
       hnv .= 0.5 .* ( hnv .+ h2v )
       unv .= 0.5 .* ( unv .+ u2v )

       if (write_results)
         ke[step]=Eₖ(un,H,dΩ)
         pe[step]=Eₚ(hn,hn,0.5*g,dΩ)
         kin_to_pot[step]=compute_kin_to_pot!(u1v,unv,divvh,hnv)
         pot_to_kin[step]=compute_pot_to_kin!(h1v,hnv,qdivu,unv)
         mass[step] = compute_total_mass!(h1v,L2MM,hnv)
         if mod(step, out_period) == 0
           println(step)
           pvd[dt*Float64(step)] = new_vtk_step(Ω,joinpath(out_dir,"n=$(step)"),hn,un)
         end
         # save global scalar snapshots
         append_to_csv(joinpath(out_dir,"wave_eq_geosciences_data.csv");
                        time = step*dt,
                        hn_dot_div_un = kin_to_pot[step],
                        un_dot_grad_hn = pot_to_kin[step],
                        mass = mass[step],
                        kinetic = ke[step],
                        potential= pe[step])
       end
    end
    if (write_results)
      pvd[dt*Float64(N)] = new_vtk_step(Ω,joinpath(out_dir,"n=$(N)"),hn,un)
      vtk_save(pvd)
      generate_energy_plots(out_dir,N,ke,pe,kin_to_pot,pot_to_kin)
    end
    un,hn
  end
  if (write_results)
    rm(out_dir,force=true,recursive=true)
    mkdir(out_dir)
    pvdfile=joinpath(out_dir,"wave_eq_ncells_$(num_cells(model))_order_$(order)_ssrk2")
    paraview_collection(run_simulation,pvdfile)
  else
    run_simulation()
  end
end

model=CubedSphereDiscreteModel(12)
g=1.0
H=1.0
T=2π
N=2000
order=0
degree=4
@time un,hn =
  solve_wave_equation_ssrk2(model,order,degree,g,H,T,N;write_results=false,out_period=10)

@test Eₖ(un,H,Measure(Triangulation(model),degree)) ≈ 1.6984501177784049e-10

end # module


text:
function PVCoordinates()
    return PVCoordinates(())
end

function PVCoordinates(arg0::FieldVector3D)
    return PVCoordinates((FieldVector3D,), arg0)
end

function PVCoordinates(arg0::PVCoordinates, arg1::PVCoordinates)
    return PVCoordinates((PVCoordinates, PVCoordinates), arg0, arg1)
end

function PVCoordinates(arg0::Vector3D, arg1::Vector3D)
    return PVCoordinates((Vector3D, Vector3D), arg0, arg1)
end

function PVCoordinates(arg0::Vector3D, arg1::Vector3D, arg2::Vector3D)
    return PVCoordinates((Vector3D, Vector3D, Vector3D), arg0, arg1, arg2)
end

function PVCoordinates(arg0::jdouble, arg1::PVCoordinates)
    return PVCoordinates((jdouble, PVCoordinates), arg0, arg1)
end

function PVCoordinates(arg0::jdouble, arg1::PVCoordinates, arg2::jdouble, arg3::PVCoordinates)
    return PVCoordinates((jdouble, PVCoordinates, jdouble, PVCoordinates), arg0, arg1, arg2, arg3)
end

function PVCoordinates(arg0::jdouble, arg1::PVCoordinates, arg2::jdouble, arg3::PVCoordinates, arg4::jdouble, arg5::PVCoordinates)
    return PVCoordinates((jdouble, PVCoordinates, jdouble, PVCoordinates, jdouble, PVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5)
end

function PVCoordinates(arg0::jdouble, arg1::PVCoordinates, arg2::jdouble, arg3::PVCoordinates, arg4::jdouble, arg5::PVCoordinates, arg6::jdouble, arg7::PVCoordinates)
    return PVCoordinates((jdouble, PVCoordinates, jdouble, PVCoordinates, jdouble, PVCoordinates, jdouble, PVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function shifted_by(obj::PVCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", PVCoordinates, (jdouble,), arg0)
end

function to_string(obj::PVCoordinates)
    return jcall(obj, "toString", JString, ())
end


text:
################################################################################

@testset verbose = true "Distance" begin

  ################################################################################

  # testing
  @testset "Euclidean distance" begin
    @test euclDist(collect(1:10), fill(5, 10)) == 9.219544457292887
    @test euclDist(collect(1:10), fill(7, 10)) == 10.246950765959598
    @test euclDist(collect(1:10), fill(1, 10)) == 16.881943016134134
  end

  ################################################################################

  # testing
  @testset "Bhattacharyya distance" begin
    @test bhattDist(collect(1:10), fill(5, 10)) == -3.936532135073928
    @test bhattDist(collect(1:10), fill(7, 10)) == -4.101741743301577
    @test bhattDist(collect(1:10), fill(1, 10)) == -3.1556496448455857
  end

  ################################################################################

  # testing
  @testset "Amplitude" begin
    @test HiddenMarkovModelReaders.amplitude(collect(1:10)) == 19.621416870348583
  end

  ################################################################################

end

################################################################################

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
# Important packages
println("== Compiling packages...")
using JLD

# Loading the array
println("== Loading data...")
count = load("backup_matrices/counts.jld")["counts"]

# Checks if the count array is symmetric
isSymmetric = Symmetric(count) == count
println("== The count array is symmetric: ", isSymmetric)
if (isSymmetric == false)
	exit
end

# Normalizing the count array
println("== Normalizing count array (ETA: 1h30 from start)...")
count /= norm(count)

# Saving the normalized array
println("== Exporting normalized array...")
save("backup_matrices/normalized.jld", "norm", count)

# Decomposing w/ eigendecomposition
println("== Eigendecomposition (ETA: 50min from start)...")
d, v = eig(count)

# Saving matrices
println("== Exporting eigenvalues and eigenvectors arrays...")
save("backup_matrices/eigenvalues.jld", "values", d)
save("backup_matrices/eigenvectors.jld", "vectors", v)


text:
"""
    cm2inch(x)

Convert length in cm to inch units: 1 inch = 2.54 cm.
"""
cm2inch(x) = x/2.54

"""
    embed_system(m;name)

Embeds system `m` as the single component of a larger system.
This helps to match the naming of the states, parameters, observables to
the namespace of the system.

```jldocstest; output=false
using ModelingToolkit, DifferentialEquations
using MTKHelpers
function samplesystem(;name)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)  # RHS is observed
    ps = @parameters τ       # parameters
    ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end
@named m = samplesystem()
@named sys = embed_system(m)

# note that keys are `m.x`,`m.τ` or `m.RHS`.
# Hence only m needs to be defined rather than all the states, parameters,
# and observables of the system.
prob = ODEProblem(sys, [m.x => 0.0], (0.0,10.0), [m.τ => 3.0])
sol = solve(prob);
dx1 = sol[m.RHS][1:3]
dx2 = sol[getproperty(m,:RHS)][1:3]  # access by symbols

# using Plots
# plot(sol, vars=[m.x,m.RHS])
# plot(sol, vars=getproperty.(Ref(m),[:x, :RHS])) # access by symbols
length(dx2) == 3
# output
true
```
"""
function embed_system(m;name)
    @named _sys_embed = ODESystem(Equation[], ModelingToolkit.get_iv(m))
    sys = compose(_sys_embed, [m]; name)
    structural_simplify(sys)
end

"""
    symbol(t)

Extract the inner symbol from a Term or Num object.
"""
function symbol(t::Term); symbol(t.f); end
symbol(num::Num) = symbol(num.val)
symbol(s) = Symbol(s)

"""
    strip_namespace(s)

Omit the part before the first dot.
"""
function strip_namespace(s::String); match(r"[^.₊]+$",s).match; end
function strip_namespace(s::Symbol); Symbol(strip_namespace(string(s))); end

"""
    statesyms(sys::ODYSystem)
    parsyms(sys::ODYSystem)

Extract the basic symbols without namespace of system states and system parameters.
"""
function statesyms(sys::ODESystem); symbol.(states(sys)); end
parsyms(sys::ODESystem) = symbol.(parameters(sys))

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


text:
"""
**PACOTE DIDÁTICO PARA CALCULAR A ÁREA E VOLUME DE CILINDROS.**

SINTAXE:

arealateral(raio, altura_cilindro)

areatotal(raio, altura_cilindro)

arealateral(raio, altura_cilindro)

volcilindro(raio, altura_cilindro)

"""
module CalCL

# O código principal do pacote

include("areaBC.jl")
include("arealateral.jl")
include("volcilindro.jl")
include("areatotal.jl")

end


text:
@testset "Global" begin

function kernel(X)
    ptr = AMDGPUnative.get_global_pointer(Val(:myglobal), Float32)
    Base.unsafe_store!(ptr, 3f0)
    nothing
end

hk = AMDGPUnative.rocfunction(kernel, Tuple{Int32})
exe = hk.mod.exe
gbl = HSARuntime.get_global(exe.exe, :myglobal)
gbl_ptr = Base.unsafe_convert(Ptr{Float32}, gbl.ptr)
@test Base.unsafe_load(gbl_ptr) == 0f0
Base.unsafe_store!(gbl_ptr, 2f0)
@test Base.unsafe_load(gbl_ptr) == 2f0

dev = AMDGPUnative.default_device()
queue = AMDGPUnative.default_queue(dev)
kern = AMDGPUnative.create_kernel(dev, exe, hk.fun.entry, (Int32(3),))
signal = AMDGPUnative.create_event()
AMDGPUnative.launch_kernel(queue, kern, signal; groupsize=(1,1,1), gridsize=(1,1,1))
wait(signal)
@test Base.unsafe_load(gbl_ptr) == 3f0

end


text:
# This file was generated, do not modify it. # hide
x1 = Wtrain.x1
x2 = Wtrain.x2

mask_1 = ytrain .== 1
mask_2 = ytrain .== 2
mask_3 = ytrain .== 3

figure(figsize=(8, 6))
plot(x1[mask_1], x2[mask_1], linestyle="none", marker="o", color="red")
plot(x1[mask_2], x2[mask_2], linestyle="none", marker="o", color="blue")
plot(x1[mask_3], x2[mask_3], linestyle="none", marker="o", color="magenta")

xlabel("PCA dimension 1", fontsize=14)
ylabel("PCA dimension 2", fontsize=14)
legend(["Class 1", "Class 2", "Class 3"], fontsize=12)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/EX-wine-pca.svg") # hide

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
#
# asian-winstonplot.jl
#
# Note:
#
# With OSX I have had some problems using the Pango as part of the Winston package
# to display titles and text etc.
#
# This seems to be a problem with Homebrew and the following worked for me
# It only needs to be run once,
#
# using Homebrew
# Homebrew.update()
# Homebrew.rm("pango")
# Homebrew.add("pango")
#

import Winston;
wn = Winston;

S0  = 100;      # Spot price
K   = 102;      # Strike price
r   = 0.05;     # Risk free rate
q   = 0.0;      # Dividend yield
v   = 0.2;      # Volatility
tma = 0.25;     # Time to maturity
T = 100;        # Number of time steps
dt = tma/T;     # Time increment
S = zeros(Float64,T)
S[1] = S0;

dW = randn(T)*sqrt(dt);
[ S[t] = S[t-1] * (1 + (r - q - 0.5*v*v)*dt + v*dW[t] + 0.5*v*v*dW[t]*dW[t]) for t=2:T ]
x = linspace(1, T, T);
p = wn.FramedPlot(title = "Random Walk, drift 5%, volatility 20%")
wn.add(p, wn.Curve(x,S,color="red"))
wn.display(p)

max_stars_repo_path: Module 3/Chapter01/asian-winstonplot.jl

max_issues_repo_path: Module 3/Chapter01/asian-winstonplot.jl

max_forks_repo_path: Module 3/Chapter01/asian-winstonplot.jl


text:
module NVML

using ..APIUtils

using ..CUDA

using CEnum

using Libdl

libnvml() = Sys.iswindows() ? :nvml : Symbol("libnvidia-ml.so.1")
has_nvml() = Libdl.dlopen(libnvml(); throw_error=false) !== nothing

# core library
include("libnvml_common.jl")
include("error.jl")
include("libnvml.jl")

# wrappers
include("system.jl")
include("device.jl")

end