text:
const drivers = AG.listdrivers()
const drivermapping = Dict(
    ".shp" => "ESRI Shapefile",
    ".gpkg" => "GPKG",
    ".geojson" => "GeoJSON",
    ".vrt" => "VRT",
    ".csv" => "CSV",
    ".fgb" => "FlatGeobuf",
    ".gml" => "GML",
    ".nc" => "netCDF",
)

# Support "exotic" types, but this needs some piracy
Base.convert(::Type{AG.OGRFieldType}, ft::Type{Bool}) = AG.OFTInteger
Base.convert(::Type{AG.OGRFieldType}, ft::Type{Int16}) = AG.OFTInteger
Base.convert(::Type{AG.OGRFieldType}, ft::Type{Float32}) = AG.OFTReal

subtypes = Dict(
    Bool => AG.OFSTBoolean,
    Int16 => AG.OFSTInt16,
    Float32 => AG.OFSTFloat32,
    )

function AG.setfield!(feature::AG.Feature, i::Integer, value::Int16)
    AG.GDAL.ogr_f_setfieldinteger(feature.ptr, i, value)
    return feature
end

function AG.setfield!(feature::AG.Feature, i::Integer, value::Bool)
    AG.GDAL.ogr_f_setfieldinteger(feature.ptr, i, value)
    return feature
end

function AG.setfield!(feature::AG.Feature, i::Integer, value::Float32)
    AG.GDAL.ogr_f_setfielddouble(feature.ptr, i, value)
    return feature
end

function read(fn::AbstractString; kwargs...)
    ds = AG.read(fn; kwargs...)
    if ArchGDAL.nlayer(ds) > 1
        @warn "This file has multiple layers, you only get the first layer by default now."
    end
    read(ds, 0)
end

function read(fn::AbstractString, layer::Union{Integer,AbstractString}; kwargs...)
    ds = AG.read(fn; kwargs...)
    read(ds, layer)
end

function read(ds, layer)
    table = AG.getlayer(ds, layer)
    if table.ptr == C_NULL
        throw(ArgumentError("Given layer id/name doesn't exist. For reference this is the dataset:\n$ds"))
    end
    df = DataFrame(table)
    "" in names(df) && rename!(df, Dict(Symbol("") => :geom, ))  # needed for now
    df
end

function write(fn::AbstractString, table; layer_name::AbstractString="data", geom_column::Symbol=:geom, crs::GFT.GeoFormat=GFT.EPSG(4326), driver::Union{Nothing,AbstractString}=nothing)
    rows = Tables.rows(table)
    sch = Tables.schema(rows)

    # Geom type can't be inferred from here
    geom_type = AG.getgeomtype(rows[1][geom_column])

    # Find driver
    _, extension = splitext(fn)
    if extension in keys(drivermapping)
        driver = AG.getdriver(drivermapping[extension])
    elseif driver !== nothing
        driver = AG.getdriver(driver)
    else
        error("Couldn't determine driver for $extension. Please provide one of $(keys(drivermapping))")
    end

    # Figure out attributes
    fields = Vector{Tuple{Symbol,DataType}}()
    for (name, type) in zip(sch.names, sch.types)
        if type != AG.IGeometry && name != geom_column
            types = Base.uniontypes(type)
            if length(types) == 1
                push!(fields, (Symbol(name), type))
            elseif length(types) == 2 && Missing in types
                push!(fields, (Symbol(name), types[2]))
            else
                error("Can't convert to GDAL type from $type")
            end
        end
    end
    AG.create(
        fn,
        driver=driver
    ) do ds
        AG.createlayer(
            name=layer_name,
            geom=geom_type,
            spatialref=AG.importCRS(crs)
        ) do layer
            for (name, type) in fields
                AG.createfielddefn(String(name), convert(AG.OGRFieldType, type)) do fd
                    AG.setsubtype!(fd, get(subtypes, type, AG.OFSTNone))
                    AG.addfielddefn!(layer, fd)
                end
            end
            for row in rows
                AG.createfeature(layer) do feature
                    AG.setgeom!(feature, getproperty(row, geom_column))
                    for (name, _) in fields
                        field = getproperty(row, name)
                        if !ismissing(field)
                            AG.setfield!(feature, AG.findfieldindex(feature, name), getproperty(row, name))
                        else
                            AG.GDAL.ogr_f_setfieldnull(feature.ptr, AG.findfieldindex(feature, name))
                        end
                    end
                end
            end
            AG.copy(layer, dataset=ds, name=layer_name)
        end
    end
    fn
end


text:
function get_follicle_from_whisker(x,ind=1)
   x_f = zeros(Float64,length(x))

   for i=1:length(x)
        try
           x_f[i] = x[i][ind]
        catch
            x_f[i] = NaN
        end
    end
    x_f
end

#=
Going to start with 1 as nearest follicle (this is OPPOSITE of janelia)
=#
function correct_follicle(x::Array{T,1},y::Array{T,1},x_f::Real,y_f::Real) where T
    d1 = sqrt((x[1]-x_f)^2 + (y[1]-y_f)^2)
    d2 = sqrt((x[end]-x_f)^2 + (y[end]-y_f)^2)

    if d2 < d1
        reverse!(x)
        reverse!(y)
    end
    nothing
end

function correct_follicle_all(xx::Array,yy::Array,x_f::Real,y_f::Real)

    for i=1:length(xx)
        try
            correct_follicle(xx[i],yy[i],x_f,y_f)
        catch
        end
    end
end

ext: jl


text:
"""
RationalQuadraticKernel([ρ=1.0[,α=2.0]])
The rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y)=(1+ρ²||x−y||²/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
"""
struct RationalQuadraticKernel{T,Tr,Tα<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::SqEuclidean
    α::Tα
    function RationalQuadraticKernel{T,Tr,Tα}(t::Tr,α::Tα) where {T,Tr,Tα<:Real}
        new{T,Tr,Tα}(t,SqEuclidean(),α)
    end
end

function RationalQuadraticKernel(ρ::T₁=1.0,α::T₂=2.0) where {T₁<:Real,T₂<:Real}
    @check_args(RationalQuadraticKernel, α, α > zero(T₂), "α > 1")
    RationalQuadraticKernel{T₁,ScaleTransform{Base.RefValue{T₁}},T₂}(ScaleTransform(ρ),α)
end

function RationalQuadraticKernel(ρ::A,α::T=2.0) where {A<:AbstractVector{<:Real},T<:Real}
    @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
    RationalQuadraticKernel{eltype(A),ScaleTransform{A},T}(ScaleTransform(ρ),α)
end

function RationalQuadraticKernel(t::Tr,α::T=2.0) where {Tr<:Transform,T<:Real}
    @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
    RationalQuadraticKernel{eltype(t),Tr,T}(t,α)
end

params(k::RationalQuadraticKernel) = (params(transform(k)),k.α)

@inline kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²/κ.α)^(-κ.α)

"""
`GammaRationalQuadraticKernel([ρ=1.0[,α=2.0[,γ=2.0]]])`
The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y)=(1+ρ^(2γ)||x−y||^(2γ)/α)^(-α)
```
where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
"""
struct GammaRationalQuadraticKernel{T,Tr,Tα<:Real,Tγ<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::SqEuclidean
    α::Tα
    γ::Tγ
    function GammaRationalQuadraticKernel{T,Tr,Tα,Tγ}(t::Tr,α::Tα,γ::Tγ) where {T,Tr,Tα<:Real,Tγ<:Real}
        new{T,Tr,Tα,Tγ}(t,SqEuclidean(),α,γ)
    end
end

function GammaRationalQuadraticKernel(ρ::T₁=1.0,α::T₂=2.0,γ::T₃=2.0) where {T₁<:Real,T₂<:Real,T₃<:Real}
    @check_args(GammaRationalQuadraticKernel, α, α > one(T₂), "α > 1")
    @check_args(GammaRationalQuadraticKernel, γ, γ >= one(T₂), "γ >= 1")
    GammaRationalQuadraticKernel{T₁,ScaleTransform{Base.RefValue{T₁}},T₂,T₃}(ScaleTransform(ρ),α,γ)
end

function GammaRationalQuadraticKernel(ρ::A,α::T₁=2.0,γ::T₂=2.0) where {A<:AbstractVector{<:Real},T₁<:Real,T₂<:Real}
    @check_args(GammaRationalQuadraticKernel, α, α > one(T₁), "α > 1")
    @check_args(GammaRationalQuadraticKernel, γ, γ >= one(T₂), "γ >= 1")
    GammaRationalQuadraticKernel{eltype(A),ScaleTransform{A},T₁,T₂}(ScaleTransform(ρ),α,γ)
end

function GammaRationalQuadraticKernel(t::Tr,α::T₁=2.0,γ::T₂=2.0) where {Tr<:Transform,T₁<:Real,T₂<:Real}
    @check_args(GammaRationalQuadraticKernel, α, α > one(T₁), "α > 1")
    @check_args(GammaRationalQuadraticKernel, γ, γ >= one(T₂), "γ >= 1")
    GammaRationalQuadraticKernel{eltype(t),Tr,T₁,T₂}(t,α,γ)
end

params(k::GammaRationalQuadraticKernel) = (params(k.transform),k.α,k.γ)

@inline kappa(κ::GammaRationalQuadraticKernel, d²::T) where {T<:Real} = (one(T)+d²^κ.γ/κ.α)^(-κ.α)


text:
#P is the phylogenetic matrix which is M x N matrix, M=num species(configuration), N=num of domains (spins)
using Printf
using LinearAlgebra

##################################################
# Evaluate hamming distance
##################################################

@inline eval_hamm(ai,bi)= ai != bi ? 1 : 0

function hamming_distance(a::AbstractArray, b::AbstractArray)
    if(length(a) != length(b))
        throw( DimensionMismatch("First array has length $(length(a)), while second array has length $(length(b))"))
    end
    if length(a) == length(b)
        dist=0.0
        @simd for I in eachindex(a,b)
            @inbounds ai=a[I]
            @inbounds bi=b[I]
            dist+=eval_hamm(ai,bi)
        end
        dist/=length(a)
    end

    return dist

end

##################################################
# compute weights
##################################################

function compute_weigths(P::Matrix{Int8}, theta::Real)

    M,N = size(P)
    W = ones(M)

    theta = Float64(theta)
    @printf("%s%0.2f ",  "theta=",theta)

    if theta == 0
        @printf("%s%d%s%d%s%d\n","M=",M," N=",N," Meff=",M)
        return W, Float64(M)
    end

    ZZ = Vector{Int8}[vec(P[i,:]) for i=1:M]
    @simd for i = 1:M
        Zi = ZZ[i]
        @simd for j = (i+1):M
            Zj=ZZ[j]
            #if sequence has distance < threshold -> increase weight
            if hamming_distance(Zi,Zj)<theta
                @inbounds W[i] += 1
                @inbounds W[j] += 1
            end
        end
        #weight
    end
    for i = 1:M
        W[i] = 1 / W[i]
    end

    Meff = sum(W)

    @printf("%s%d%s%d%s%0.2f\n","M=",M," N=",N," Meff=",Meff)

    return W, Meff
end

##################################################
# compute frequencies
##################################################

# W is the vector of weights
function compute_freqs(P::Matrix{Int8}, W::Vector{Float64}, Meff::Float64)
    M,N = size(P)
    ZZ = Vector{Int8}[vec(P[:,i]) for i=1:N]

    Pij = zeros(N,N)
    Pi = zeros(N)
    for i = 1:N
        Zi = ZZ[i]
        for j = i:N
            Zj = ZZ[j]
            for k = 1:M
                Pij[i,j] += W[k]*Zi[k]*Zj[k]
            end
        end
    end

    for i = 1:N
        Pij[i,i] /= Meff
        for j = (i+1):N
            Pij[i,j] /= Meff
            Pij[j,i] = Pij[i,j]
        end
    end

     Pi = diag(Pij)
    return Pi,Pij
end

function compute_new_frequencies(P::Matrix{Int8},theta::Real)

    W,Meff = compute_weigths(P,theta)
    Pi_t,Pij_t = compute_freqs(P,W,Meff)

    return Pi_t,Pij_t,Meff,W
end

##################################################
# add pseudocount
##################################################
function add_pseudocount(Pi_true::Vector{Float64},Pij_true::Matrix{Float64},Meff::Float64,pc::Float64)

    N = length(Pi_true)
    Pij = (1-pc) * Pij_true .+ pc/4
    for i = 1:N
        Pij[i,i] += pc/2
    end

    Pi = ( 1 - pc ) * Pi_true .+ pc/2
    return Pi, Pij
end

##################################################
# compute empirical correlation matrix C
##################################################
compute_C(Pi::Vector{Float64}, Pij::Matrix{Float64}) = Pij - Pi * Pi'

##################################################
# APC correction (minimize backgroung influence like phylogeny and site entropy)
##################################################
function correct_APC(S::Matrix)
    N = size(S,1)
    Si = sum(S,dims=1)
    Sj = sum(S,dims=2)
    Sa = sum(S) * (1 - 1/N)

    S -= (Sj * Si) / Sa
    return S
end

max_stars_repo_name: clayfos/PhyDCA

max_issues_repo_name: clayfos/PhyDCA

max_forks_repo_name: clayfos/PhyDCA


text:
using Pkg
for p in ("Knet","ArgParse")
    haskey(Pkg.installed(),p) || Pkg.add(p)
end

"""

julia bilstm-tagger.jl # to use with default options on CPU
julia bilstm-tagger.jl --usegpu # to use with default options on GPU
julia bilstm-tagger.jl -h # to see all options with default values

This example implements a named entity tagger built on top of a BiLSTM
neural network similar to the model defined in 'Bidirectional LSTM-CRF Models
for Sequence Tagging', Zhiheng Huang, Wei Xu, Kai Yu, arXiv technical report
1508.01991, 2015. Originally, this model implemented for dynet-benchmarks.

* Paper url: https://arxiv.org/pdf/1508.01991.pdf
* DyNet report: https://arxiv.org/abs/1701.03980
* Benchmark repo: https://github.com/neulab/dynet-benchmark

"""
module Tagger
using Knet, ArgParse, Dates, Random, Printf

include(Knet.dir("data","wikiner.jl"))
t00 = now()

function main(args)
    s = ArgParseSettings()
    s.description = "Bidirectional LSTM Tagger in Knet"

    @add_arg_table s begin
        ("--usegpu"; action=:store_true; help="use GPU or not")
        ("--embed"; arg_type=Int; default=128; help="word embedding size")
        ("--hidden"; arg_type=Int; default=50; help="LSTM hidden size")
        ("--mlp"; arg_type=Int; default=32; help="MLP size")
        ("--timeout"; arg_type=Int; default=600; help="max timeout (in seconds)")
        ("--epochs"; arg_type=Int; default=100; help="number of training epochs")
        ("--minoccur"; arg_type=Int; default=6; help="word min occurence limit")
        ("--report"; arg_type=Int; default=500; help="report period in iters")
        ("--valid"; arg_type=Int; default=10000; help="valid period in iters")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:seed] > 0 && Knet.setseed(o[:seed])
    atype = o[:atype] = !o[:usegpu] ? Array{Float32} : KnetArray{Float32}
    datadir = abspath(joinpath(@__DIR__, "../data/tags"))
    datadir = isdir(datadir) ? datadir : WIKINER_DIR

    # load WikiNER data
    data = WikiNERData(datadir, o[:minoccur])

    # build model
    nwords = length(data.w2i); ntags = data.ntags
    w, srnn = initweights(
        atype, o[:hidden], nwords, ntags, o[:mlp], o[:embed], o[:usegpu])
    opt = optimizers(w, Adam)

    # train bilstm tagger
    nwords = data.nwords; ntags = data.ntags
    println("nwords=$nwords, ntags=$ntags"); flush(stdout)
    println("startup time: ", Int((now()-t00).value)*0.001); flush(stdout)
    t0 = now()
    all_time = dev_time = all_tagged = this_tagged = this_loss = 0
    o[:timeout] = o[:timeout] <= 0 ? Inf : o[:timeout]
    for epoch = 1:o[:epochs]
        shuffle!(data.trn)
        for k = 1:length(data.trn)
            iter = (epoch-1)*length(data.trn) + k
            if o[:report] > 0 && iter % o[:report] == 0
                @printf("%f\n", this_loss/this_tagged); flush(stdout)
                all_tagged += this_tagged
                this_loss = this_tagged = 0
                all_time = Int((now()-t0).value)*0.001
            end

            if all_time > o[:timeout] || o[:valid] > 0 && iter % o[:valid] == 0
                dev_start = now()
                good_sent = bad_sent = good = bad = 0.0
                for sent in data.dev
                    seq = make_input(sent, data.w2i)
                    nwords = length(sent)
                    ypred,_ = predict(w, seq, srnn)
                    ypred = map(
                        x->data.i2t[x], mapslices(argmax,Array(ypred),dims=1))
                    ygold = map(x -> x[2], sent)
                    same = true
                    for (y1,y2) in zip(ypred, ygold)
                        if y1 == y2
                            good += 1
                        else
                            bad += 1
                            same = false
                        end
                    end
                    if same
                        good_sent += 1
                    else
                        bad_sent += 1
                    end
                end
                dev_time += Int((now()-dev_start).value)*0.001
                train_time = Int((now()-t0).value)*0.001-dev_time

                @printf(
                    "tag_acc=%.4f, sent_acc=%.4f, time=%.4f, word_per_sec=%.4f\n",
                    good/(good+bad), good_sent/(good_sent+bad_sent), train_time,
                    all_tagged/train_time); flush(stdout)

                all_time > o[:timeout] && return
            end

            # train on minibatch
            x = make_input(data.trn[k], data.w2i)
            y = make_output(data.trn[k], data.t2i)
            batch_loss = train!(w,x,y,srnn,opt)
            this_loss += batch_loss
            this_tagged += length(data.trn[k])
        end
        @printf("epoch %d finished\n", epoch-1); flush(stdout)
    end
end

function make_input(sample, w2i)
    nwords = length(sample)
    x = map(i->get(w2i, sample[i][1], w2i[UNK]), [1:nwords...])
    x = reshape(x,1,length(x))
    x = convert(Array{Int64}, x)
end

function make_output(sample, t2i)
    nwords = length(sample)
    y = map(i->t2i[sample[i][2]], [1:nwords...])
    y = reshape(y,1,length(y))
    y = convert(Array{Int64}, y)
end

# w[1]   => weight/bias params for forward LSTM network
# w[2:5] => weight/bias params for MLP+softmax network
# w[6]   => word embeddings
function initweights(atype, hidden, words, tags, embed, mlp, usegpu, winit=0.01)
    w = Array{Any}(undef,6)
    input = embed
    srnn, wrnn = rnninit(input, hidden; bidirectional=true, usegpu=usegpu)
    w[1] = wrnn
    w[2] = convert(atype, winit*randn(mlp, 2*hidden))
    w[3] = convert(atype, zeros(mlp, 1))
    w[4] = convert(atype, winit*randn(tags, mlp))
    w[5] = convert(atype, winit*randn(tags, 1))
    w[6] = convert(atype, winit*randn(embed, words))
    return w, srnn
end

# loss function
function loss(w, x, ygold, srnn, h=nothing, c=nothing)
    py, _ = predict(w,x,srnn,h,c)
    return nll(py,ygold)
end

function predict(ws,xs,srnn,hx=nothing,cx=nothing)
    wx = ws[6]
    r = srnn; wr = ws[1]
    wmlp = ws[2]; bmlp = ws[3];
    wy = ws[4]; by = ws[5]
    x = wx[:,xs]
    y, hy, cy = rnnforw(r,wr,x)
    y2 = reshape(y,size(y,1),size(y,2)*size(y,3))
    y3 = wmlp * y2 .+ bmlp
    return wy*y3.+by, hy, cy
end

lossgradient = gradloss(loss)

function train!(w,x,y,srnn,opt,h=nothing,c=nothing)
    gloss, lossval = lossgradient(w, x, y, srnn)
    update!(w,gloss,opt)
    return lossval*size(x,2)
end

splitdir(PROGRAM_FILE)[end] == "bilstm-tagger.jl" && main(ARGS)
end # module


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license
using FFTW
using FFTW: fftw_provider
using AbstractFFTs: Plan, plan_inv
using Test
using LinearAlgebra

@info "FFTW provider: $(fftw_provider)"

@testset "colmajorstrides" begin
    @test @inferred(FFTW.colmajorstrides((3, 5, 7))) == (1, 3, 15)
end

m4 = [16.    2     3    13;
      5    11    10     8;
      9     7     6    12;
      4    14    15     1]

true_fft_m4 = [
    34.            34.            34.            34.;
    7. - 1.0im  -5. + 3.0im  -3. + 5.0im   1. - 7.0im;
    16.           -16.           -16.            16.;
    7. + 1.0im  -5. - 3.0im  -3. - 5.0im   1. + 7.0im ]

true_fftn_m4 = [
    136.        0          0         0 ;
    0.       20          8 + 8im   0 - 12im ;
    0.       32 + 32im   0        32 - 32im ;
    0.        0 + 12im   8 - 8im  20 ]

true_fftd2_m4 = [
    34.   13 + 11im    4   13 - 11im ;
    34.   -5 -  3im   -4   -5 +  3im ;
    34.    3 +  5im   -4    3 -  5im ;
    34.  -11 - 13im    4  -11 + 13im ]

b = rand(17,14)
b[3:6,9:12] = m4
sm4 = view(b,3:6,9:12)

m3d = map(Float32,copy(reshape(1:5*3*2, 5, 3, 2)))
true_fftd3_m3d = Array{Float32}(undef, 5, 3, 2)
true_fftd3_m3d[:,:,1] = 17:2:45
true_fftd3_m3d[:,:,2] .= -15

@testset "fft" begin
    a = rand(8) + im*rand(8)
    @test norm(ifft(fft(a)) - a) < 1e-8
    @test norm(ifft(fft(a,1),1) - a) < 1e-8
    @test norm(ifft(fft(a,[1]),[1]) - a) < 1e-8
    @test norm(ifft(fft(a,(1,)),(1,)) - a) < 1e-8
    a = rand(-10:10, 8) + im*rand(-10:10, 8)
    @test norm(ifft(fft(a)) - a) < 1e-8

    # use invoke to force usage of CTPlan versions even if FFTW is present
    for A in (Array,SubArray)
        for f in (:fft,:ifft,:plan_fft,:plan_ifft)
            f_ = Symbol(f, "_")
            @eval begin
                $f_(x::$A{T,N}) where {T,N} = invoke($f, Tuple{AbstractArray{T,N}}, x)
                $f_(x::$A{T,N},r::R) where {T,N,R} = invoke($f,Tuple{AbstractArray{T,N},R},x,r)
            end
        end
    end

    for (f,fi,pf,pfi) in ((fft,ifft,plan_fft,plan_ifft),
                          (fft_,ifft_,plan_fft_,plan_ifft_))
        pm4 = pf(m4,1)

        fft_m4 = f(m4,1)
        fftd2_m4 = f(m4,2)
        ifft_fft_m4 = fi(f(m4,1),1)
        fftn_m4 = f(m4)
        ifftn_fftn_m4 = fi(f(m4))

        fft!_m4 = complex(m4); fft!(fft!_m4,1)
        fft!d2_m4 = complex(m4); fft!(fft!d2_m4,2)
        ifft!_fft_m4 = f(m4,1); ifft!(ifft!_fft_m4,1)
        fft!n_m4 = complex(m4); fft!(fft!n_m4)
        ifft!n_fftn_m4 = f(m4); ifft!(ifft!n_fftn_m4)

        pfft_m4 = pf(m4,1)*m4
        pfftd2_m4 = pf(m4,2)*m4
        pifft_fft_m4 = pfi(fft_m4,1)*fft_m4
        pfftn_m4 = pf(m4)*m4
        pifftn_fftn_m4 = pfi(fftn_m4)*fftn_m4

        pfft!_m4 = complex(m4); plan_fft!(pfft!_m4,1)*pfft!_m4
        pfft!d2_m4 = complex(m4); plan_fft!(pfft!d2_m4,2)*pfft!d2_m4
        pifft!_fft_m4 = f(m4,1); plan_ifft!(pifft!_fft_m4,1)*pifft!_fft_m4
        pfft!n_m4 = complex(m4); plan_fft!(pfft!n_m4)*pfft!n_m4
        pifft!n_fftn_m4 = f(m4); plan_ifft!(pifft!n_fftn_m4)*pifft!n_fftn_m4

        sfftn_m4 = f(sm4)
        psfftn_m4 = pf(sm4)*sm4
        sfft!n_b = map(Complex{Float64},b)
        sfft!n_m4 = view(sfft!n_b,3:6,9:12); fft!(sfft!n_m4)
        psfft!n_b = map(Complex{Float64},b)
        psfft!n_m4 = view(psfft!n_b,3:6,9:12); plan_fft!(psfft!n_m4)*psfft!n_m4

        for i = 1:length(m4)
            @test fft_m4[i] ≈ true_fft_m4[i]
            @test fftd2_m4[i] ≈ true_fftd2_m4[i]
            @test ifft_fft_m4[i] ≈ m4[i]
            @test fftn_m4[i] ≈ true_fftn_m4[i]
            @test ifftn_fftn_m4[i] ≈ m4[i]

            @test fft!_m4[i] ≈ true_fft_m4[i]
            @test fft!d2_m4[i] ≈ true_fftd2_m4[i]
            @test ifft!_fft_m4[i] ≈ m4[i]
            @test fft!n_m4[i] ≈ true_fftn_m4[i]
            @test ifft!n_fftn_m4[i] ≈ m4[i]

            @test pfft_m4[i] ≈ true_fft_m4[i]
            @test pfftd2_m4[i] ≈ true_fftd2_m4[i]
            @test pifft_fft_m4[i] ≈ m4[i]
            @test pfftn_m4[i] ≈ true_fftn_m4[i]
            @test pifftn_fftn_m4[i] ≈ m4[i]

            @test pfft!_m4[i] ≈ true_fft_m4[i]
            @test pfft!d2_m4[i] ≈ true_fftd2_m4[i]
            @test pifft!_fft_m4[i] ≈ m4[i]
            @test pfft!n_m4[i] ≈ true_fftn_m4[i]
            @test pifft!n_fftn_m4[i] ≈ m4[i]

            @test sfftn_m4[i] ≈ true_fftn_m4[i]
            @test sfft!n_m4[i] ≈ true_fftn_m4[i]
            @test psfftn_m4[i] ≈ true_fftn_m4[i]
            @test psfft!n_m4[i] ≈ true_fftn_m4[i]
        end

        ifft!(sfft!n_m4)
        plan_ifft!(psfft!n_m4)*psfft!n_m4
        @test norm(sfft!n_m4 - m4) < 1e-8
        @test norm(psfft!n_m4 - m4) < 1e-8

        # The following capabilities are FFTW only.
        # They are not available in MKL, and hence do not test them.
        if fftw_provider == "fftw"
            @testset "FFTW-specific" begin
                ifft3_fft3_m3d = fi(f(m3d))

                fftd3_m3d = f(m3d,3)
                ifftd3_fftd3_m3d = fi(fftd3_m3d,3)

                fft!d3_m3d = complex(m3d); fft!(fft!d3_m3d,3)
                ifft!d3_fftd3_m3d = copy(fft!d3_m3d); ifft!(ifft!d3_fftd3_m3d,3)

                pfftd3_m3d = pf(m3d,3)*m3d
                pifftd3_fftd3_m3d = pfi(fftd3_m3d,3)*fftd3_m3d

                pfft!d3_m3d = complex(m3d); plan_fft!(pfft!d3_m3d,3)*pfft!d3_m3d
                pifft!d3_fftd3_m3d = copy(fft!d3_m3d); plan_ifft!(pifft!d3_fftd3_m3d,3)*pifft!d3_fftd3_m3d

                @test isa(fftd3_m3d, Array{Complex{Float32},3})
                @test isa(ifftd3_fftd3_m3d, Array{Complex{Float32},3})
                @test isa(fft!d3_m3d, Array{Complex{Float32},3})
                @test isa(ifft!d3_fftd3_m3d, Array{Complex{Float32},3})
                @test isa(pfftd3_m3d, Array{Complex{Float32},3})
                @test isa(pifftd3_fftd3_m3d, Array{Complex{Float32},3})
                @test isa(pfft!d3_m3d, Array{Complex{Float32},3})
                @test isa(pifft!d3_fftd3_m3d, Array{Complex{Float32},3})

                for i = 1:length(m3d)
                    @test fftd3_m3d[i] ≈ true_fftd3_m3d[i]
                    @test ifftd3_fftd3_m3d[i] ≈ m3d[i]
                    @test ifft3_fft3_m3d[i] ≈ m3d[i]

                    @test fft!d3_m3d[i] ≈ true_fftd3_m3d[i]
                    @test ifft!d3_fftd3_m3d[i] ≈ m3d[i]

                    @test pfftd3_m3d[i] ≈ true_fftd3_m3d[i]
                    @test pifftd3_fftd3_m3d[i] ≈ m3d[i]
                    @test pfft!d3_m3d[i] ≈ true_fftd3_m3d[i]
                    @test pifft!d3_fftd3_m3d[i] ≈ m3d[i]
                end
            end
        end  # if fftw_provider == "fftw"

    end
end

@testset "rfft/rfftn" begin
    rfft_m4 = rfft(m4,1)
    rfftd2_m4 = rfft(m4,2)
    rfftn_m4 = rfft(m4)

    prfft_m4 = plan_rfft(m4,1)*m4
    prfftd2_m4 = plan_rfft(m4,2)*m4
    prfftn_m4 = plan_rfft(m4)*m4

    srfftn_m4 = rfft(sm4)
    psrfftn_m4 = plan_rfft(sm4)*sm4

    for i = 1:3, j = 1:4
        @test rfft_m4[i,j] ≈ true_fft_m4[i,j]
        @test rfftd2_m4[j,i] ≈ true_fftd2_m4[j,i]
        @test rfftn_m4[i,j] ≈ true_fftn_m4[i,j]

        @test prfft_m4[i,j] ≈ true_fft_m4[i,j]
        @test prfftd2_m4[j,i] ≈ true_fftd2_m4[j,i]
        @test prfftn_m4[i,j] ≈ true_fftn_m4[i,j]

        @test srfftn_m4[i,j] ≈ true_fftn_m4[i,j]
        @test psrfftn_m4[i,j] ≈ true_fftn_m4[i,j]
    end

    irfft_rfft_m4 = irfft(rfft_m4,size(m4,1),1)
    irfft_rfftd2_m4 = irfft(rfftd2_m4,size(m4,2),2)
    irfftn_rfftn_m4 = irfft(rfftn_m4,size(m4,1))

    pirfft_rfft_m4 = plan_irfft(rfft_m4,size(m4,1),1)*rfft_m4
    pirfft_rfftd2_m4 = plan_irfft(rfftd2_m4,size(m4,2),2)*rfftd2_m4
    pirfftn_rfftn_m4 = plan_irfft(rfftn_m4,size(m4,1))*rfftn_m4

    for i = 1:length(m4)
        @test irfft_rfft_m4[i] ≈ m4[i]
        @test irfft_rfftd2_m4[i] ≈ m4[i]
        @test irfftn_rfftn_m4[i] ≈ m4[i]

        @test pirfft_rfft_m4[i] ≈ m4[i]
        @test pirfft_rfftd2_m4[i] ≈ m4[i]
        @test pirfftn_rfftn_m4[i] ≈ m4[i]
    end

    @testset "with preallocated array" begin
        prfft_m4_prealloc = zero(prfft_m4)
        prfftd2_m4_prealloc = zero(prfftd2_m4)
        prfftn_m4_prealloc = zero(prfftn_m4)

        pirfft_rfft_m4_prealloc = zero(pirfft_rfft_m4)
        pirfft_rfftd2_m4_prealloc = zero(pirfft_rfftd2_m4)
        pirfftn_rfftn_m4_prealloc = zero(pirfftn_rfftn_m4)

        mul!(prfft_m4_prealloc,plan_rfft(m4,1),m4)
        mul!(prfftd2_m4_prealloc,plan_rfft(m4,2),m4)
        mul!(prfftn_m4_prealloc,plan_rfft(m4),m4)

        for i = 1:3, j = 1:4
            @test prfft_m4_prealloc[i,j] ≈ true_fft_m4[i,j]
            @test prfftd2_m4_prealloc[j,i] ≈ true_fftd2_m4[j,i]
            @test prfftn_m4_prealloc[i,j] ≈ true_fftn_m4[i,j]
        end

        mul!(pirfft_rfft_m4_prealloc,plan_irfft(rfft_m4,size(m4,1),1),prfft_m4_prealloc)
        mul!(pirfft_rfftd2_m4_prealloc,plan_irfft(rfftd2_m4,size(m4,2),2),prfftd2_m4_prealloc)
        mul!(pirfftn_rfftn_m4_prealloc,plan_irfft(rfftn_m4,size(m4,1)),prfftn_m4_prealloc)

        for i = 1:length(m4)
            @test pirfft_rfft_m4_prealloc[i] ≈ m4[i]
            @test pirfft_rfftd2_m4_prealloc[i] ≈ m4[i]
            @test pirfftn_rfftn_m4_prealloc[i] ≈ m4[i]
        end
    end

    if fftw_provider == "fftw"
        @testset "FFTW-specific" begin
            rfftn_m3d = rfft(m3d)
            rfftd3_m3d = rfft(m3d,3)
            @test size(rfftd3_m3d) == size(true_fftd3_m3d)
            irfft_rfftd3_m3d = irfft(rfftd3_m3d,size(m3d,3),3)
            irfftn_rfftn_m3d = irfft(rfftn_m3d,size(m3d,1))
            for i = 1:length(m3d)
                @test rfftd3_m3d[i] ≈ true_fftd3_m3d[i]
                @test irfft_rfftd3_m3d[i] ≈ m3d[i]
                @test irfftn_rfftn_m3d[i] ≈ m3d[i]
            end

            fftn_m3d = fft(m3d)
            @test size(fftn_m3d) == (5,3,2)
            rfftn_m3d = rfft(m3d)
            @test size(rfftn_m3d) == (3,3,2)
            for i = 1:3, j = 1:3, k = 1:2
                @test rfftn_m3d[i,j,k] ≈ fftn_m3d[i,j,k]
            end
        end # !mkl
    end
end

# FFT self-test algorithm (for unscaled 1d forward FFTs):
#   Funda Ergün, "Testing multivariate linear functions: Overcoming
#   the generator bottleneck," Proc. 27th ACM Symposium on the Theory
#   of Computing, pp. 407-416 (1995).
# Check linearity, impulse-response, and time-shift properties.
function fft_test(p::Plan{T}, ntrials=4,
                  tol=1e5 * eps(real(T))) where T<:Complex
    ndims(p) == 1 || throw(ArgumentError("not a 1d FFT"))
    n = length(p)
    twopi_i = (-2 * convert(real(T), π)/n * (0:n-1)) * im
    for trial = 1:ntrials
        # linearity:
        x = rand(T, n)
        y = rand(T, n)
        α = rand(T)
        β = rand(T)
        X = p * (α*x + β*y)
        err = norm(α * (p*x) + β * (p*y) - X, Inf) / norm(X, Inf)
        err <= tol || error("linearity error $err in $p")

        # impulse-response:
        z = zeros(T, n)
        i = rand(0:n-1)
        z[i+1] = 1
        X = exp.(twopi_i*i)
        err = norm(p*z - X, Inf) / norm(X, Inf)
        err <= tol || error("impulse-response error $err in $p")

        # time-shift:
        if n > 1
            s = rand(1:n-1)
            X = (p*x).*exp.(twopi_i*s)
            err = norm(p*circshift(x,s) - X, Inf) / norm(X, Inf)
            err <= tol || error("time-shift error $err in $p")
        end
    end
end

@testset "Self-test: $(T)" for T in (Complex{Float32}, Complex{Float64})
    for n in [1:100; 121; 143; 1000; 1024; 1031; 2000; 2048]
        x = zeros(T, n)
        fft_test(plan_fft(x))
        fft_test(plan_fft_(x))
    end
end

@testset "Inversion, scaling, and pre-allocated variants: $(T)" for T in (Complex{Float32}, Complex{Float64})
    for x in (T[1:100;], copy(reshape(T[1:200;], 20,10)))
        y = similar(x)
        for planner in (plan_fft, plan_fft_, plan_ifft, plan_ifft_)
            p = planner(x)
            pi = inv(p)
            p3 = 3*p
            p3i = inv(p3)
            @test eltype(p) == eltype(pi) == eltype(p3) == eltype(p3i) == T
            @test norm(x - p3i * (p * 3x)) < eps(real(T)) * 10000
            @test norm(3x - pi * (p3 * x)) < eps(real(T)) * 10000
            mul!(y, p, x)
            @test y == p * x
            ldiv!(y, p, x)
            @test y == p \ x
        end
    end
end

let
    plan32 = plan_fft([1.0:2048.0;])
    plan64 = plan_fft([1f0:2048f0;])
    FFTW.flops(plan32)
    FFTW.flops(plan64)
end

@testset "Base Julia issue #9772, with size $(size(x))" for x in (randn(10),randn(10,12))
    z = complex(x)
    y = rfft(x)
    @inferred rfft(x)

    # See Julia issue #23063
    if ndims(x) == 2
        @test_broken @inferred brfft(x,18)
    end

    @inferred brfft(y,10)
    for f in (plan_bfft!, plan_fft!, plan_ifft!,
              plan_bfft, plan_fft, plan_ifft,
              fft, bfft, fft_, ifft)
        p = @inferred f(z)
        if isa(p, Plan)
            @inferred plan_inv(p)
        end
    end
    for f in (plan_bfft, plan_fft, plan_ifft,
              plan_rfft, fft, bfft, fft_, ifft)
        # More of #23063 (why does plan_rfft work and the others don't)?
        if ndims(x) == 2 && f != plan_rfft
            @test_broken @inferred f(x)
            @test_broken @inferred plan_inv(f(x))
            continue
        end
        p = @inferred f(x)
        if isa(p, Plan)
            @inferred plan_inv(p)
        end
    end
    # note: inference doesn't work for plan_fft_ since the
    #       algorithm steps are included in the CTPlan type
end

@testset "Base Julia issue #17896" begin
    a = rand(5)
    @test  fft(a) ==  fft(view(a,:)) ==  fft(view(a, 1:5)) ==  fft(view(a, [1:5;]))
    @test rfft(a) == rfft(view(a,:)) == rfft(view(a, 1:5)) == rfft(view(a, [1:5;]))
    a16 = convert(Vector{Float16}, a)
    @test  fft(a16) ==  fft(view(a16,:)) ==  fft(view(a16, 1:5)) ==  fft(view(a16, [1:5;]))
    @test rfft(a16) == rfft(view(a16,:)) == rfft(view(a16, 1:5)) == rfft(view(a16, [1:5;]))
end

if fftw_provider == "fftw"
    @testset "Discrete cosine transform (DCT) tests" begin
        a = rand(8,11) + im*rand(8,11)
        @test norm(idct(dct(a)) - a) < 1e-8

        X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4)
        Y = rand(17,14)
        Y[3:5,9:12] = X
        sX = view(Y,3:5,9:12)

        true_Xdct = [  13.856406460551018  -3.863239728836245   2.886751345948129  -0.274551994240164; -2.828427124746190  -2.184015211898548  -4.949747468305834   3.966116180118245; 4.898979485566356  -0.194137576915510  -2.857738033247041   2.731723009609389 ]

        true_Xdct_1 = [    5.773502691896258   4.618802153517007   6.350852961085884  10.969655114602890; -4.242640687119286  -2.121320343559643   4.242640687119286  -3.535533905932738; 1.632993161855452   2.041241452319315   5.715476066494083   0.408248290463863 ]

        true_Xdct_2 = [    8.  -3.854030797826254  -3.0  3.761176226848022;
                           4.0  -2.071929829606556   4.0  -2.388955165168770; 12.  -0.765366864730179   4.0  -1.847759065022573 ]

        Xdct = dct(X)
        Xdct! = float(X); dct!(Xdct!)
        Xdct_1 = dct(X,1)
        Xdct!_1 = float(X); dct!(Xdct!_1,1)
        Xdct_2 = dct(X,2)
        Xdct!_2 = float(X); dct!(Xdct!_2,2)

        Xidct = idct(true_Xdct)
        Xidct! = copy(true_Xdct); idct!(Xidct!)
        Xidct_1 = idct(true_Xdct_1,1)
        Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1)
        Xidct_2 = idct(true_Xdct_2,2)
        Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2)

        pXdct = plan_dct(X)*(X)
        pXdct! = float(X); plan_dct!(pXdct!)*(pXdct!)
        pXdct_1 = plan_dct(X,1)*(X)
        pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)*(pXdct!_1)
        pXdct_2 = plan_dct(X,2)*(X)
        pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)*(pXdct!_2)

        pXidct = plan_idct(true_Xdct)*(true_Xdct)
        pXidct! = copy(true_Xdct); plan_idct!(pXidct!)*(pXidct!)
        pXidct_1 = plan_idct(true_Xdct_1,1)*(true_Xdct_1)
        pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)*(pXidct!_1)
        pXidct_2 = plan_idct(true_Xdct_2,2)*(true_Xdct_2)
        pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)*(pXidct!_2)

        sXdct = dct(sX)
        psXdct = plan_dct(sX)*(sX)
        sYdct! = copy(Y); sXdct! = view(sYdct!,3:5,9:12); dct!(sXdct!)
        psYdct! = copy(Y); psXdct! = view(psYdct!,3:5,9:12); plan_dct!(psXdct!)*(psXdct!)

        for i = 1:length(X)
            @test Xdct[i] ≈ true_Xdct[i]
            @test Xdct![i] ≈ true_Xdct[i]
            @test Xdct_1[i] ≈ true_Xdct_1[i]
            @test Xdct!_1[i] ≈ true_Xdct_1[i]
            @test Xdct_2[i] ≈ true_Xdct_2[i]
            @test Xdct!_2[i] ≈ true_Xdct_2[i]

            @test pXdct[i] ≈ true_Xdct[i]
            @test pXdct![i] ≈ true_Xdct[i]
            @test pXdct_1[i] ≈ true_Xdct_1[i]
            @test pXdct!_1[i] ≈ true_Xdct_1[i]
            @test pXdct_2[i] ≈ true_Xdct_2[i]
            @test pXdct!_2[i] ≈ true_Xdct_2[i]

            @test Xidct[i] ≈ X[i]
            @test Xidct![i] ≈ X[i]
            @test Xidct_1[i] ≈ X[i]
            @test Xidct!_1[i] ≈ X[i]
            @test Xidct_2[i] ≈ X[i]
            @test Xidct!_2[i] ≈ X[i]

            @test pXidct[i] ≈ X[i]
            @test pXidct![i] ≈ X[i]
            @test pXidct_1[i] ≈ X[i]
            @test pXidct!_1[i] ≈ X[i]
            @test pXidct_2[i] ≈ X[i]
            @test pXidct!_2[i] ≈ X[i]

            @test sXdct[i] ≈ true_Xdct[i]
            @test psXdct[i] ≈ true_Xdct[i]
            @test sXdct![i] ≈ true_Xdct[i]
            @test psXdct![i] ≈ true_Xdct[i]
        end
    end
end # fftw_provider == "fftw"

@testset "UNALIGNED flag" begin
    A = rand(Float32, 35)
    Ac = rand(Complex{Float32}, 35)
    Y = Array{Complex{Float32}}(undef, 20)
    Yc = Array{Complex{Float32}}(undef, 35)
    planr = plan_rfft(Array{Float32}(undef, 32), flags=FFTW.UNALIGNED)
    planc = plan_fft(Array{Complex{Float32}}(undef, 32), flags=FFTW.UNALIGNED)
    for ioff in 0:3
        ii = 1+ioff:32+ioff
        @test planr * view(A, ii) ≈ planr * A[ii] ≈ rfft(view(A, ii)) ≈ rfft(A[ii])
        @test planc * view(Ac, ii) ≈ planc * Ac[ii] ≈ fft(view(Ac, ii)) ≈ fft(Ac[ii])
        for ooff in 0:3
            io = 1+ooff:17+ooff
            FFTW.mul!(view(Y, io), planr, view(A, ii))
            @test Y[io] ≈ rfft(A[ii])
            FFTW.mul!(view(Y, io), planr, A[ii])
            @test Y[io] ≈ rfft(A[ii])
            io = 1+ooff:32+ooff
            FFTW.mul!(view(Yc, io), planc, view(Ac, ii))
            @test Yc[io] ≈ fft(Ac[ii])
            FFTW.mul!(view(Yc, io), planc, Ac[ii])
            @test Yc[io] ≈ fft(Ac[ii])
        end
    end

    # check whether FFTW on this architecture has nontrivial alignment requirements
    nontrivial_alignment = FFTW.fftw_provider == "fftw" && ccall((:fftwf_alignment_of, FFTW.libfftw3f[]), Int32, (Int,), 8) != 0
    if nontrivial_alignment
        @test_throws ArgumentError plan_rfft(Array{Float32}(undef, 32)) * view(A, 2:33)
        @test_throws ArgumentError plan_fft(Array{Complex{Float32}}(undef, 32)) * view(Ac, 2:33)
    end
end

@testset "Base Julia issue #19892" begin
    a = randn(10^5,1)
    FFTW.set_num_threads(1)
    p1 = plan_rfft(a, flags=FFTW.ESTIMATE)
    FFTW.set_num_threads(2)
    p2 = plan_rfft(a, flags=FFTW.ESTIMATE)
    @test p1*a ≈ p2*a
    # make sure threads are actually being used for p2
    # (tests #21163).
    if FFTW.has_sprint_plan
        @test !occursin("dft-thr", string(p1))
        @test occursin("dft-thr", string(p2))
    end
end


text:
#==
All differentials need to define + and *.
That happens here.

We just use @eval to define all the combinations for AbstractDifferential
subtypes, as we know the full set that might be encountered.
Thus we can avoid any ambiguities.

Notice:
    The precedence goes:
    `DoesNotExist, Zero, One, AbstractThunk, Composite, Any`
    Thus each of the @eval loops create most definitions of + and *
    defines the combination this type with all types of  lower precidence.
    This means each eval loops is 1 item smaller than the previous.
==#

Base.:+(::DoesNotExist, ::DoesNotExist) = DoesNotExist()
Base.:*(::DoesNotExist, ::DoesNotExist) = DoesNotExist()
for T in (:One, :AbstractThunk, :Composite, :Any)
    @eval Base.:+(::DoesNotExist, b::$T) = b
    @eval Base.:+(a::$T, ::DoesNotExist) = a

    @eval Base.:*(::DoesNotExist, ::$T) = DoesNotExist()
    @eval Base.:*(::$T, ::DoesNotExist) = DoesNotExist()
end
# `DoesNotExist` and `Zero` have special relationship,
# DoesNotExist wins add, Zero wins *. This is (in theory) to allow `*` to be used for
# selecting things.
Base.:+(::DoesNotExist, ::Zero) = DoesNotExist()
Base.:+(::Zero, ::DoesNotExist) = DoesNotExist()
Base.:*(::DoesNotExist, ::Zero) = Zero()
Base.:*(::Zero, ::DoesNotExist) = Zero()

Base.muladd(::Zero, x, y) = y
Base.muladd(x, ::Zero, y) = y
Base.muladd(x, y, ::Zero) = x*y

Base.muladd(::Zero, ::Zero, y) = y
Base.muladd(x, ::Zero, ::Zero) = Zero()
Base.muladd(::Zero, x, ::Zero) = Zero()

Base.muladd(::Zero, ::Zero, ::Zero) = Zero()

Base.:+(::Zero, ::Zero) = Zero()
Base.:*(::Zero, ::Zero) = Zero()
for T in (:One, :AbstractThunk, :Composite, :Any)
    @eval Base.:+(::Zero, b::$T) = b
    @eval Base.:+(a::$T, ::Zero) = a

    @eval Base.:*(::Zero, ::$T) = Zero()
    @eval Base.:*(::$T, ::Zero) = Zero()
end

Base.:+(a::One, b::One) = extern(a) + extern(b)
Base.:*(::One, ::One) = One()
for T in (:AbstractThunk, :Composite, :Any)
    if T != :Composite
        @eval Base.:+(a::One, b::$T) = extern(a) + b
        @eval Base.:+(a::$T, b::One) = a + extern(b)
    end

    @eval Base.:*(::One, b::$T) = b
    @eval Base.:*(a::$T, ::One) = a
end

Base.:+(a::AbstractThunk, b::AbstractThunk) = unthunk(a) + unthunk(b)
Base.:*(a::AbstractThunk, b::AbstractThunk) = unthunk(a) * unthunk(b)
for T in (:Composite, :Any)
    @eval Base.:+(a::AbstractThunk, b::$T) = unthunk(a) + b
    @eval Base.:+(a::$T, b::AbstractThunk) = a + unthunk(b)

    @eval Base.:*(a::AbstractThunk, b::$T) = unthunk(a) * b
    @eval Base.:*(a::$T, b::AbstractThunk) = a * unthunk(b)
end

function Base.:+(a::Composite{P}, b::Composite{P}) where P
    data = elementwise_add(backing(a), backing(b))
    return Composite{P, typeof(data)}(data)
end
function Base.:+(a::P, d::Composite{P}) where P
    net_backing = elementwise_add(backing(a), backing(d))
    if debug_mode()
        try
            return construct(P, net_backing)
        catch err
            throw(PrimalAdditionFailedException(a, d, err))
        end
    else
        return construct(P, net_backing)
    end
end
Base.:+(a::Composite{P}, b::P) where P = b + a

# We intentionally do not define, `Base.*(::Composite, ::Composite)` as that is not meaningful
# In general one doesn't have to represent multiplications of 2 differentials
# Only of a differential and a scaling factor (generally `Real`)
for T in (:Any,)
    @eval Base.:*(s::$T, comp::Composite) = map(x->s*x, comp)
    @eval Base.:*(comp::Composite, s::$T) = map(x->x*s, comp)
end


text:
using RecyclePermanent
using JLD2

include("../models/gaussianModel.jl")

@load "data/gaussianData.jld" ys θ
ys = ys[1:10]
model = makeGaussianModel(θ, ys)

N = 40 * model.n

println("Single-threaded")
rpio = RPIO{model.Tx}(N, model.n, 1)
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)

println("Multi-threaded")
rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)


text:
using FlexLinearAlgebra

export random_average_height, true_average_height

"""
`random_average_height(P,reps=100)` gives the average height of the elements of `P`
in a random linear extension. See `random_linear_extension`.
"""
function random_average_height(P::SimplePoset{T}, reps::Int=100)::Dict{T,Float64} where T
    total = FlexVector{Float64}(elements(P))
    n = card(P)
    for j=1:reps
        L = random_linear_extension(P)
        for i=1:n
            x = L[i]
            total[x] += (i-1)
        end
    end
    d = Dict{T,Float64}()
    for x in elements(P)
        d[x] = total[x]/reps
    end
    return d
end

"""
`true_average_height(P)` gives the average height of the elements of `P`
averaged over all linear extensions of `P`. See `all_linear_extensions`.
"""
function true_average_height(P::SimplePoset{T})::Dict{T,Float64} where T
    LX_list = collect(all_linear_extensions(P))  # VERY EXPENSIVE!
    nx = length(LX_list)
    n  = card(P)

    total = FlexVector{Float64}(elements(P))

    for L in LX_list
        for i=1:n
            x = L[i]
            total[x] += (i-1)
        end
    end

    d = Dict{T,Float64}()
    for x in elements(P)
        d[x] = total[x]/nx
    end
    return d
end


text:
# push!(LOAD_PATH,"../src/")
using Documenter, ClusterAnalysis

DocMeta.setdocmeta!(ClusterAnalysis, :DocTestSetup, :(using ClusterAnalysis); recursive=true)

makedocs(
    modules=[ClusterAnalysis],
    authors="AugustoCL <augustoleal72@gmail.com>, Elias Carvalho <eliascarvdev@gmail.com>",
    repo="https://github.com/AugustoCL/ClusterAnalysis.jl",
    sitename="ClusterAnalysis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://AugustoCL.github.io/ClusterAnalysis.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Algorithms" => [
            "K-Means" => "algorithms/kmeans.md",
            "DBSCAN" => "algorithms/dbscan.md"
        ],
        "API" => [
            "Functions" => "API/functions.md",
            "Types" => "API/types.md"
        ]
    ]
)

deploydocs(
    repo="github.com/AugustoCL/ClusterAnalysis.jl.git",
    devbranch = "main"
)


text:
mutable struct DESPOTSolver{S,A,O,B,RS} <: POMDPs.Solver
    belief::DESPOTBelief{S,A,O}
    bounds::B
    random_streams::RS
    root::VNode{S,A,O,B}
    node_count::Int64
    config::DESPOTConfig
    rng::AbstractRNG
    #preallocated for simulations
    curr_reward::DESPOTReward
    next_state::S
    curr_obs::O
    print_bounds::Bool

    # default constructor
    function DESPOTSolver{S,A,O,B,RS}(  ;
                            bounds::B = B(), #TODO: fix
                            rng::AbstractRNG = Base.GLOBAL_RNG,
                            search_depth::Int64 = 90,
                            main_seed::UInt32 = convert(UInt32, 42),
                            time_per_move::Float64 = 1.0,                 # sec
                            n_particles::Int64 = 500,
                            pruning_constant::Float64 = 0.0,
                            eta::Float64 = 0.95,
                            sim_len::Int64 = -1,
                            approximate_ubound::Bool = false,
                            tiny::Float64 = 1e-6,
                            max_trials::Int64 = -1,
                            rand_max::Int64 = 2147483647,
                            debug::Int64 = 0,
                            random_streams::RS = RandomStreams(n_particles, search_depth, main_seed),
                            next_state = S(),
                            curr_obs = O(),
                            print_bounds = false
                           ) where {S,A,O,B,RS}

        this = new()

        # supplied variables
        this.bounds = bounds

        # Instantiate and initialize config
        this.config = DESPOTConfig()
        this.config.search_depth = search_depth
        this.config.main_seed = main_seed
        this.config.time_per_move = time_per_move
        this.config.n_particles = n_particles
        this.config.pruning_constant = pruning_constant
        this.config.eta = eta
        this.config.sim_len = sim_len
        this.config.approximate_ubound = approximate_ubound
        this.config.tiny = tiny
        this.config.max_trials = max_trials
        this.config.rand_max = rand_max
        this.config.debug = debug
        this.rng = rng
        this.next_state = next_state
        this.curr_obs = curr_obs
        this.curr_reward = 0.0
        this.random_streams = random_streams
        this.print_bounds = print_bounds
        return this
    end
end

function init_solver{S,A,O,B}(solver::DESPOTSolver{S,A,O,B}, pomdp::POMDPs.POMDP{S,A,O})

    # Instantiate random streams

    fill_random_streams!(solver.random_streams, solver.config.rand_max)
    init_bounds(solver.bounds, pomdp, solver.config)

    return nothing
end

function new_root{S,A,O,B}(solver::DESPOTSolver{S,A,O,B},
                  pomdp::POMDP{S,A,O},
                  belief::DESPOTBelief{S})

    solver.belief = belief
    solver.root = VNode{S,A,O,B}(
                        pomdp,
                        belief.particles,
                        solver.bounds,
                        0,
                        1.0,
                        false,
                        solver.config)
    return nothing
end

function search{S,A,O,B}(solver::DESPOTSolver{S,A,O,B}, pomdp::POMDP{S,A,O})
    n_trials::Int64 = 0
    start_time::Float64 = time()
    stop_now::Bool = false

    if solver.print_bounds
        @printf("Before: lBound = %.10f, uBound = %.10f\n", solver.root.lbound, solver.root.ubound)
    end

    while ((excess_uncertainty(solver.root.lbound,
                                solver.root.ubound,
                                solver.root.lbound,
                                solver.root.ubound,
                                0,
                                solver.config.eta,
                                discount(pomdp)) > 1e-6)
                                && !stop_now)

        trial(solver, pomdp, solver.root, n_trials)
        n_trials += 1
        if ((solver.config.max_trials > 0) && (n_trials >= solver.config.max_trials)) ||
        ((solver.config.time_per_move > 0) && ((time() - start_time) >= solver.config.time_per_move))
            stop_now = true
        end
    end

    if solver.print_bounds
        @printf("After:  lBound = %.10f, uBound = %.10f\n", solver.root.lbound, solver.root.ubound)
        @printf("Number of trials: %d\n", n_trials)
    end

    if solver.config.pruning_constant != 0.0
        total_pruned = prune(solver.root) # Number of non-child belief nodes pruned
        act = solver.root.pruned_action
        default = default_action(solver.bounds, pomdp, solver.root.particles, solver.config)
        return (act == -1 ? default : act), n_trials #TODO: fix actions
    elseif !solver.root.in_tree
        println("Root not in tree")
        default = default_action(solver.bounds, pomdp, solver.root.particles, solver.config)
        return default, n_trials
    else
        return get_lb_action(solver.root, solver.config, discount(pomdp)), n_trials
    end
    return nothing
end

function trial{S,A,O,B}(solver::DESPOTSolver{S,A,O,B}, pomdp::POMDP{S,A,O}, node::VNode{S,A,O,B}, n_trials::Int64)

    n_nodes_added::Int64 = 0
    ubound::Float64 = 0.0

    if (node.depth >= solver.config.search_depth) || isterminal(pomdp, node.particles[1].state)
      return 0 # nodes added
    end

    if isempty(node.q_nodes)
        expand_one_step(solver, pomdp, node, solver.bounds)
    end

    a_star::A = node.best_ub_action

    i_star::Int, weighted_eu_star::Float64 =
                               get_best_weuo(node.q_nodes[a_star],
                                             solver.root,
                                             solver.config,
                                             discount(pomdp)) # it's an array!

    o_star, vnode_star = node.q_nodes[a_star].obs_and_nodes[i_star]

    if weighted_eu_star > 0.0
        add(solver.belief.history, a_star, o_star)
        n_nodes_added = trial(solver,
                            pomdp,
                            vnode_star,
                            n_trials)
        remove_last(solver.belief.history)
    end
    node.n_tree_nodes += n_nodes_added

    # Backup
    potential_lbound = node.q_nodes[a_star].first_step_reward +
                    discount(pomdp) * get_lower_bound(node.q_nodes[a_star])
    node.lbound = max(node.lbound, potential_lbound)

    # As the upper bound of a_star may become smaller than the upper bound of
    # another action, we need to check all actions - unlike the lower bound.
    node.ubound = -Inf

    for a in iterator(actions(pomdp))
        ubound = node.q_nodes[a].first_step_reward +
            discount(pomdp) * get_upper_bound(node.q_nodes[a])
        if ubound > node.ubound
            node.ubound = ubound
            node.best_ub_action = a
        end
    end

    # Sanity check
    if (node.lbound > node.ubound + solver.config.tiny)
        println("depth = $(node.depth)")
        warn("Lower bound ($(node.lbound)) is higher than upper bound ($(node.ubound))")
    end

    if !node.in_tree
        node.in_tree = true
        node.n_tree_nodes += 1
        n_nodes_added += 1
    end
    return n_nodes_added
end

function expand_one_step{S,A,O,B}(solver::DESPOTSolver{S,A,O}, pomdp::POMDP{S,A,O}, node::VNode{S,A,O,B}, bounds::B)

    q_star::Float64 = -Inf
    first_step_reward::Float64 = 0.0
    remaining_reward::Float64 = 0.0
    rng = create_rng(solver.random_streams)

    for curr_action in iterator(actions(pomdp))
        obs_to_particles = Dict{O,Vector{DESPOTParticle{S}}}()

        for p in node.particles

            set_rng_state!(rng, solver.random_streams, p.id+1, node.depth+1)

            step(
                solver,
                pomdp,
                p.state,
                rng,
                curr_action)

            if !haskey(obs_to_particles, solver.curr_obs)
                obs_to_particles[solver.curr_obs] = DESPOTParticle{S}[]
            end

            push!(obs_to_particles[solver.curr_obs],
                  DESPOTParticle(solver.next_state,
                  p.id,
                  p.weight))
            first_step_reward += solver.curr_reward * p.weight
        end

        first_step_reward /= node.weight

        new_qnode = QNode{S,A,O,B}(
                          pomdp,
                          bounds,
                          obs_to_particles,
                          node.depth,
                          curr_action,
                          first_step_reward,
                          solver.belief.history,
                          solver.config)
        node.q_nodes[curr_action] = new_qnode
        remaining_reward = get_upper_bound(new_qnode)

        if (first_step_reward + discount(pomdp)*remaining_reward) > (q_star + solver.config.tiny)
            q_star = first_step_reward + discount(pomdp) * remaining_reward
            node.best_ub_action = curr_action
        end

        first_step_reward = 0.0
    end # for a
    return node
end

# fill in pre-allocated variables
function step{S,A,O,B}(solver::DESPOTSolver{S,A,O,B},
              pomdp::POMDPs.POMDP{S,A,O},
              state::S,
              rng::AbstractRNG,
              action::A)

    solver.next_state, solver.curr_obs, solver.curr_reward =
        POMDPs.generate_sor(pomdp, state, action, rng)

    return nothing
end


text:
"""
    difference(ε::Real, y::T, x::T) where {T}

Computes `dx` where `dx` is defined s.t.
```julia
y = x + ε * dx
```
where `dx` is a valid tangent type for `x`.

If `(y - x) / ε` is defined, then this operation is equivalent to doing that. For functions
where these operations aren't defined, `difference` can still be defined without commiting
type piracy while `-` and `/` cannot.
"""
difference(::Real, ::T, ::T) where {T<:Symbol} = DoesNotExist()
difference(::Real, ::T, ::T) where {T<:AbstractChar} = DoesNotExist()
difference(::Real, ::T, ::T) where {T<:AbstractString} = DoesNotExist()
difference(::Real, ::T, ::T) where {T<:Integer} = DoesNotExist()

difference(ε::Real, y::T, x::T) where {T<:Number} = (y - x) / ε

difference(ε::Real, y::T, x::T) where {T<:StridedArray} = difference.(ε, y, x)

function difference(ε::Real, y::T, x::T) where {T<:Tuple}
    return Composite{T}(difference.(ε, y, x)...)
end

function difference(ε::Real, ys::T, xs::T) where {T<:NamedTuple}
    return Composite{T}(; map((y, x) -> difference(ε, y, x), ys, xs)...)
end

function difference(ε::Real, y::T, x::T) where {T}
    if !isstructtype(T)
        throw(ArgumentError("Non-struct types are not supported by this fallback."))
    end

    field_names = fieldnames(T)
    if length(field_names) > 0
        tangents = map(field_names) do field_name
            difference(ε, getfield(y, field_name), getfield(x, field_name))
        end
        return Composite{T}(; NamedTuple{field_names}(tangents)...)
    else
        return NO_FIELDS
    end
end


text:
"""
Mot fileformat

Militairy opening times for usage with route choices

See EUROCONTROL NEST Manual for Mot fileformat description

Output:

Mot fileformat example
1 MIL1 1000 1200
1 MIL1 1400 1700
1 MIL2 0900 1300
1 MIL3 1100 1300
1 MIL4 0900 1800

2 MIL1 1000 1200 055 245
2 MIL1 1400 1700 085 285
2 MIL2 0900 1300 000 600
2 MIL3 1100 1300 125 245
2 MIL4 0900 1800 000 285

3 ED 2200 0500 EX
3 LF 2300 0600 EX
3 LI 2200 0500 E
"""
module Mot

using ..util
using Dates

export read

const hhmm = DateFormat("HHMM")

struct OnlyTime
    begintime::Union{Time, Missing}
    endtime::Union{Time, Missing}
end

struct TimeLevel
    begintime::Union{Time, Missing}
    endtime::Union{Time, Missing}
    lowerFL::Int64
    upperFL::Int64
end

struct TimeType
    begintime::Union{Time, Missing}
    endtime::Union{Time, Missing}
    type::String
end

struct OpeningTimes
    onlytimes::Vector{OnlyTime}
    timelevels::Vector{TimeLevel}
    timetypes::Vector{TimeType}
end

function read(file)
    militairyzones = Dict{String, OpeningTimes}()
    onlytimes = Vector{OnlyTime}()
    timelevels = Vector{TimeLevel}()
    timetypes = Vector{TimeType}()
    for line in eachline(file)
        line_elements = split(line)
        if length(line_elements) > 3
            begintime = format_time.(line_elements[3], hhmm)
            endtime = format_time.(line_elements[4], hhmm)
            if haskey(militairyzones, line_elements[2])
                openingtimes  = militairyzones[line_elements[2]]
                onlytimes = openingtimes.onlytimes
                timelevels = openingtimes.timelevels
                timetypes = openingtimes.timetypes
            else
                onlytimes = Vector{OnlyTime}()
                timelevels = Vector{TimeLevel}()
                timetypes = Vector{TimeType}()
            end

            if line_elements[1] == "1"
                onlytimes = vcat(onlytimes, OnlyTime(begintime, endtime))
            elseif line_elements[1] == "2"
                lowerFL = parse(Int64, line_elements[5])
                upperFL = parse(Int64, line_elements[6])
                timelevels = vcat(timelevels, TimeLevel(begintime, endtime, lowerFL, upperFL))
            elseif line_elements[1] == "3"
                timetypes = vcat(timetypes, TimeType(begintime, endtime, line_elements[5]))
            else
            end
            militairyzones[line_elements[2]] = OpeningTimes(onlytimes, timelevels, timetypes)
        end
    end
    return militairyzones
end

end # module


text:
#=

  https://www.allendowney.com/blog/2018/10/18/how-tall-is-a/
  """
  Here are a series of problems I posed in my Bayesian statistics class:

  1) Suppose you meet an adult resident of the U.S. who is 170 cm tall.
     What is the probability that they are male?

  2) Suppose I choose two U.S. residents at random and A is taller than B.
     How tall is A?

  3) In a room of 10 randomly chosen U.S. residents, A is the second tallest.
     How tall is A?
     And what is the probability that A is male?

  As background: For adult male residents of the US, the mean and standard deviation of
  height are 178 cm and 7.7 cm. For adult female residents the corresponding stats
  are 163 cm and 7.3 cm.  And 51% of the adult population is female.

  """

  Cf ~/blog/how_tall_is_a.blog
     ~/webppl/how_tall_is_a.wppl

=#

using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
Turing.setadbackend(:forwarddiff)
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:tracker)
# Turing.setadbackend(:zygote)

include("jl_utils.jl")

#=
  Model 1:

  1) Suppose you meet an adult resident of the U.S. who is 170 cm tall.
  What is the probability that they are male?

  According to https://www.allendowney.com/blog/2018/10/18/how-tall-is-a/
    female 0.5432225483131837
    male 0.45677745168681627

  This model:
  Distributions of variable gender (num:0)
    2.00000 =>    5406  (0.540600)
    1.00000 =>    4594  (0.459400)
=#
@model function how_tall_is_1(height=170)
    male = 1
    female = 2
    gender ~ Categorical(simplex([0.49,0.51]))
    height ~ gender == male ? Normal(178,7.7) : Normal(163,7.3)

    # true ~ Dirac(height == 170.0)

end

#=
println("Model 1")
model = how_tall_is_1(170)
# chns = sample(model, Prior(), 10_000)
chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 1_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)

#chns = sample(model, NUTS(), 1000)
# chns = sample(model, HMC(0.1, 5), 1000)
# chns = sample(model, HMCDA(0.15, 0.65), 1000)
# chns = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)

display(chns)
println("male: 1 female: 2")
show_var_dist_pct(chns,:gender)
# show_var_dist_pct(chns,:height,20)
=#

#=

  2) Suppose I choose two U.S. residents at random and A is taller than B.
  How tall is A?

  Solution from https://www.allendowney.com/blog/2018/10/18/how-tall-is-a/
  A: 176.67506663725212
  B: 164.05298129722925

  This model:
    ...
    height[1]   176.4384    9.0135     0.2850    0.1777   824.6929    1.0007      179.5152
    height[2]   164.5469    8.4415     0.2669    0.2017   828.7247    0.9998      180.3928

=#

@model function how_tall_is_2()
    male = 1
    female = 2
    gender = tzeros(2)
    height = tzeros(2)
    for p in 1:2
        gender[p] ~ Categorical(simplex([0.49,0.51]))
        height[p] ~ gender[p] == male ? Normal(178,7.7) : Normal(163,7.3)
    end
    true ~ Dirac(height[1] > height[2])

end

#=
println("\n\nModel 2")
model = how_tall_is_2()
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(15), 1_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
println("male: 1 female: 2")
show_var_dist_pct(chns,Symbol("gender[1]"))
show_var_dist_pct(chns,Symbol("gender[2]"))
=#

#=

  3) In a room of 10 randomly chosen U.S. residents, A is the second tallest.
  How tall is A?
  And what is the probability that A is male?

    Solution from https://www.allendowney.com/blog/2018/10/18/how-tall-is-a/
    A: 181.60660153115973

=#

# Note:
# This model does NOT work as expected. It's fine for n=2 (see how_tall_is_2), but not the
# generalized version of n=10. See how_tall_is_3_2() for a working model (thanks to Seth Axen).
# See the discussion here:
#  https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765
#
@model function how_tall_is_3()
    male = 1
    female = 2
    n = 10
    gender = tzeros(n)
    height = tzeros(n)
    for p in 1:n
        gender[p] ~ Categorical(simplex([0.49,0.51]))
        height[p] ~ gender[p] == male ? Normal(178,7.7) : Normal(163,7.3)
    end

    for p in 1:n-1
        # true ~ Dirac(height[p] > height[p+1])
        height[p] > height[p+1] || begin Turing.@addlogprob! -Inf; end
        # height[p] > height[p+1] || begin Turing.@addlogprob! -2; end
    end

    a_is_male ~ Dirac(gender[2] == male)
end

#=
println("\n\nModel 3")
model = how_tall_is_3()
# chns = sample(model, Prior(), 10_000)
chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)
display(chns)
println("male: 1 female: 2")
show_var_dist_pct(chns,:a_is_male)
=#

#
# Thanks to Seth Axen for this approach.
# From https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765/2
# Note that the heights[i] are all about 170cm; it's the height_ordered that we are
# interested in:
#=
   height_ordered[1]   186.706
   height_ordered[2]   181.5520
   height_ordered[3]   177.8786
   height_ordered[4]   174.6669
   height_ordered[5]   171.6923
   height_ordered[6]   168.7933
   height_ordered[7]   165.8899
   height_ordered[8]   162.8638
   height_ordered[9]   159.4338
  height_ordered[10]   154.5771
            height_A   181.5520
          isfemale_A     0.0893
=#
@model function how_tall_is_3_2()
    male = 1
    female = 2
    n = 10
    # This is from Seth Axen:
    isfemale ~ filldist(Bernoulli(0.51), n)
    height = tzeros(n) # Vector{T}(undef, n)
    for i in 1:n
        if isfemale[i]
            height[i] ~ Normal(163, 7.3)
        else
            height[i] ~ Normal(178, 7.7)
        end
    end
    # hakank: I added this since I actually wanted the
    #         weights of all people.
    height_ordered = tzeros(n) # Vector{T}(undef, n)
    perm = sortperm(height)
    for i in 1:n
        height_ordered[i] ~ Dirac(height[perm[n-i+1]])
    end

    # From Seth:
    id_A = perm[n-1]
    height_A ~ Dirac(height[id_A])
    isfemale_A ~ Dirac(isfemale[id_A])

end

model = how_tall_is_3_2()
# Comment from Seth:
# "no inference necessary! We just sample directly from the prior-predictive distribution"
chns = sample(model, Prior(), 10_000)
display(chns)
show_var_dist_pct(chns,:isfemale_A)
# Testing with IS(), same result
# chns = sample(model, IS(), 10_000)
# display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function assign_scenario(s::BasePriceScenario, model::JuMP.Model, config::OptionsV2016R2, params::ParametersV2016, vars::VariablesV2016)
    JuMP.set_value(params.ψ₂, 0.000001);
    for i in 1:config.N
        JuMP.set_upper_bound(vars.CPRICE[i], 10000.0);
    end
    optimize!(model);

    photel = value.(vars.CPRICE);

    for i in 1:config.N
        if i <= config.tnopol
            JuMP.set_upper_bound(vars.CPRICE[i], max(photel[i],params.cpricebase[i]));
        end
    end
    JuMP.set_value(params.ψ₂, config.ψ₂);
end

function assign_scenario(s::BasePriceScenario, model::Model, config::OptionsV2016R, params::ParametersV2016, vars::VariablesV2016)
    JuMP.set_value(params.ψ₂, 0.0);
    optimize!(model);

    photel = value.(vars.CPRICE);

    for i in 1:config.N
        if i <= config.tnopol
            JuMP.set_upper_bound(vars.CPRICE[i], max(photel[i],params.cpricebase[i]));
        end
    end
    JuMP.set_value(params.ψ₂, config.ψ₂);
end

function assign_scenario(s::OptimalPriceScenario, model::Model, config::Options, params::ParametersV2016, vars::VariablesV2016)
    JuMP.fix(vars.μ[1], config.μ₀; force=true);
end

function assign_scenario(s::Scenario, model::JuMP.Model, config::Options, params::ParametersV2016, vars::VariablesV2016)
    if typeof(config) <: OptionsV2016R
        error("$(s) is not a valid scenario for v2016R beta");
    else
        error("$(s) is not a valid scenario for v2016R2");
    end
end


text:
@bf read = [right!, read!, right!, inc!]
@bf write = [left!, write!, right!, drop]

@bf print = [[write], eachr]
@bf println = [print, '\n', write]

@bf readln = [nil, [head, '\n', !=], [read, push], doloop, pop, drop]

@bf panic = ["PANIC: ", print, println, halt]

@bf prompt = [print, readln]

@bf whileempty_ = [[isempty], dip, swap,
                   [[drop], dip, [call], keep, whileempty_], [drop], iff]
@bf whileempty = [nil, swap, whileempty_]

# Printing
@bf charstr = ['\'', swap, '\'', 3]

# Parsing
@bf isdigit = [['0', >=], ['9', <=], bi, and]
@bf digiterr = [charstr, " is not a digit", cat]
@bf assertdigit = [dup, isdigit, [digiterr, panic], unless]
@bf toint = [assertdigit, '0', -]
@bf parseint = [reverse, 0, [[toint], [10, *], bi_, +], fold]

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
description(i::AnchoredInterval{P}) where P = description(i, P > zero(P) ? "B" : "E")

function description(interval::AnchoredInterval{P, T}, s::String) where {P, T}
    return string(
        first(inclusivity(interval)) ? '[' : '(',
        description(anchor(interval), abs(P), s),
        last(inclusivity(interval)) ? ']' : ')',
    )
end

function description(interval::AnchoredInterval{P, ZonedDateTime}, s::String) where P
    return string(
        first(inclusivity(interval)) ? '[' : '(',
        description(anchor(interval), abs(P), s),
        anchor(interval).zone.offset,
        last(inclusivity(interval)) ? ']' : ')',
    )
end

function description(dt::Date, p::Period, suffix::String)
    ds = @sprintf("%04d-%02d-%02d", year(dt), month(dt), day(dt))
    return "$(prefix(p))$suffix $ds"
end

function description(dt::AbstractDateTime, p::Period, suffix::String)
    ts = time_string(dt, p)

    # If we would display only HE00, display HE24 for the previous day instead
    if ts == "00" && suffix == "E"
        P, max_val = coarserperiod(typeof(p))
        dt -= oneunit(P)
        ts = string(max_val)
    end

    @static if VERSION < v"0.7.0-DEV.3216"
        isa(dt, Time) && return string(prefix(p), suffix, " ", ts)
    end

    ds = @sprintf("%04d-%02d-%02d", year(dt), month(dt), day(dt))

    if p isa TimePeriod
        return "$ds $(prefix(p))$suffix$ts"
    else
        ts = (ts == "00:00:00" && !isa(dt, ZonedDateTime)) ? "" : " $ts"
        return "$(prefix(p))$suffix $ds$ts"
    end
end

function time_string(dt::AbstractDateTime, p::Period)
    t = (hour(dt), minute(dt), second(dt), millisecond(dt))
    if p isa Hour && all(t[2:end] .== 0)
        return @sprintf("%02d", t[1])
    elseif p isa Minute && all(t[3:end] .== 0)
        return @sprintf("%02d:%02d", t[1:2]...)
    elseif !isa(p, Millisecond) && t[4] == 0
        return @sprintf("%02d:%02d:%02d", t[1:3]...)
    else
        return @sprintf("%02d:%02d:%02d.%03d", t...)
    end
end

function prefix(p::Period)
    string(value(p) == 1 ? "" : value(p), prefix(typeof(p)))
end

prefix(::Type{Year}) = "Y"
prefix(::Type{Month}) = "Mo"
prefix(::Type{Day}) = "D"
prefix(::Type{Hour}) = "H"
prefix(::Type{Minute}) = "M"
prefix(::Type{Second}) = "S"
prefix(::Type{Millisecond}) = "ms"


text:
export StateOverriddenEnv

"""
    StateOverriddenEnv(f, env)

Apply `f` to override `state(env)`.

!!! note
    If the meaning of state space is changed after apply `f`, one should
    manually redefine the `RLBase.state_space(env::YourSpecificEnv)`.
"""
struct StateOverriddenEnv{F,E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    f::F
end

StateOverriddenEnv(f) = env -> StateOverriddenEnv(f, env)

RLBase.state(env::StateOverriddenEnv, args...; kwargs...) =
    env.f(state(env.env, args...; kwargs...))


text:
# This file is a part of UnsafeArrays.jl, licensed under the MIT License (MIT).

"""
    uview(A::AbstractArray, I...)

Unsafe equivalent of `view`. May return an `UnsafeArray`, a standard
`SubArray` or `A` itself, depending on `I...` and the type of `A`.

As `uview` may return an `UnsafeArray`, `A` itself and it's contents *must* be
protected from garbage collection (e.g. via `GC.@preserve` on Julia > v0.6)
and memory reallocation while the view is in use.

Use `uviews(f::Function, As::AbstractArray...)` or `@uviews A ... expr` to
use unsafe views of one or multiple arrays with automatically GC protection.

```
uviews(A, B, ...) do (A_u, B_u, ...)
    # Do something with the unsafe views A_u, B_u, ...
    # Code here must not resize/append/etc. A, B, ...
end
```

To provide support for `uview` for custom array types, add methods to
function `UnsafeArrays.unsafe_uview`.
"""
function uview end
export uview

Base.@propagate_inbounds uview(A::AbstractArray) = unsafe_uview(A)

Base.@propagate_inbounds function uview(A::AbstractArray, idx, I...)
    J = Base.to_indices(A, (idx, I...))
    @boundscheck checkbounds(A, J...)
    unsafe_uview(A, J...)
end

Base.@propagate_inbounds uview(A::UnsafeArray) = A

"""
    UnsafeArray.unsafe_uview(A::AbstractArray, I::Vararg{Base.ViewIndex,N})
    UnsafeArray.unsafe_uview(A::AbstractArray, i::Base.ViewIndex)
    UnsafeArray.unsafe_uview(A::AbstractArray)

To support `uview` for custom array types, add methods to `unsafe_uview`
instead of `uview`. Implementing
`UnsafeArray.unsafe_uview(A::CustomArrayType)` will often be sufficient.

It may be necessary to provide a custom implementation of `Base.deepcopy`
for types with a custom implementation of `unsafe_uview`, to ensure that
the result of `deepcopy` does *NOT* contain unsafe views.
"""
function unsafe_uview end

Base.@propagate_inbounds unsafe_uview(A::AbstractArray{T,N}, I::Vararg{Base.ViewIndex,N}) where {T,N} =
    Base.unsafe_view(unsafe_uview(A), I...)

Base.@propagate_inbounds unsafe_uview(A::AbstractArray{T,N}, i) where {T,N} =
    Base.unsafe_view(unsafe_uview(A), i::Base.ViewIndex)

Base.@propagate_inbounds unsafe_uview(A::AbstractArray{T,N}) where {T,N} = A

Base.@propagate_inbounds unsafe_uview(A::UnsafeArray{T,N}) where {T,N} = A

Base.@propagate_inbounds unsafe_uview(A::DenseArray{T,N}) where {T,N} =
    _maybe_unsafe_uview(Val{isbitstype(T)}(), A)

Base.@propagate_inbounds unsafe_uview(A::SubArray{T,N}) where {T,N} =
    _maybe_unsafe_uview(Val{isbitstype(T) && Base.iscontiguous(typeof(A))}(), A)

Base.@propagate_inbounds function _maybe_unsafe_uview(unsafe_compatible::Val{true}, A::AbstractArray{T,N}) where {T,N}
    @boundscheck _require_one_based_indexing(A)
    UnsafeArray{T,N}(unsafe_compatible, pointer(A), size(A))
end

Base.@propagate_inbounds _maybe_unsafe_uview(unsafe_compatible::Val{false}, A::AbstractArray{T,N}) where {T,N} = A

"""
    uviews(f::Function, As::AbstractArray...)

Equivalent to `f(map(uview, As)...)`. Automatically protects the array(s)
`As` from garbage collection during execution of `f`.

Example:

```
uviews(A, B, ...) do (A_u, B_u, ...)
    # Do something with the unsafe views A_u, B_u, ...
    # Code here must not resize/append/etc. A, B, ...
end
```

In many cases, it may be preferable to use `@uviews` insted of `views`.
"""
function uviews end
export uviews

@inline function uviews(f::Function, As::AbstractArray...)
    GC.@preserve(As, f(map(uview, As)...))
end

"""
    @uviews A B ... expr

Replace arrays `A`, `B`, ... by uview(`A`), uview(`B`), ... during execution
of `expr`, while protecting the original arrays from garbage collection.

Equivalent to

```
GC.@preserve A B ... begin
    let A = uview(A), B = uview(B), ...
        expr
    end
end
```

The unsafe views must not be allowed to escape the scope of `expr`. The
original arrays must not be resized/appended/etc. during the execution of
`expr`.
"""
macro uviews(args...)
    syms = args[1:end-1]
    expr = args[end]

    binds = Expr(:block)
    for s in syms
        s isa Symbol || error("@uviews targets must be a symbols")
        push!(binds.args, :($s = UnsafeArrays.uview($s)))
    end
    let_expr = Expr(:let, binds, expr)
    esc(:(GC.@preserve $(syms...) $(let_expr)))
end

export @uviews

function Base.mightalias(A::UnsafeArray, B::UnsafeArray)
    pfA = pointer(A, firstindex(A))
    plA = pointer(A, lastindex(A))
    pfB = pointer(B, firstindex(B))
    plB = pointer(B, lastindex(B))

    (pfA <= pfB <= plA) || (pfA <= plB <= plA) || (pfB <= pfA <= plB)
end

function Base.mightalias(A::UnsafeArray, B::AbstractArray)
    @uviews B begin
        if typeof(B) <: UnsafeArray
            Base.mightalias(A, B)
        else
            false
        end
    end
end

Base.mightalias(A::AbstractArray, B::UnsafeArray) = Base.mightalias(B, A)

Base.mightalias(A::SubArray{T,N,<:UnsafeArray}, B::AbstractArray) where {T,N} =
    Base.mightalias(parent(A), B)

Base.mightalias(A::SubArray{T1,N1,<:UnsafeArray}, B::SubArray{T2,N2,<:UnsafeArray}) where {T1,N1,T2,N2} =
    Base.mightalias(parent(A), parent(B))

Base.mightalias(A::SubArray{T,N,<:UnsafeArray}, B::UnsafeArray) where {T,N} =
    Base.mightalias(parent(A), B)

Base.mightalias(A::AbstractArray, B::SubArray{T,N,<:UnsafeArray}) where {T,N} =
    Base.mightalias(B, A)

Base.mightalias(A::UnsafeArray, B::SubArray{T,N,<:UnsafeArray}) where {T,N} =
    Base.mightalias(B, A)


text:
module ExactDiagonalizationTools

using QuantumModels

using LinearAlgebra
using SparseArrays
#using DataStructures
using SymTensors
using Random
using QuadGK
RLorCX = Union{Float64, ComplexF64}
⊗ = kron

import LinearAlgebra: normalize!

import Base: +, *
include("model.jl")

include("ket.jl")
include("xxz.jl")
include("longxxz.jl")
include("qitf.jl")
include("hopping.jl")
include("correlationmatrix.jl")

export KetState
export opextend
export normalize!
export randket
export apply1site
export measure1point

export xxz_hamiltonian
export xxz_longrange
export j1j2_explicit
export qitf_hamiltonian, qitf_bondtensor
export qitf_energy_exact

export hopping_chain
export nnhoppingchain
export correlationmatrix

# Model stuff to be separated
export Point2D, Bond2D, UnitCell2D
export triangular_unitcell
export enlargeunitcell
export makehamiltonian
export makemodel
export makemodelJW

export makemodelJW_u1sym

export generatebdg

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
module history

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
@doc "`Faces(t)` returns an iterator for faces from representation of topology `t`" Faces
@doc "`Edges(t)` returns an iterator for edges from representation of topology `t`" Edges

Faces(t::PlainDS) = t
Edges(t::PlainDS) = filter(x->x[1]<x[2],decompose(Face{2,Int},t)) ### Hopefully works

start(iter::FaceRing{PlainDS}) = find_triangle_vertex(iter.v,iter.t)
done(iter::FaceRing{PlainDS},ti::Int) = ti<=length(iter.t) ? false : true

function next(iter::FaceRing{PlainDS},i::Int)
    v = iter.v
    nexti = find_triangle_vertex(iter.v,iter.t[i+1:end]) + i  # possible botleneck here
    return i, nexti
end

Base.iterate(iter::FaceRing{PlainDS}) = next(iter,start(iter))
function Base.iterate(iter::FaceRing{PlainDS},ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end

start(iter::EdgeRing{PlainDS}) = find_triangle_vertex(iter.v,iter.t)
done(iter::EdgeRing{PlainDS},ti::Int) = ti<=length(iter.t) ? false : true

function next(iter::EdgeRing{PlainDS},i::Int)

    v = iter.v
    nexti = find_triangle_vertex(iter.v,iter.t[i+1:end]) + i  # possible botleneck here
    face = iter.t[i]
    w = face.==v
    cw = w[[3,1,2]]
    ccw = w[[2,3,1]]

    return (face[cw]...,face[ccw]...), nexti
end

Base.iterate(iter::EdgeRing{PlainDS}) = next(iter,start(iter))
function Base.iterate(iter::EdgeRing{PlainDS},ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end

start(iter::VertexRing{PlainDS}) = find_triangle_vertex(iter.v,iter.t)
done(iter::VertexRing{PlainDS},ti::Int) = ti<=length(iter.t) ? false : true

function next(iter::VertexRing{PlainDS},ti::Int)

    v = iter.v
    faces = iter.t
    face = faces[ti]
    cw = (face.==v)[[3,1,2]]
    vi, = face[cw]

    nexti = find_triangle_vertex(iter.v,iter.t[ti+1:end]) + ti  # possible botleneck
    return vi, nexti
end

Base.iterate(iter::VertexRing{PlainDS}) = next(iter,start(iter))
function Base.iterate(iter::VertexRing{PlainDS},ti::Int)
    if done(iter,ti)
        return nothing
    else
        return next(iter,ti)
    end
end


text:
__precompile__(true)

module SpaceExParser

#=========================================
Load dependencies
=========================================#

using DataStructures: OrderedDict
using EzXML, HybridSystems, LazySets, MathematicalSystems, Reexport, SymEngine

#=========================
Input/Output functions
==========================#
export readsxmodel
include("io.jl")

#=========================
Parsing SpaceExParser files
==========================#
include("parse.jl")

#================================================
Converting symbolic expressions and systems types
=================================================#
include("convert.jl")
include("symbolic.jl")

#==========================
Writing systems to new file
===========================#
export writesxmodel
include("write.jl")

end # module


text:
module StackCollections

struct Unsafe end
const unsafe = Unsafe()

abstract type AbstractStackSet <: AbstractSet{Int} end

@static if VERSION >= v"1.4"
    Base.hasfastin(::Type{AbstractStackSet}) = true
end

@noinline throw_empty_err() = throw(ArgumentError("collection must be non-empty"))

include("stackvector.jl")
include("digitset.jl")
include("stackset.jl")
include("onehotvector.jl")

export StackVector, DigitSet, StackSet, setindex, push, pop, delete,
       complement, isdisjoint, OneHotVector, reverse

end


text:
function read_fasta(filepath::String)
    @info "Reading FASTA file $filepath"
    open(FASTA.Reader, filepath) do f
        nread = 0
        record = FASTA.Record()
        seqs = Dict{String, LongSequence{DNAAlphabet{4}}}()
        while !eof(f)
            read!(f, record)
            nread += 1
            recname = FASTA.identifier(record)
            seqs[recname] = FASTA.sequence(LongSequence{DNAAlphabet{4}}, record)
        end
        if nread > length(seqs)
            error("Only retained $(length(seqs)) of $nread records read. This is commonly caused by sequences with duplicate names. Check your FASTA file.")
        end
        @info "Done"
        return seqs
    end
end


text:
#=
Basic thermodynamic properties
=#
using LinearAlgebra

include("solvecubic.jl")

@inline ncomponents(mix::BrusilovskyEoSMixture) = length(mix.components)
@inline Base.length(mix::BrusilovskyEoSMixture) = length(mix.components)

"""
    a_coef(substance, RT)

Return EoS coefficient ``a(T, Ψ) = a_c ϕ(T, Ψ)`` of `substance` at given `RT`.
```math
ϕ(T, Ψ) = [ 1 + Ψ (1 - T_r^0.5) ]^2
```
Reference: Brusilovsky2002[section: 5.5.2 (algorithm step 3), eq: (4.34), see p.142 and p.164]
"""
function a_coef(substance::BrusilovskyEoSComponent, RT::Real)
    psi = substance.Psi
    phi = ( 1.0 + psi * (1.0 - sqrt(RT / substance.RTc)) )^2
    return substance.ac * phi
end

"""
    pressure(substance, υ, RT)

Compute pressure (Pa) of `substance` at given molar volume `υ` (m³ mol⁻¹) and
thermal energy `RT` (J mol⁻¹).
"""
function pressure(substance::BrusilovskyEoSComponent, υ::Real, RT::Real)
    acoeff = a_coef(substance, RT)
    b = substance.b
    c = substance.c
    d = substance.d
    return RT / (υ - b) - acoeff / ((υ + c) * (υ + d))
end

"""
    pressure(substance, nmol, V, RT)

Compute pressure (Pa) of `substance` at given number of moles `nmol` (mol),
total volume `V` (m³) and thermal energy `RT` (J mol⁻¹).
"""
function pressure(substance::BrusilovskyEoSComponent, nmol::Real, V::Real, RT::Real)
    acoeff = a_coef(substance, RT)
    b = substance.b
    c = substance.c
    d = substance.d
    υ = V / nmol
    return RT / (υ - b) - acoeff / ((υ + c) * (υ + d))
end

"""
    pressure(substance; nmol = 1, volume, temperature)

Compute pressure (Pa) of `substance` at given number of moles `nmol` (mol),
total volume (m³) and temperature (K).
"""
function pressure(
    substance::BrusilovskyEoSComponent;
    nmol::Real = 1,
    volume::Real,
    temperature::Real,
)
    RT = GAS_CONSTANT_SI * temperature
    return pressure(substance, nmol, volume, RT)
end

"""
    wilson_saturation_pressure(substance, RT)

Return approximate saturation pressure of `substance` at `RT` (J mol⁻¹).

Reference: Brusilovsky2002 [p 272, eq 5.4]
"""
function wilson_saturation_pressure(substance::BrusilovskyEoSComponent, RT::Real)
    return wilson_saturation_pressure(
        substance.Pc, substance.RTc, substance.acentric_factor, RT
    )
end

"""
    wilson_saturation_pressure(Pc::Real, RTc::Real, acentric_factor::Real, RT::Real)

Approximate saturation pressure at `RT` using Wilson correlation.

# Arguments
- `Pc::Real` - critical pressure
- `RTc::Real` - gas constant * critical temperature
- `acentric_factor::Real` - acentric factor
- `RT::Real` - gas constant * temperature

Reference: Brusilovsky2002 [p 272, eq 5.4], Mikyska2010 DOI 10.1002/aic.12387 [Algorithm step 1]
"""
function wilson_saturation_pressure(Pc::Real, RTc::Real, acentric_factor::Real, RT::Real)
    return Pc * exp(5.373 * (1.0 + acentric_factor) * (1.0 - RTc / RT))
end

"""
    eos_parameters(mixture::BrusilovskyEoSMixture, nmol, RT[; buf])

Return `Am`, `Bm`, `Cm`, `Dm` coefficients and `aij` matrix of EoS for `mixture` with
composition `nmol` and thermal energy `RT`. Allocations may be avoided by passing `buf`.

# Arguments
- `nmol::AbstractVector`: Vector of composition
- `RT::Real`: Thermal energy (J mol⁻¹)

# Keywords
- `buf`: see ?pressure(mixture)
"""
function eos_parameters(
    mixture::BrusilovskyEoSMixture,
    nmol::AbstractVector{<:Real},
    RT::Real;
    buf = NamedTuple(),
)
    return __eos_parameters_impl__(mixture, nmol, RT, buf)
end

# When `buf` is BrusilovskyThermoBuffer
function __eos_parameters_impl__(
    mixture::BrusilovskyEoSMixture{T},
    nmol::AbstractVector{<:Real},
    RT::Real,
    buf::BrusilovskyThermoBuffer{T},
) where {T}
    aij = buf.matr
    ai = buf.vec1
    return __eos_parameters_impl__(mixture, nmol, RT, ai, aij)
end

# When `buf` is mapping-like object
function __eos_parameters_impl__(
    mixture::BrusilovskyEoSMixture{T},
    nmol::AbstractVector{<:Real},
    RT::Real,
    buf::Union{NamedTuple, AbstractDict},
) where {T}
    nc = length(nmol)
    aij = haskey(buf, :aij) ? buf[:aij] : Matrix{T}(undef, nc, nc)
    ai = haskey(buf, :ai) ? buf[:ai] : Vector{T}(undef, nc)
    return __eos_parameters_impl__(mixture, nmol, RT, ai, aij)
end

# Core function
# See Brusilovsky2002 [p 187, eqs 4.159, 4.163-4.165] and [p 189, eq 4.168]
function __eos_parameters_impl__(
    mixture::BrusilovskyEoSMixture{T},
    nmol::AbstractVector{<:Real},
    RT::Real,
    auxv::AbstractVector{<:Real},
    auxm::AbstractMatrix{<:Real},
) where {T}
    ai, aij = auxv, auxm

    comp = mixture.components
    ai .= a_coef.(comp, RT)

    Bm = Cm = Dm = zero(T)
    @inbounds for i in eachindex(nmol, comp)
        Bm += nmol[i] * comp[i].b
        Cm += nmol[i] * comp[i].c
        Dm += nmol[i] * comp[i].d
    end

    temp = RT / GAS_CONSTANT_SI - 273.16
    eij, gij, hij = mixture.eij, mixture.gij, mixture.hij
    aij .= (one(T) .- (eij .+ temp .* (gij .+ temp .* hij))) .* sqrt.(ai .* ai')
    Am = dot(nmol, aij, nmol)  # Am = nmolᵀ aij nmol
    return Am, Bm, Cm, Dm, aij
end

"""
    pressure(mixture, nmol, volume, RT[; buf])

Return pressure (Pa) of `mixture` at given

- `nmol::AbstractVector`: composition (molar parts) or number of moles (mol)
- `volume::Real`: molar volume (m³ mol⁻¹) or volume (m³)
- `RT::Real`: thermal energy (J mol⁻¹)

Allocations may be avoided by passing `buf`.

# Keywords
- `buf::Union{BrusilovskyThermoBuffer,NamedTuple,AbstractDict}`: Buffer for intermediate
    calculations. In case of `NamedTuple` and `AbstractDict` `buf` should contain `buf[:ai]`
    `NC = ncomponents(mixture)` vector and `buf[:aij]` NCxNC matrix.

See also: [`thermo_buffer`](@ref)
"""
function pressure(
    mixture::BrusilovskyEoSMixture,
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf = NamedTuple(),
)
    Am, Bm, Cm, Dm, _ = eos_parameters(mixture, nmol, RT; buf = buf)
    return sum(nmol) * RT / (volume - Bm) - Am/((volume + Cm) * (volume + Dm))
end

"""
    compressibility(mixture, χ, P, RT[, phase='g'][; buf])

Compute compressibility (z-factor) of `mixture` in `phase` at given

- composition in moles `χ` (mol)
- pressure `P` (Pa)
- thermal energy `RT` (J mol⁻¹)

# Optional arguments

- `phase::AbstractChar='g'`: specifies phase of `mixture` (`'g'` for gas, `'l'` for liquid)

# Keywords:
- `buf::Union{BrusilovskyThermoBuffer,NamedTuple,AbstractDict}`: buffer for intermediate
    calculations (see [`pressure`](@ref))

See also: [`pressure`](@ref)
"""
function compressibility(
    mix::BrusilovskyEoSMixture,
    χ::AbstractVector,
    P,
    RT,
    phase::AbstractChar='g';
    buf = NamedTuple(),
)

    phase in ('g', 'l') || throw(DomainError(repr(phase), "Phase must be 'g' or 'l'"))

    ntotal = sum(χ)
    am, bm, cm, dm = eos_parameters(mix, χ, RT; buf = buf)
    prt = P / (RT * ntotal)
    am *= prt / (RT * ntotal)
    bm *= prt
    cm *= prt
    dm *= prt
    zroots = solve_cubic(
        1.0,
        cm + dm - bm - 1.0,
        am - bm * cm + cm * dm - bm * dm - dm - cm,
        -bm * cm * dm - cm * dm - am * bm,
    )

    if phase == 'g'
        return maximum(z for z in zroots if z > 0.0) # NaN is filtered
    else # phase == 'l'
        return minimum(z for z in zroots if z > 0.0)
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     comment_magics: false
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

# ## Initialize a distribution

using GEMPIC

@doc SumCosGaussian

# +
using Plots

k, α = 0.5, 0.1
nx, nv = 64, 128
xmin, xmax = 0, 2π/k
vmin, vmax = -6, 6
xg = LinRange(xmin, xmax, nx+1)[1:end-1]
vg = LinRange(vmin, vmax, nv+1)[1:end-1];

# -

df_1d1v = SumCosGaussian{1,1}( [[k]] , [α] , [[1.0]] , [[0.0]], [1.0])
f = zeros(Float64, (nx,nv))
for  j in eachindex(vg), i in eachindex(xg)
    f[i,j] = df_1d1v( xg[i], vg[j])
end
surface(xg, vg, f')

two_gaussians = SumCosGaussian{1,1}( [[k]] , [0.0] , [[0.2],[0.2]] , [[-1.0],[1.0]], [0.5, 0.5])
f = zeros(Float64, (nx,nv))
for  j in eachindex(vg), i in eachindex(xg)
    f[i,j] = two_gaussians( xg[i], vg[j])
end
surface(xg, vg, f')

# +
fxv(x, v) = ( 1 + α * cos(k * x )) / sqrt(2π) * exp(- (v^2)/ 2)

surface(xg, vg, fxv)
# -

@doc CosSumGaussian

df = CosSumGaussian{1,1}( [[k]], [0.1], [[1.0]], [[0.0]], [1.0] )

plot( xg, [GEMPIC.eval_x_density(df, x) for x in xg])

plot( vg, [GEMPIC.eval_v_density(df, v) for v in vg])


text:
#from https://github.com/JuliaText/Embeddings.jl/blob/master/test/runtests.jl
"""
    tempdatadeps(fun)
Run the function and delete all created datadeps afterwards
"""
function tempdatadeps(fun)
    tempdir = mktempdir()
    try
        @info "sending all datadeps to $tempdir"
        withenv("DATADEPS_LOAD_PATH"=>tempdir) do
            fun()
        end
    finally
        try
            @info "removing $tempdir"
            rm(tempdir, recursive=true, force=true)
        catch err
            @warn "Something went wrong with removing tempdir" tempdir exception=err
        end
    end
end

"""
@testset_nokeep_data
Use just like @testset,
but know that it deletes any downloaded data dependencies when it is done.
"""
macro testset_nokeep_data(name, expr)
    quote
        tempdatadeps() do
            @testset $name $expr
        end
    end |> esc
end

@testset "Pretrain" begin
  using Transformers.Pretrain
  using DataDeps
  using BytePairEncoding
  ENV["DATADEPS_ALWAYS_ACCEPT"] = true
  ENV["DATADEPS_NO_STANDARD_LOAD_PATH"] = true
  ENV["DATADEPS_PROGRESS_UPDATE_PERIOD"] = "Inf" #should be able remove one day
  ENV["DATADEP_PROGRESS_UPDATE_PERIOD"] = "Inf"

  @test_nowarn pretrains()
  model_list = map(l->join([l[1],l[2]], "-"), pretrains().content[1].rows[2:end])

  for model ∈ model_list
    @testset_nokeep_data "$model" begin
      GC.gc()
      if startswith(lowercase(model), "bert") && model[6:end] != "uncased_L-12_H-768_A-12"
        model_name = model[6:end]
        @test_nowarn Transformers.Pretrain.@datadep_str "BERT-$model_name/$model_name.tfbson"
      else
        @test_nowarn x = Transformers.Pretrain.@pretrain_str model
      end
    end
  end

end