text:
using LinearAlgebra: lmul!
using LinearAlgebra.BLAS: gemm!
using NNlib

"""

    conv4(w, x; kwargs...)

Execute convolutions or cross-correlations using filters specified
with `w` over tensor `x`.

Currently KnetArray{Float32/64,4/5} and Array{Float32/64,4} are
supported as `w` and `x`.  If `w` has dimensions `(W1,W2,...,I,O)` and
`x` has dimensions `(X1,X2,...,I,N)`, the result `y` will have
dimensions `(Y1,Y2,...,O,N)` where

    Yi=1+floor((Xi+2*padding[i]-Wi)/stride[i])

Here `I` is the number of input channels, `O` is the number of output
channels, `N` is the number of instances, and `Wi,Xi,Yi` are spatial
dimensions.  `padding` and `stride` are keyword arguments that can be
specified as a single number (in which case they apply to all
dimensions), or an array/tuple with entries for each spatial
dimension.

# Keywords

* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `dilation=1`: dilation factor for each dimension.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.

"""
function conv4(w::KnetArray{T},x::KnetArray{T}; handle=cudnnhandle(), alpha=1,
               o...) where {T} # padding=0, stride=1, dilation=1, mode=0
    beta=0 # nonzero beta does not make sense when we create y
    y = similar(x, cdims(w,x;o...))
    (algo,workSpace) = conv4_algo(w, x, y; handle=handle, o...)
    @cudnn(cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,   UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(T(alpha)),TD(x),x,FD(w),w,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),TD(y),y)
    return y
end

function conv4x(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T}; handle=cudnnhandle(), alpha=1,
                   o...) where {T} # padding=0, stride=1, dilation=1, mode=0
    beta = 0
    dx = similar(x)
    (algo,workSpace) = conv4x_algo(w,x,dy,dx; handle=handle, o...)
    if cudnnVersion >= 4000
        @cudnn(cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),TD(dx),dx)
    elseif cudnnVersion >= 3000
        @cudnn(cudnnConvolutionBackwardData_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),TD(dx),dx)
    else
        @cudnn(cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,       Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),Ref(T(beta)),TD(dx),dx)
    end
    return dx
end

function conv4w(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T}; handle=cudnnhandle(), alpha=1,
                   o...) where {T} # padding=0, stride=1, dilation=1, mode=0
    beta = 0
    dw = similar(w)
    (algo,workSpace) = conv4w_algo(w,x,dy,dw;handle=handle,o...)
    if cudnnVersion >= 4000
        @cudnn(cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),FD(dw),dw)
    elseif cudnnVersion >= 3000
        @cudnn(cudnnConvolutionBackwardFilter_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),FD(dw),dw)
    else
        @cudnn(cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,       Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),Ref(T(beta)),FD(dw),dw)
    end
    return dw
end

@primitive conv4(w,x; o...),dy  conv4w(w,x,dy;o...)  conv4x(w,x,dy;o...)
@zerograd  conv4x(w,x,dy;o...)
@zerograd  conv4w(w,x,dy;o...)

"""

    pool(x; kwargs...)

Compute pooling of input values (i.e., the maximum or average of
several adjacent values) to produce an output with smaller height
and/or width.

Currently 4 or 5 dimensional KnetArrays with `Float32` or `Float64`
entries are supported.  If `x` has dimensions `(X1,X2,...,I,N)`, the
result `y` will have dimensions `(Y1,Y2,...,I,N)` where

    Yi=1+floor((Xi+2*padding[i]-window[i])/stride[i])

Here `I` is the number of input channels, `N` is the number of
instances, and `Xi,Yi` are spatial dimensions.  `window`, `padding`
and `stride` are keyword arguments that can be specified as a single
number (in which case they apply to all dimensions), or an array/tuple
with entries for each spatial dimension.

# Keywords:

* `window=2`: the pooling window size for each dimension.
* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=window`: the number of elements to slide to reach the next pooling window.
* `mode=0`: 0 for max, 1 for average including padded values, 2 for average excluding padded values.
* `maxpoolingNanOpt=0`: Nan numbers are not propagated if 0, they are propagated if 1.
* `alpha=1`: can be used to scale the result.
* `handle`: Handle to a previously created cuDNN context. Defaults to a Knet allocated handle.

"""
function pool(x::KnetArray{T}; handle=cudnnhandle(), alpha=1,
                 o...) where {T} # window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    y = similar(x, pdims(x; o...))
    beta = 0
    @cudnn(cudnnPoolingForward,
          (Cptr, Cptr,      Ptr{T},    Cptr,Ptr{T},Ptr{T},   Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(T(alpha)),TD(x),x,    Ref(T(beta)),TD(y),y)
    return y
end

function poolx(x::KnetArray{T},y::KnetArray{T},dy::KnetArray{T}; handle=cudnnhandle(), alpha=1, mode=0,
                  o...) where {T} # window=2, padding=0, stride=window, maxpoolingNanOpt=0
    dx = similar(x)
    beta = 0
    @cudnn(cudnnPoolingBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,PD(x;mode=mode,o...),Ref(T(alpha)),TD(y),y,TD(dy),dy,TD(x),x,Ref(T(beta)),TD(dx),dx)
    return dx
end

@primitive pool(x;o...),dy,y  poolx(x,y,dy;o...)
@zerograd  poolx(x,y,dy;o...)

### CPU convolution using NNlib

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

function conv4(w::AbstractArray{T,N}, x::AbstractArray{T,N};
               padding=0, stride=1, dilation=1, mode=0, alpha=1) where {T,N}
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    y = conv(x, w, cdims)
    alpha == 1 ? y : lmul!(alpha, y)
end

function conv4w(w::AbstractArray{T,N},x::AbstractArray{T,N},dy::AbstractArray{T,N};
                padding=0, stride=1, dilation=1, mode=0, alpha=1) where {T,N}
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    dw = ∇conv_filter(x, dy, cdims)
    alpha == 1 ? dw : lmul!(alpha, dw)
end

function conv4x(w::AbstractArray{T,N},x::AbstractArray{T,N},dy::AbstractArray{T,N};
                padding=0, stride=1, dilation=1, mode=0, alpha=1) where {T,N}
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    dx = ∇conv_data(dy, w, cdims)
    alpha == 1 ? dx : lmul!(alpha, dx)
end

# TODO: handle alpha, maxpoolingNanOpt, mode
# TODO: pool(ka(x), stride=1) gives error

function pool(x::AbstractArray{T,N}; handle=nothing,
              alpha=1, mode=0, window=2, padding=0, stride=window, maxpoolingNanOpt=0) where {T,N}
    window = expand(Val(N-2), window)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    pdims = PoolDims(x, window; padding = padding, stride = stride)
    y = (mode == 0 ? maxpool(x, pdims) :
         mode == 1 ? meanpool(x, pdims) :
         # mode == 2 ? meanpool(x, pdims) : ## CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING is missing in NNlib Issue #218
         error("mode=$mode is not supported for CPU pool."))
    alpha == 1 ? y : lmul!(alpha, y)
end

function poolx(x::AbstractArray{T,N},y::AbstractArray{T,N},dy::AbstractArray{T,N}; handle=nothing,
               alpha=1, mode=0, window=2, padding=0, stride=window, maxpoolingNanOpt=0) where {T,N}
    if alpha != 1
        y = y ./ T(alpha)
    end
    window = expand(Val(N-2), window)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    pdims = PoolDims(x, window; padding = padding, stride = stride)
    dx = (mode == 0 ? ∇maxpool(dy, y, x, pdims) :
          mode == 1 ? ∇meanpool(dy, y, x, pdims) :
          # mode == 2 ? ∇meanpool(dy, y, x, pdims) : ## CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING is missing in NNlib Issue #218
          error("mode=$mode is not supported for CPU pool."))
    alpha == 1 ? dx : lmul!(alpha, dx)
end

"""

Unpooling; `reverse` of pooling.

TODO: Does not work correctly for every window, padding, mode combination. Test before use.

    x == pool(unpool(x;o...); o...)

"""
function unpool(x; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,x))
    y = similar(x,updims(x; window=window, o...))
    poolx(y,x,x.*w; o..., window=window, mode=1, alpha=1/alpha)
end

function unpoolx(dy; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,dy))
    pool(dy; o..., window=window, mode=1, alpha=1/alpha) * w
end

# @primitive unpool(x;o...),dy,y -pool(-dy;o...)
@primitive  unpool(x;o...),dy,y  unpoolx(dy;o...)

"""

    y = deconv4(w, x; kwargs...)

Simulate 4-D deconvolution by using _transposed convolution_ operation. Its forward pass is
equivalent to backward pass of a convolution (gradients with respect to input
tensor). Likewise, its backward pass (gradients with respect to input tensor) is equivalent
to forward pass of a convolution. Since it swaps forward and backward passes of convolution
operation, padding and stride options belong to output tensor. See [this
report](https://arxiv.org/abs/1603.07285) for further explanation.

Currently KnetArray{Float32/64,4} and Array{Float32/64,4} are supported as `w` and `x`.  If
`w` has dimensions `(W1,W2,...,O,I)` and `x` has dimensions `(X1,X2,...,I,N)`, the result
`y` will have dimensions `(Y1,Y2,...,O,N)` where

Yi = Wi+stride[i]*(Xi-1)-2*padding[i]

Here I is the number of input channels, O is the number of output channels, N is the number
of instances, and Wi,Xi,Yi are spatial dimensions. padding and stride are keyword arguments
that can be specified as a single number (in which case they apply to all dimensions), or an
array/tuple with entries for each spatial dimension.

# Keywords

* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.

"""
function deconv4(w,x; o...)
    y = similar(x,dcdims(w,x;o...))
    return conv4x(w,y,x;o...)
end

function deconv4w(w,x,dy; o...)
    return conv4w(w,dy,x;o...)
end

function deconv4x(w,x,dy; o...)
    return conv4(w,dy;o...)
end

@primitive deconv4(w,x; o...),dy,y  deconv4w(w,x,dy; o...)  deconv4x(w,x,dy; o...)
@zerograd deconv4w(w,x,dy; o...)
@zerograd deconv4x(w,x,dy; o...)

# cudnn descriptors

mutable struct TD; ptr; end
TD(a::KnetArray{T}) where {T} = TD(T,size(a))
TD(T::Type, dims::Integer...) = TD(T, dims)
function TD(T::Type, dims)
    d = Cptr[0]
    @cudnn(cudnnCreateTensorDescriptor,(Ptr{Cptr},),d)
    n = length(dims)
    sz = [Cint(dims[i]) for i=n:-1:1]
    st = similar(sz); st[n] = 1
    for i=(n-1):-1:1; st[i] = st[i+1] * sz[i+1]; end
    @cudnn(cudnnSetTensorNdDescriptor,
          (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint}),
          d[1], DT(T), n, sz, st)
    td = TD(d[1])
    finalizer(x->@cudnn(cudnnDestroyTensorDescriptor,(Cptr,),x.ptr), td)
    return td
end

mutable struct FD; ptr; end
FD(a::KnetArray{T}) where {T}=FD(T,size(a))
FD(T::Type, dims::Integer...) = FD(T,dims)
function FD(T::Type, dims)
    d = Cptr[0]
    @cudnn(cudnnCreateFilterDescriptor,(Ptr{Cptr},),d)
    n = length(dims)
    sz = [Cint(dims[i]) for i=n:-1:1]
    if cudnnVersion >= 5000
        @cudnn(cudnnSetFilterNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
              d[1], DT(T), 0,     n,   sz)
    elseif cudnnVersion >= 4000
        @cudnn(cudnnSetFilterNdDescriptor_v4,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
              d[1], DT(T), 0,     n,   sz)
    else
        @cudnn(cudnnSetFilterNdDescriptor,
              (Cptr,UInt32,Cint,Ptr{Cint}),
              d[1], DT(T),    n,   sz)
    end
    fd = FD(d[1])
    finalizer(x->@cudnn(cudnnDestroyFilterDescriptor,(Cptr,),x.ptr), fd)
    return fd
end

mutable struct CD; ptr
    function CD(w::KnetArray,x::KnetArray; padding=0, stride=1, dilation=1, mode=0, upscale=nothing)
        upscale !== nothing && error("upscale is deprecated, please use dilation instead.")
        d = Cptr[0]
        @cudnn(cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        if cudnnVersion >= 4000
            @cudnn(cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(dilation,nd),mode,DT(x))
        elseif cudnnVersion > 3000 # does not work when cudnnVersion==3000
            @cudnn(cudnnSetConvolutionNdDescriptor_v3,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(dilation,nd),mode,DT(x))
        else
            @cudnn(cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(dilation,nd),mode)
        end
        cd = new(d[1])
        finalizer(x->@cudnn(cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr),cd)
        return cd
    end
end

mutable struct PD; ptr
    function PD(x::KnetArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
        d = Cptr[0]
        @cudnn(cudnnCreatePoolingDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        if cudnnVersion >= 5000
            @cudnn(cudnnSetPoolingNdDescriptor,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        elseif cudnnVersion >= 4000
            @cudnn(cudnnSetPoolingNdDescriptor_v4,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        else
            @cudnn(cudnnSetPoolingNdDescriptor,
                  (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        end
        pd = new(d[1])
        finalizer(x->@cudnn(cudnnDestroyPoolingDescriptor,(Cptr,),x.ptr), pd)
        return pd
    end
end

import Base: unsafe_convert
unsafe_convert(::Type{Cptr}, td::TD)=td.ptr
unsafe_convert(::Type{Cptr}, fd::FD)=fd.ptr
unsafe_convert(::Type{Cptr}, cd::CD)=cd.ptr
unsafe_convert(::Type{Cptr}, pd::PD)=pd.ptr

# fill and reverse Cint array with padding etc. for cudnn calls
function cdsize(w, nd)
    if isa(w,Number)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch("$w $nd"))
    end
end

# convert padding etc. size to an Int array of the right dimension
function psize(p, x)
    nd = ndims(x)-2
    if isa(p,Number)
        fill(Int(p),nd)
    elseif length(p)==nd
        collect(Int,p)
    else
        throw(DimensionMismatch("psize: $p $nd"))
    end
end

DT(::KnetArray{Float32})=Cint(0)
DT(::KnetArray{Float64})=Cint(1)
DT(::KnetArray{Float16})=Cint(2)
DT(::Type{Float32}) = Cint(0)
DT(::Type{Float64}) = Cint(1)
DT(::Type{Float16}) = Cint(2)

# outputDim = 1 + ( inputDim + 2*pad - (((filterDim-1)*dilation)+1) )/convolutionStride;

function cdims(w,x; padding=0, stride=1, dilation=1, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            di = (if isa(dilation,Number); dilation; else dilation[i]; end)
            1 + (size(x,i) + 2*pi - (((size(w,i)-1)*di)+1)) ÷ si
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

function dcdims(w,x; padding=0, stride=1, dilation=1, o...)
    # TODO: handle dilation here
    dilation != 1 && error("deconv4 cannot handle dilation!=1 yet.")
    N = ndims(x)
    @assert size(x,N-1) == size(w,N)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            # di = (if isa(dilation,Number); dilation; else dilation[i]; end)
            si*(size(x,i)-1) + size(w,i) - 2*pi
        elseif i == N-1
            size(w,N-1)
        else
            size(x,N)
        end
    end
end

function pdims(x; window=2, padding=0, stride=window, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            wi = (if isa(window,Number); window; else window[i]; end)
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) + 2*pi - wi, si)
        else
            size(x,i)
        end
    end
end

function updims(x; window=2, padding=0, stride=window, o...)
    window = psize(window,x)
    stride = psize(stride,x)
    padding = psize(padding,x)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            (size(x,i)-1)*stride[i]+window[i]-2*padding[i]
        else
            size(x,i)
        end
    end
end

# convolution padding size that preserves the input size when filter size is odd and stride=1
padsize(w)=ntuple(i->div(size(w,i)-1,2), ndims(w)-2)

## Utilities to find a fast algorithm

struct cudnnConvolutionFwdAlgoPerf_t
    algo::Cint
    status::Cint
    time::Cfloat
    memory::Csize_t
    determinism::Cint
    mathType::Cint
    r1::Cint; r2::Cint; r3::Cint
end

const CUDNN_MAX_FIND = 100      # How many times can we call FindAlgorithm
const requestedAlgoCount = 10
const returnedAlgoCount = Cint[0]
const perfResults = Array{cudnnConvolutionFwdAlgoPerf_t}(undef,requestedAlgoCount)
bytes(x::KnetArray{T}) where {T}=length(x)*sizeof(T)

const conv4_algos = Dict()
function conv4_algo(w::KnetArray{T}, x::KnetArray{T}, y::KnetArray{T}; handle=cudnnhandle(), o...) where {T}
    global conv4_algos, requestedAlgoCount, returnedAlgoCount, perfResults
    key = (T,size(w),size(x),o...)
    if haskey(conv4_algos, key)
        p = conv4_algos[key]
        return (p.algo, cudnnWorkSpace(p.memory))
    elseif length(conv4_algos) >= CUDNN_MAX_FIND
        return (0, cudnnWorkSpace())
    else
        # Knet.gc(); @dbg print('*')
        # TODO: not sure how much space to allocate here without risking OOM, will figure out later.
        # workSpace, workSpaceSizeInBytes = C_NULL, 0
        workSpace = KnetArray{UInt8}(undef, gpufree() ÷ 5) # Try with max %20 of available memory
        wd, xd, yd, cd = FD(w), TD(x), TD(y), CD(w,x;o...)
        @cudnn(cudnnFindConvolutionForwardAlgorithmEx,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr,Cptr,Csize_t),
              handle,xd,x,wd,w,cd,yd,y,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        workSpace = nothing; Knet.gc(); GC.gc()
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4_algos[key] = p
        return (p.algo, cudnnWorkSpace(p.memory))
    end
end

const conv4w_algos = Dict()
function conv4w_algo(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T},dw::KnetArray{T}; handle=cudnnhandle(), o...) where {T}
    global conv4w_algos, requestedAlgoCount, returnedAlgoCount, perfResults
    key = (T,size(w),size(x),o...)
    if haskey(conv4w_algos, key)
        p = conv4w_algos[key]
        return (p.algo, cudnnWorkSpace(p.memory))
    elseif length(conv4w_algos) >= CUDNN_MAX_FIND
        return (0, cudnnWorkSpace())
    else
        # Knet.gc(); @dbg print('*')
        # workSpace, workSpaceSizeInBytes = C_NULL, 0
        workSpace = KnetArray{UInt8}(undef, gpufree() ÷ 5) # Try with max %20 of available memory
        wd, xd, yd, cd = FD(dw), TD(x), TD(dy), CD(w,x;o...)
        @cudnn(cudnnFindConvolutionBackwardFilterAlgorithmEx,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr,Cptr,Csize_t),
              handle,xd,x,yd,dy,cd,wd,dw,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        workSpace = nothing; Knet.gc(); GC.gc()
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4w_algos[key] = p
        return (p.algo, cudnnWorkSpace(p.memory))
    end
end

const conv4x_algos = Dict()
function conv4x_algo(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T},dx::KnetArray{T}; handle=cudnnhandle(), o...) where {T}
    global conv4x_algos, requestedAlgoCount, returnedAlgoCount, perfResults
    key = (T,size(w),size(x),o...)
    if haskey(conv4x_algos, key)
        p = conv4x_algos[key]
        return (p.algo, cudnnWorkSpace(p.memory))
    elseif length(conv4x_algos) >= CUDNN_MAX_FIND
        return (0, cudnnWorkSpace())
    else
        # Knet.gc(); @dbg print('*')
        # workSpace, workSpaceSizeInBytes = C_NULL, 0
        workSpace = KnetArray{UInt8}(undef, gpufree() ÷ 5) # Try with max %20 of available memory
        wd, xd, yd, cd = FD(w), TD(dx), TD(dy), CD(w,x;o...)
        @cudnn(cudnnFindConvolutionBackwardDataAlgorithmEx,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr,Cptr,Csize_t),
              handle,wd,w,yd,dy,cd,xd,dx,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        workSpace = nothing; Knet.gc(); GC.gc()
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4x_algos[key] = p
        return (p.algo, cudnnWorkSpace(p.memory))
    end
end

CUDNN_WORKSPACE_MAXSIZE = 0     # Will be set to 20% of available gpu memory
function perfChoose(ps, n)
    global CUDNN_WORKSPACE_MAXSIZE
    if n==ps
        warn("returnedAlgoCount==requestedAlgoCount")
    end
    if CUDNN_WORKSPACE_MAXSIZE == 0
        CUDNN_WORKSPACE_MAXSIZE = gpufree()÷5
    end
    (ibest,mbest,tbest) = (0,CUDNN_WORKSPACE_MAXSIZE,Inf)
    for i = 1:n
        if ps[i].status == 0 && ps[i].memory < mbest && ps[i].time < tbest * 1.1
            (ibest,mbest,tbest) = (i,ps[i].memory,ps[i].time)
        end
    end
    if ibest == 0; error("No good algo found."); end
    return ps[ibest]
end

# This assumes one workspace per gpu. TODO: What about streams/contexts/handles etc?
const CUDNN_WORKSPACE = []
function cudnnWorkSpace_fast_but_risky(len=0;dev=gpu())
    global CUDNN_WORKSPACE
    if dev==-1; error("No cudnnWorkSpace for CPU"); end
    i = dev+2
    if isempty(CUDNN_WORKSPACE); resize!(CUDNN_WORKSPACE,gpuCount()+1); end
    if !isassigned(CUDNN_WORKSPACE,i) || length(CUDNN_WORKSPACE[i]) < len
        CUDNN_WORKSPACE[i]=KnetArray{UInt8}(undef,len);
    end
    return CUDNN_WORKSPACE[i]
end

# Fresh workspace for every op is safer:
cudnnWorkSpace(len=0)=KnetArray{UInt8}(undef,len)


text:
workspace()

@everywhere include("EratosthenesExampleSoftware.jl")

module Temp

using Eratosthenes # Bring in the top-level sim environment.

using ..EratosthenesExampleSoftware

scenario = setup("scenarios/basic-cube.yaml")
results = mc(scenario, 5)

display(results)

end

ext: jl


text:
# This script generates catalogues
# of objects in disk galaxies

println('\n', " "^4, "> Loading the packages...")

using Distributions
using LaTeXStrings
using Plots
using Random

# Use the PGFPlotsX backend for plots
# (this `@eval` trick is necessary to fool `DaemonMode`)
@eval Main begin
    using Plots
    pgfplotsx()
end

# Change some of the default parameters for plots
default(fontfamily = "Computer Modern", dpi = 300)

# Define the paths to output directories
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
PLOTS_DIR = joinpath(ROOT_DIR, "plots")

# Make sure the needed directories exist
mkpath(PLOTS_DIR)

# Define the number of objects in a catalogue
N = 1000

# Define probability density functions for some of the distributions
normal_pdf(x, μ, σ) = 1 / (σ * √(2 * π)) * exp(-0.5 * ((x - μ) / σ)^2)
uniform_pdf(x, a, b) = a < x < b ? 1 / (b - a) : 0
uniform_radius_pdf(x, r) = 0 < x < r ? 2 * x / r^2 : 0
exponential_pdf(x, b) = x >= 0 ? 1 / b * exp(-x / b) : 0
exponential_radius_pdf(r, b) = r >= 0 ? 2 * π / b * r * exp(-π / b * r^2) : 0
laplace_pdf(x, μ, b) = 1 / (2 * b) * exp(-abs(x - μ) / b)

"Return a range of ticks that are dividable by π / denom, along with the labels"
function piticks(start, stop, denom)
    # Compute the integer bounds for a range
    # of values dividable by π / denom
    a = Int(cld(start, π / denom))
    b = Int(fld(stop, π / denom))
    # Compute the range
    ticks = range(a * π / denom, b * π / denom; step = π / denom)
    # Define the labels for each tick
    ticklabels = piticklabel.((a:b) .// denom)
    return ticks, ticklabels
end

"Return a label for a rational number which is a multiple of π"
function piticklabel(x::Rational)
    # Return zero if it's a zero
    iszero(x) && return L"0"
    # Define the sign of the rational
    sign = x < 0 ? "-" : ""
    # Get the numerator and the denominator
    n, d = abs(numerator(x)), denominator(x)
    # Add a multiple to π if necessary
    N = n == 1 ? "" : repr(n)
    # Return the label
    d == 1 ? L"%$sign%$N\pi" : L"%$sign\frac{%$N\pi}{%$d}"
end

function plot_distribution(sample::AbstractVector, pdf::Function, dir::AbstractString)
    println(" "^4, "> Plotting the $dir distribution...")

    # Define the output directory
    output_dir = joinpath(PLOTS_DIR, dir)

    # Make sure the directory exists
    mkpath(output_dir)

    # Plot the histogram of the distribution
    histogram(
        sample;
        label = "",
        title = "",
        xlabel = L"x",
        ylabel = L"f(x)",
        normalize = true,
        size = (400, 400),
        yminorticks = 5,
    )

    # Plot the probability function over the histogram
    plot!(pdf, label = "", lw = 1.5)

    # Save the figure as PDF and PNG
    savefig(joinpath(output_dir, "Histogram.pdf"))
    savefig(joinpath(output_dir, "Histogram.png"))
end

# Set a seed for the random number generator
Random.seed!(2)

# Plot a sample from the normal distribution with μ=1 and σ=2
plot_distribution(rand(Normal(1, 2), N), (x) -> normal_pdf(x, 1, 2), "normal")

# Plot a sample from the uniform distribution with a=1 and b=2
plot_distribution(rand(Uniform(1, 2), N), (x) -> uniform_pdf(x, 1, 2), "uniform")

# Plot a sample from the exponential distribution with b=0.5
plot_distribution(rand(Exponential(0.5), N), (x) -> exponential_pdf(x, 0.5), "exponential")

# Plot a sample from the Laplace distribution with μ = 0 and b=0.5
plot_distribution(rand(Laplace(0, 0.5), N), (x) -> laplace_pdf(x, 0, 0.5), "laplace")

"Plot a histogram of the radius sample"
function plot_histogram_radius(r, pdf::Function, dir)
    # Plot the histogram of the radius sample
    histogram(
        r;
        label = "",
        title = "",
        xlabel = L"r",
        ylabel = L"f(r)",
        normalize = true,
        size = (400, 400),
        yminorticks = 5,
    )

    # Plot the probability function over the histogram
    plot!(pdf, label = "", lw = 1.5)

    # Save the figure as PDF and PNG
    savefig(joinpath(dir, "Histogram R.pdf"))
    savefig(joinpath(dir, "Histogram R.png"))

end

"Plot a histogram of the angle sample"
function plot_histogram_angle(θ, dir)
    # Plot the histogram of the angle sample
    histogram(
        θ;
        label = "",
        title = "",
        xlabel = L"\theta",
        ylabel = L"f(\theta)",
        normalize = true,
        bins = range(0, stop = 2 * π, step = π / 4),
        size = (400, 400),
        xtick = piticks(0, 2 * π, 4),
        yminorticks = 5,
    )

    # Plot the probability function over the histogram
    plot!((x) -> uniform_pdf(x, 0, 2 * π), label = "", lw = 1.5)

    # Save the figure as PDF and PNG
    savefig(joinpath(dir, "Histogram θ.pdf"))
    savefig(joinpath(dir, "Histogram θ.png"))
end

"Plot a histogram of the height sample"
function plot_histogram_height(z, b, dir)
    # Plot the histogram of the height sample
    histogram(
        z;
        label = "",
        title = "",
        xlabel = L"z",
        ylabel = L"f(z)",
        normalize = true,
        size = (400, 400),
        yminorticks = 5,
    )

    # Plot the probability function over the histogram
    plot!((x) -> laplace_pdf(x, 0, b), label = "", lw = 1.5)

    # Save the figure as PDF and PNG
    savefig(joinpath(dir, "Histogram Z.pdf"))
    savefig(joinpath(dir, "Histogram Z.png"))
end

"Plot a heatmap of the symmetry plane"
function plot_heatmap_symmetry_plane(r, θ, dir; rₘ = 0.)
    # Switch to GR for heatmaps
    # (see https://github.com/JuliaPlots/Plots.jl/issues/3266)
    gr()

    # Plot the heatmap of the symmetry plane
    histogram2d(
        θ,
        r;
        proj = :polar,
        label = "",
        title = "",
        xlabel = L"\theta",
        yaxis = false,
        normalize = :probability,
        show_empty_bins = true,
        size = (525, 400),
        right_margin = 7.5 * Plots.mm,
        colorbar = :best,
        legend = false,
        nbins = 15,
    )

    # Add a radius error to show the scale
    θ₀ = 59 * π / 32
    if rₘ == 0
      rₘ = maximum(r)
    end
    xₘ = 1.1 * cos(θ₀ + π / 64)
    yₘ = 1.1 * sin(θ₀ + π / 64)
    plot!([0, θ₀], [0, rₘ], lw = 1.5, color = palette(:default)[3], arrow = true)
    annotate!(xₘ, yₘ, ("$(round(rₘ, digits = 4))", 10))

    # Save the figure as PNG
    savefig(joinpath(dir, "Heatmap Symmetry Plane.png"))

    # Switch back to PGFPlotsX
    pgfplotsx()
end

"Plot a heatmap of the vertical cut"
function plot_heatmap_vertical_cut(r, z, dir)
    # Plot the heatmap of a vertical cut
    histogram2d(
        r,
        z;
        label = "",
        title = "",
        xlabel = L"r",
        ylabel = L"z",
        normalize = :probability,
        size = (400, 400),
        nbins = 40,
    )

    # Save the figure as PDF and PNG
    savefig(joinpath(dir, "Heatmap Vertical Cut.pdf"))
    savefig(joinpath(dir, "Heatmap Vertical Cut.png"))
end

"""
Generate a catalogue for a disk with uniform distribution in
the symmetry plane and Laplace distribution on the vertical axis
"""
function plot_uniform_disk(seed, R::F = 2.0, b::F = 0.5) where {F <: Real}
    println(" "^4, "> Plotting the uniform disk (seed = $seed) distribution...")

    # Set a seed for the random generator
    Random.seed!(seed)

    # Define the output directory
    output_dir = joinpath(PLOTS_DIR, "uniform_disk", "$seed")

    # Make sure the directory exists
    mkpath(output_dir)

    # Generate samples from a uniformly distributed
    # circle with exponentially distributed heights
    # (see https://stackoverflow.com/a/50746409 and
    # https://stats.stackexchange.com/a/378957 for
    # explanation)
    r = sqrt.(rand(Uniform(0, 1), N)) * R
    θ = rand(Uniform(0, 2 * π), N)
    z = rand(Laplace(0, b), N)

    # Plot the histogram of the radius sample
    plot_histogram_radius(r, (x) -> uniform_radius_pdf(x, R), output_dir)

    # Plot the histogram of the angle sample
    plot_histogram_angle(θ, output_dir)

    # Plot the histogram of the height sample
    plot_histogram_height(z, b, output_dir)

    # Plot the heatmap of the symmetry plane
    plot_heatmap_symmetry_plane(r, θ, output_dir, rₘ = R)

    # Plot the heatmap of the vertical cut
    plot_heatmap_vertical_cut(r, z, output_dir)
end

# Plot samples from the uniform disk distribution
plot_uniform_disk(2)

"""
Generate a catalogue for a disk with exponential distribution in
the symmetry plane and Laplace distribution on the vertical axis
"""
function plot_exponential_disk(seed, b₁::F = 2.6, b₂::F = 0.9) where {F <: Real}
    println(" "^4, "> Plotting the exponential disk (seed = $seed) distribution...")

    # Set a seed for the random generator
    Random.seed!(seed)

    # Define the output directory
    output_dir = joinpath(PLOTS_DIR, "exponential_disk", "$seed")

    # Make sure the directory exists
    mkpath(output_dir)

    # Generate samples from an exponentially distributed
    # circle with exponentially distributed heights
    # (see https://stackoverflow.com/a/50746409 and
    # https://stats.stackexchange.com/a/378957 for
    # explanation)
    r = sqrt.(-b₁ / π .* log.(1 .- rand(Uniform(0, 1), N)))
    θ = rand(Uniform(0, 2 * π), N)
    z = rand(Laplace(0, b₂), N)

    # Plot the histogram of the radius sample
    plot_histogram_radius(r, (x) -> exponential_radius_pdf(x, b₁), output_dir)

    # Save the figure as PDF and PNG
    savefig(joinpath(output_dir, "Histogram R.pdf"))
    savefig(joinpath(output_dir, "Histogram R.png"))

    # Plot the histogram of the angle sample
    plot_histogram_angle(θ, output_dir)

    # Plot the histogram of the height sample
    plot_histogram_height(z, b₂, output_dir)

    # Plot the heatmap of the symmetry plane
    plot_heatmap_symmetry_plane(r, θ, output_dir)

    # Plot the heatmap of the vertical cut
    plot_heatmap_vertical_cut(r, z, output_dir)
end

# Plot samples from the exponential disk distribution
plot_exponential_disk(1)
plot_exponential_disk(2)
plot_exponential_disk(3)

println()


text:
using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff

import LinearAlgebra

# n = Int(1e1)
n = Int(1e4)
k = 5 * Int(1e3)
number_nonzero = 40

xpi = rand(n);
total = sum(xpi);
const xp = xpi # ./ total;

f(x) = LinearAlgebra.norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

lmo = FrankWolfe.KSparseLMO(number_nonzero, 1.0);
## alternative lmo
# lmo = FrankWolfe.ProbabilitySimplexOracle(1)
x0 = FrankWolfe.compute_extreme_point(lmo, ones(n));

@time x, v, primal, dual_gap, trajectorylazy, active_set = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    epsilon=1e-5,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    lazy=true,
);

@time x, v, primal, dual_gap, trajectoryAFW, active_set = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    epsilon=1e-5,
    emphasis=FrankWolfe.memory,
    verbose=true,
    away_steps=true,
    trajectory=true,
);

@time x, v, primal, dual_gap, trajectoryFW = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    epsilon=1e-5,
    trajectory=true,
    away_steps=false,
);

data = [trajectorylazy, trajectoryAFW, trajectoryFW]
label = ["LAFW" "AFW" "FW"]

FrankWolfe.plot_trajectories(data, label, xscalelog=true)


text:
############################################
#
#    dolfin.cpp.geometry
# https://fenicsproject.org/docs/dolfin/2018.1.0/python/_autogenerated/dolfin.cpp.geometry.html
#
#    ufl.geometry
#
############################################

@pyclass dolfin Point
@pyclass dolfin BoundingBoxTree
@pyclass dolfin MeshPointIntersection

@pyclass ufl FacetNormal Expression
@pyclass ufl CellDiameter Expression

@pyclass ufl CellNormal
@pyclass ufl CellVolume
@pyclass ufl Circumradius
@pyclass ufl MaxCellEdgeLength
@pyclass ufl MaxFacetEdgeLength
@pyclass ufl MinCellEdgeLength
@pyclass ufl MinFacetEdgeLength
@pyclass ufl SpatialCoordinate

export Point, BoundingBoxTree, MeshPointIntersection, FacetNormal, CellDiameter,
       CellNormal, CellVolume, Circumradius, MaxCellEdgeLength, MaxFacetEdgeLength,
       MinCellEdgeLength, MinFacetEdgeLength, SpatialCoordinate


text:
module Thermic

"""
    run_thermic( dem_file::AbstractString, source_file::AbstractString, river_file::AbstractString, source_temperature::Real, source_flow_rate::Real, river_temperature::Real, river_flow_rate::Real, output_path::AbstractString=".\\")
"""
function run_thermic( dem_file::AbstractString, source_file::AbstractString, river_file::AbstractString, source_temperature::Real, source_flow_rate::Real, river_temperature::Real, river_flow_rate::Real, output_path::AbstractString=".\\")

    src_geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(src_geom) != 0
        throw(DomainError(source_file, "The shapefile must contain a point."))
    end

    river_layer = agd.getlayer(agd.read(river_file), 0)

 # IL CONTROLLO SULLA GEOMETRIA PUO' ESSERE FATTO SOLO SU FEATURES
    if agd.geomdim(river_layer) != 1
        throw(DomainError(river_file, "The river shapefile geometry in not valid."))
    end

    refsys = agd.getspatialref(source)

    dem = agd.read(dem_file)

    if agd.importWKT(agd.getproj(dem)) !=  refsys ||  agd.getspatialref(layer) != refsys
        throw(DomainError("The reference systems are not uniform. Aborting analysis." ))
    end

    demband = agd.getband(dem, 1)

    #   start_time = time.time()

 """ CREDO STIA CONSIDERANDO UN'AREA DI 5 METRI(?) INTORNO ALLA SORGENTE
    buffer_sorgente = processing.run("native:buffer", {"INPUT": self.source, "DISTANCE": 5, "OUTPUT":"memory:"})
    sorgente = buffer_sorgente["OUTPUT"]
 """
    # Consider a 5m area around the source
    source_area = agd.buffer(source, 5.0)
    # Find the portion of the river that intersects said area; there could be no intersection
    for feature in layer
        segment = agd.getgeom(feature)
        if agd.intersects(source_area, segment)
            intersection = agd.intersection(source_area, segment)

         """ IN QUESTA PARTE SETTA I VALORI DI TEMPERATURA E PORTATA DEL FIUME A QUELLI DELLA FUNZIONANTE
            nuova_portata = { a.fieldNameIndex(self.text_river_q): f[self.text_source_q]}
            nuova_temp = { a.fieldNameIndex(self.text_river_t): f[self.text_source_t]}
            self.river.dataProvider().changeAttributeValues({a.id(): nuova_portata})
            self.river.dataProvider().changeAttributeValues({a.id(): nuova_temp})
         """
            # BISOGNA VEDERE COME GESTIRE L'INPUT DI source E river
            new_flow = Dict("Portata" => source_flow_rate)
            new_temperature = Dict("Temperatura" => source_temperature)
            flow_idx, temp_idx = agd.findfieldindex( Ref(feature), [:portata, :temperatura] )
            agd.setfield!.( Ref(feature), [flow_idx, temp_idx], [new_flow, new_temperature] )
            break
        end
    end

 """ STA PRENDENDO L'INTERSEZIONE CON river DI QUALCOSA, MA NON SO COSA """
 # FORSE STA PRENDENDO LE CELLE DAL RASTER CHE CORRISPONDONO A river
    param_point_intersect_all = Dict( :Input_Fields => [], :Output => :memory, :Intersect => river, :Intersect_Fields => [], :Input => river )
    point_intersect_all_result = processing.run( "native:lineintersections", param_point_intersect_all )
    point_intersect_all = point_intersect_all_result[:Output]

    point_intersect_result = processing.run( "qgis:deleteduplicategeometries", Dict( :Output => :memory, :Input => point_intersect_all ) )
    point_intersect = point_intersect_result[:Output]
 """                                                                 """
    # NON E' CORRETTO, BISOGNA VEDERE COME OTTENIAMO point_intersect
    features = agd.getfeatures(point_intersect)

    count = 1
    frmain = 0
    dict = Dict()
    for feature in features
        geom = agd.getgeom(feature)
        x = agd.getx(geom, 0)
        y = agd.gety(geom, 0)
        r, c = toIndexes(dtm, x, y)
        z = demband[r, c]

        res = agd.getfield.( Ref(feature), [:portata_river, :temperatura_river, :portata_river_2, :temperatura_river_2] )
        if z in keys(dict)
            dict[z] = res
        else
            push!( dict, z => res )
        end
    end

    flow_rate = 0.0
    temperature = 0.0
    for key in keys(dict)
        if count == 1
            if dict[key][1] >= dict[key][3]
                frmain = dict[key][1]
                flow_rate = dict[key][1]
                temperature = dict[key][2]
            else
                frmain = dict[key][3]
                flow_rate = diz[key][3]
                temperature = diz[key][4]
            end
        end
        if dict[key][1] == frmain
            flow_rate1 = flow_rate
            temperature1 = temperature
            flow_rate2 = dict[key][3]
            temperature2 = dict[key][4]
        else
            flow_rate1 = flow_rate
            temperature1 = temperature
            flow_rate2 = dict[key][1]
            temperature2 = dict[key][2]
        end

        flow_rate = flow_rate1 + flow_rate2
        temperature = ( ( flow_rate1 * temperature1 ) + ( flow_rate2 * temperature2 ) ) / ( flow_rate1 + flow_rate2 )
        count += 1
    end

    #    tempoanalisi = time.time() - start_time
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
using Pkg
Pkg.add("JSON")
Pkg.add("DataFrames")
Pkg.add("CSV")
using JSON
using DataFrames
using CSV

for file in readdir()
    if endswith(file, ".json")
        f = JSON.parse(read(open(file, "r"), String))
        df = DataFrame(f[1])
            for dict in f[2:end]
                push!(df, dict)
            end
        println(df)
        name = replace(file, r".json" => s".csv")
        CSV.write(name, df)
    end
end


text:
function make_logpdf(π::AbstractStatelessPolicy)
    return (a)->logpdf(π, a)
end

function make_logpdf(π::AbstractPolicy)
    return (s,a)->logpdf(π, s, a)
end

function ope(π, stat, ism::AbstractImportanceSampling; α=0.0)

    if α > 0.0
        J = (D, π)->stat(estimate_returns_withentropy(ism, π, D, α))
    else
        J = (D, π)->stat(estimate_returns(ism, π, D))
    end

    return J
end


text:
###
### **acb_hypgeom.h** -- hypergeometric functions of complex variables
###

### Convergent series
arbcall"void acb_hypgeom_pfq_bound_factor(mag_t C, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, ulong n)"
arbcall"slong acb_hypgeom_pfq_choose_n(acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_fme(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t w, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_sum_invz(acb_t s, acb_t t, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, const acb_t w, slong n, slong prec)"
arbcall"void acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)"
#ni arbcall"void acb_hypgeom_pfq_series_sum_forward(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)"
#ni arbcall"void acb_hypgeom_pfq_series_sum_bs(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)"
#ni arbcall"void acb_hypgeom_pfq_series_sum_rs(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)"
#ni arbcall"void acb_hypgeom_pfq_series_sum(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)"
#ni arbcall"void acb_hypgeom_pfq_series_direct(acb_poly_t res, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)"

### Asymptotic series
arbcall"void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong n, slong prec)"
arbcall"int acb_hypgeom_u_use_asymp(const acb_t z, slong prec)"

### Generalized hypergeometric function
arbcall"void acb_hypgeom_pfq(acb_poly_t res, acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, int regularized, slong prec)"

### Confluent hypergeometric functions
arbcall"void acb_hypgeom_u_1f1_series(acb_poly_t res, const acb_poly_t a, const acb_poly_t b, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_m_asymp(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_m_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_0f1_asymp(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_0f1_direct(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_0f1(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec)"

### Error functions and Fresnel integrals
arbcall"void acb_hypgeom_erf_propagated_error(mag_t re, mag_t im, const acb_t z)"
arbcall"void acb_hypgeom_erf_1f1a(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_erf_1f1b(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_erf_asymp(acb_t res, const acb_t z, int complementary, slong prec, slong prec2)"
arbcall"void acb_hypgeom_erf(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_erf_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_erf_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_erfc(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_erfc_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_erfc_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_erfi(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_erfi_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_erfi_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_fresnel(acb_t res1, acb_t res2, const acb_t z, int normalized, slong prec)"
arbcall"void _acb_hypgeom_fresnel_series(acb_ptr res1, acb_ptr res2, acb_srcptr z, slong zlen, int normalized, slong len, slong prec)"
arbcall"void acb_hypgeom_fresnel_series(acb_poly_t res1, acb_poly_t res2, const acb_poly_t z, int normalized, slong len, slong prec)"

### Bessel functions
arbcall"void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_y(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_jy(acb_t res1, acb_t res2, const acb_t nu, const acb_t z, slong prec)"

### Modified Bessel functions
arbcall"void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)"
arbcall"void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)"
arbcall"void acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)"
arbcall"void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t nu, const acb_poly_t z, int scaled, slong len, slong prec)"
arbcall"void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)"
arbcall"void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)"

### Airy functions
arbcall"void acb_hypgeom_airy_direct(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_airy_asymp(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong n, slong prec)"
arbcall"void acb_hypgeom_airy_bound(mag_t ai, mag_t ai_prime, mag_t bi, mag_t bi_prime, const acb_t z)"
arbcall"void acb_hypgeom_airy(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_airy_jet(acb_ptr ai, acb_ptr bi, const acb_t z, slong len, slong prec)"
arbcall"void _acb_hypgeom_airy_series(acb_ptr ai, acb_ptr ai_prime, acb_ptr bi, acb_ptr bi_prime, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_airy_series(acb_poly_t ai, acb_poly_t ai_prime, acb_poly_t bi, acb_poly_t bi_prime, const acb_poly_t z, slong len, slong prec)"

### Coulomb wave functions
arbcall"void acb_hypgeom_coulomb(acb_t F, acb_t G, acb_t Hpos, acb_t Hneg, const acb_t l, const acb_t eta, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)"
arbcall"void _acb_hypgeom_coulomb_series(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_coulomb_series(acb_poly_t F, acb_poly_t G, acb_poly_t Hpos, acb_poly_t Hneg, const acb_t l, const acb_t eta, const acb_poly_t z, slong len, slong prec)"

### Incomplete gamma and beta functions
arbcall"void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_gamma_upper_singular(acb_t res, slong s, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)"
arbcall"void _acb_hypgeom_gamma_upper_series(acb_ptr res, const acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)"
arbcall"void acb_hypgeom_gamma_upper_series(acb_poly_t res, const acb_t s, const acb_poly_t z, int regularized, slong n, slong prec)"
arbcall"void acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)"
arbcall"void _acb_hypgeom_gamma_lower_series(acb_ptr res, const acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)"
arbcall"void acb_hypgeom_gamma_lower_series(acb_poly_t res, const acb_t s, const acb_poly_t z, int regularized, slong n, slong prec)"
arbcall"void acb_hypgeom_beta_lower(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)"
arbcall"void _acb_hypgeom_beta_lower_series(acb_ptr res, const acb_t a, const acb_t b, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)"
arbcall"void acb_hypgeom_beta_lower_series(acb_poly_t res, const acb_t a, const acb_t b, const acb_poly_t z, int regularized, slong n, slong prec)"

### Exponential and trigonometric integrals
arbcall"void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_ei_asymp(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_ei_2f2(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_ei(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_ei_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_ei_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_si_asymp(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_si_1f2(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_si(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_si_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_si_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_ci_asymp(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_ci(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_ci_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_ci_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_shi(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_shi_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_shi_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_chi_asymp(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_chi(acb_t res, const acb_t z, slong prec)"
arbcall"void _acb_hypgeom_chi_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_hypgeom_chi_series(acb_poly_t res, const acb_poly_t z, slong len, slong prec)"
arbcall"void acb_hypgeom_li(acb_t res, const acb_t z, int offset, slong prec)"
arbcall"void _acb_hypgeom_li_series(acb_ptr res, acb_srcptr z, slong zlen, int offset, slong len, slong prec)"
arbcall"void acb_hypgeom_li_series(acb_poly_t res, const acb_poly_t z, int offset, slong len, slong prec)"

### Gauss hypergeometric function
arbcall"void acb_hypgeom_2f1_continuation(acb_t res0, acb_t res1, const acb_t a, const acb_t b, const acb_t c, const acb_t z0, const acb_t z1, const acb_t f0, const acb_t f1, slong prec)"
arbcall"void acb_hypgeom_2f1_series_direct(acb_poly_t res, const acb_poly_t a, const acb_poly_t b, const acb_poly_t c, const acb_poly_t z, int regularized, slong len, slong prec)"
arbcall"void acb_hypgeom_2f1_direct(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, slong prec)"
arbcall"void acb_hypgeom_2f1_transform(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int flags, int which, slong prec)"
arbcall"void acb_hypgeom_2f1_transform_limit(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, int which, slong prec)"
arbcall"void acb_hypgeom_2f1_corner(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, slong prec)"
arbcall"int acb_hypgeom_2f1_choose(const acb_t z)"
arbcall"void acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int flags, slong prec)"

### Orthogonal polynomials and functions
arbcall"void acb_hypgeom_chebyshev_t(acb_t res, const acb_t n, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_chebyshev_u(acb_t res, const acb_t n, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_jacobi_p(acb_t res, const acb_t n, const acb_t a, const acb_t b, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_gegenbauer_c(acb_t res, const acb_t n, const acb_t m, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_laguerre_l(acb_t res, const acb_t n, const acb_t m, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_hermite_h(acb_t res, const acb_t n, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_legendre_p(acb_t res, const acb_t n, const acb_t m, const acb_t z, int type, slong prec)"
arbcall"void acb_hypgeom_legendre_q(acb_t res, const acb_t n, const acb_t m, const acb_t z, int type, slong prec)"
arbcall"void acb_hypgeom_legendre_p_uiui_rec(acb_t res, ulong n, ulong m, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_spherical_y(acb_t res, slong n, slong m, const acb_t theta, const acb_t phi, slong prec)"

### Dilogarithm
arbcall"void acb_hypgeom_dilog_zero_taylor(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_dilog_zero(acb_t res, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_dilog_transform(acb_t res, const acb_t z, int algorithm, slong prec)"
arbcall"void acb_hypgeom_dilog_continuation(acb_t res, const acb_t a, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_dilog_bitburst(acb_t res, acb_t z0, const acb_t z, slong prec)"
arbcall"void acb_hypgeom_dilog(acb_t res, const acb_t z, slong prec)"


text:
# Primary context management
#
# This is meant for interoperability with the CUDA runtime API

export
    CuPrimaryContext, unsafe_reset!, isactive, flags, setflags!

"""
    CuPrimaryContext(dev::CuDevice)

Create a primary CUDA context for a given device.

Each primary context is unique per device and is shared with CUDA runtime API. It is meant
for interoperability with (applications using) the runtime API.
"""
struct CuPrimaryContext
    dev::CuDevice

    function CuPrimaryContext(dev::CuDevice)
        pctx = new(dev)
        get!(pctx_instances, pctx, Set{WeakRef}())
        return pctx
    end
end

# keep a list of the contexts derived from a primary context.
# these need to be invalidated when we reset the primary context forcibly
# (as opposed to waiting for all derived contexts going out of scope).
const pctx_instances = Dict{CuPrimaryContext,Set{WeakRef}}()

"""
    CuContext(pctx::CuPrimaryContext)
    CuContext(f::Function, pctx::CuPrimaryContext)

Retain the primary context on the GPU, returning a context compatible with the driver API.
The primary context will be released when the returned driver context is finalized.

Contrary to the CUDA API, this call does push the context on the current CPU thread.
"""
function CuContext(pctx::CuPrimaryContext)
    handle = Ref{CuContext_t}()
    @apicall(:cuDevicePrimaryCtxRetain, (Ptr{CuContext_t}, CuDevice_t,), handle, pctx.dev)
    ctx = CuContext(handle[], false)    # CuContext shouldn't manage this ctx
    push!(CuContext, ctx)
    finalizer((ctx)->begin
        @assert isvalid(ctx)    # not owned by CuContext, so shouldn't have been invalidated
        @apicall(:cuDevicePrimaryCtxRelease, (CuDevice_t,), pctx.dev)
        invalidate!(ctx)
        delete!(pctx_instances[pctx], WeakRef(ctx))
    end, ctx)
    push!(pctx_instances[pctx], WeakRef(ctx))
    return ctx
end

function state(pctx::CuPrimaryContext)
    flags = Ref{Cuint}()
    active = Ref{Cint}()
    @apicall(:cuDevicePrimaryCtxGetState, (CuDevice_t, Ptr{Cuint}, Ptr{Cint}),
             pctx.dev, flags, active)
    return (flags[], active[] == one(Cint))
end

"""
    unsafe_reset!(pctx::CuPrimaryContext, [checked::Bool=true])

Explicitly destroys and cleans up all resources associated with a device's primary context
in the current process. Note that this forcibly invalidates all contexts derived from this
primary context, and as a result outstanding resources might become invalid.

It is normally unnecessary to call this function, as resource are automatically freed when
contexts go out of scope. In the case of primary contexts, they are collected when all
contexts derived from that primary context have gone out of scope.

The `checked` argument determines whether to verify that the primary context has become
inactive after resetting the derived driver contexts. This may not be possible, eg. if the
CUDA runtime API itself has retained an additional context instance.
"""
function unsafe_reset!(pctx::CuPrimaryContext, checked::Bool=true)
    for ref in pctx_instances[pctx]
        ctx = ref.value
        destroy!(ctx)
        finalize(ctx)
    end
    @assert isempty(pctx_instances[pctx])

    # NOTE: we don't support/perform the actual call to cuDevicePrimaryCtxReset, because of
    #       what's probably a bug in CUDA: calling cuDevicePrimaryCtxReset makes CUDA ignore
    #       all future calls to cuDevicePrimaryCtxRelease, even if those would be necessary
    #       to make the refcount of derived contexts instantiated through
    #       cuDevicePrimaryCtxRetain drop to 0. As a result, calling cuDevicePrimaryCtxReset
    #       makes that future activated primary contexts remain active indefinitely.
    #
    #       However, we don't _need_ cuDevicePrimaryCtxReset because we already forced
    #       finalization (and hence cuDevicePrimaryCtxRelease) on all derived contexts
    #       through the GC, and asserted that there's no derived contexts left.
    #@apicall(:cuDevicePrimaryCtxReset, (CuDevice_t,), pctx.dev)

    if checked
        # NOTE: having finalized all derived contexts doesn't mean the primary context is
        #       inactive now, because external consumers (eg. CUDArt) might hold another
        #       instance. That doesn't mean this logic isn't required, for correctness &
        #       ease of debugging wrt. data structures managed by CUDAdrv tied to contexts
        #       allocated via CUDAdrv. However, it requires us to still perform eg. a
        #       cudaDeviceReset call in CUDArt.
        @assert !isactive(pctx)
    end

    return
end

"""
    isactive(pctx::CuPrimaryContext)

Query whether a primary context is active.
"""
isactive(pctx::CuPrimaryContext) = state(pctx)[2]

"""
    flags(pctx::CuPrimaryContext)

Query the flags of a primary context.
"""
flags(pctx::CuPrimaryContext) = CUctx_flags(state(pctx)[1])

"""
    setflags!(pctx::CuPrimaryContext)

Set the flags of a primary context.
"""
setflags!(pctx::CuPrimaryContext, flags::CUctx_flags) =
    @apicall(:cuDevicePrimaryCtxSetFlags, (CuDevice_t, Cuint), pctx.dev, flags)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of highlighters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Highlighters" begin
    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\
    2 & true & 2.0 & 2 \\\\
    \\color{red}{3} & false & \\color{red}{3.0} & \\color{red}{3} \\\\
    4 & true & 4.0 & 4 \\\\
    5 & false & 5.0 & 5 \\\\
    6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    hl = LatexHighlighter((data, i, j)->data[i, j] == 3, "color{red}")
    result = pretty_table(
        String,
        data;
        backend = Val(:latex),
        highlighters = hl
    )

    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\
    2 & true & 2.0 & 2 \\\\
    \\textbf{\\color{red}{3}} & false & \\textbf{\\color{red}{3.0}} & \\textbf{\\color{red}{3}} \\\\
    4 & true & 4.0 & 4 \\\\
    5 & false & 5.0 & 5 \\\\
    6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    hl = LatexHighlighter((data, i, j)->data[i,j] == 3, ["textbf", "color{red}"])
    result = pretty_table(
        String,
        data;
        backend = Val(:latex),
        highlighters = hl
    )

    @test result == expected
end


text:
## Strings

const H5TYPE_VLEN_UTF8 = VariableLengthDatatype(DT_VARIABLE_LENGTH, 0x11, 0x01, 0x00,
                                                odr_sizeof(Vlen{UInt8}),
                                                FixedPointDatatype(1, false))

h5fieldtype(::JLDFile, ::Type{String}, ::Type{String}, ::Initialized) =
    H5TYPE_VLEN_UTF8
function h5fieldtype(f::JLDFile, writeas::Type{String},
                     readas::Type, init::Initialized)
    @lookup_committed f readas
    commit(f, h5fieldtype(f, writeas, writeas, init), writeas, readas)
end

fieldodr(::Type{String}, ::Bool) = Vlen{String}

# Stored as variable-length strings
odr_sizeof(x::FixedLengthString{String}) = x.length

h5type(f::JLDFile, writeas::Type{String}, x::String) =
    StringDatatype(typeof(x), jlsizeof(x))
h5type(f::JLDFile, writeas::Type{String}, x) =
    h5fieldtype(f, writeas, typeof(x), Val{true})
odr(::Type{String}) = fieldodr(String, true)
objodr(x::String) = FixedLengthString{String}(jlsizeof(x))

function jltype(f::JLDFile, dt::BasicDatatype)
    if dt.class == DT_STRING
        if (dt.bitfield1 == 0x01 || dt.bitfield1 == 0x11) && dt.bitfield2 == 0x00 && dt.bitfield3 == 0x00
            return FixedLengthString{String}(dt.size)
        else
            throw(UnsupportedFeatureException())
        end
    elseif dt.class == DT_OPAQUE
        error("attempted to read a bare (non-committed) opaque datatype")
    elseif dt.class == DT_REFERENCE
        return ReadRepresentation{Any,RelOffset}()
    else
        throw(UnsupportedFeatureException())
    end
end

function jltype(f::JLDFile, dt::VariableLengthDatatype)
    if dt == H5TYPE_VLEN_UTF8
        return ReadRepresentation{String,Vlen{String}}()
    else
        throw(UnsupportedFeatureException())
    end
end

function h5convert!(out::Pointers, fls::FixedLengthString, f::JLDFile, x, ::JLDWriteSession)
    fls.length == jlsizeof(x) || throw(InvalidDataException())
    (unsafe_copyto!(convert(Ptr{UInt8}, out), pointer(x), fls.length); nothing)
end
h5convert!(out::Pointers, ::Type{Vlen{String}}, f::JLDFile, x, wsession::JLDWriteSession) =
    store_vlen!(out, UInt8, f, unsafe_wrap(Vector{UInt8}, x), wsession)

jlconvert(::ReadRepresentation{String,Vlen{String}}, f::JLDFile, ptr::Ptr, ::RelOffset) =
    String(jlconvert(ReadRepresentation{UInt8,Vlen{UInt8}}(), f, ptr, NULL_REFERENCE))
function jlconvert(rr::FixedLengthString{String}, ::JLDFile, ptr::Ptr, ::RelOffset)
    data = Vector{UInt8}(undef, rr.length)
    unsafe_copyto!(pointer(data), convert(Ptr{UInt8}, ptr), rr.length)
    String(data)
end

# Used only for custom serialization
constructrr(::JLDFile, ::Type{String}, dt::VariableLengthDatatype{FixedPointDatatype}, ::Vector{ReadAttribute}) =
    dt == H5TYPE_VLEN_UTF8 ?
        (ReadRepresentation{String,Vlen{String}}(), true) :
        throw(UnsupportedFeatureException())

## Symbols

function h5fieldtype(f::JLDFile, ::Type{Symbol}, readas::Type, ::Initialized)
    @lookup_committed f readas
    commit(f, H5TYPE_VLEN_UTF8, Symbol, readas)
end
fieldodr(::Type{Symbol}, ::Bool) = Vlen{String}

h5type(f::JLDFile, ::Type{Symbol}, x) = h5fieldtype(f, Symbol, typeof(x), Val{true})
odr(::Type{Symbol}) = Vlen{String}

function h5convert!(out::Pointers, ::Type{Vlen{String}}, f::JLDFile, x::Symbol, ::JLDWriteSession)
    s = String(x)
    store_vlen!(out, UInt8, f, unsafe_wrap(Vector{UInt8}, s), f.datatype_wsession)
end

constructrr(::JLDFile, ::Type{Symbol}, dt::VariableLengthDatatype{FixedPointDatatype}, ::Vector{ReadAttribute}) =
    dt == H5TYPE_VLEN_UTF8 ? (ReadRepresentation{Symbol,Vlen{String}}(), true) :
                             throw(UnsupportedFeatureException())

jlconvert(::ReadRepresentation{Symbol,Vlen{String}}, f::JLDFile, ptr::Ptr, ::RelOffset) =
    Symbol(jlconvert(ReadRepresentation{UInt8,Vlen{UInt8}}(), f, ptr, NULL_REFERENCE))

## BigInts and BigFloats

writeas(::Union{Type{BigInt},Type{BigFloat}}) = String
wconvert(::Type{String}, x::BigInt) = string(x, base = 62)
wconvert(::Type{String}, x::BigFloat) = string(x)
rconvert(::Type{BigInt}, x::String) = parse(BigInt, x, base = 62)
rconvert(::Type{BigFloat}, x::String) = parse(BigFloat, x)

## Pointers

# Previously it was disallowed to serialize pointers.
# Due to popular demand and in particular to not error on serializing complex structures
# that contain non-essential pointers this has been changed to instead
# return null pointers.
writeas(::Type{Ptr{T}}) where {T} = Nothing
rconvert(::Type{Ptr{T}}, ::Nothing) where {T} = Ptr{T}(0)

## Arrays

h5fieldtype(::JLDFile, ::Type{T}, ::Type{T}, ::Initialized) where {T<:Array} =
    ReferenceDatatype()
fieldodr(::Type{T}, ::Bool) where {T<:Array} = RelOffset

@inline function odr(::Type{Array{T,N}}) where {T,N}
    writtenas = writeas(T)
    CustomSerialization(writtenas, T, fieldodr(writtenas, false))
end

# This is all so that when you define a writeas method to write something as an
# array, it writes a reference to the actual array where the datatype is
# committed and has a written_type attribute.
function h5fieldtype(f::JLDFile, ::Type{T}, readas::DataType,
                     ::Initialized) where T<:Array
    @lookup_committed f readas
    commit(f, ReferenceDatatype(), T, readas)
end
h5type(f::JLDFile, ::Type{T}, x) where {T<:Array} =
    h5fieldtype(f, T, typeof(x), Val{true})

function h5type(f::JLDFile, ::Type{T}, ::T) where T<:Array
    if T <: Array{Union{}}
        return ReferenceDatatype()
    end
    ty = T.parameters[1]
    writtenas = writeas(ty)
    if !hasfielddata(writtenas)
        # This is a hacky way to generate an instance of ty
        # the instance isn't actually needed for anything except that inside
        # h5type ty is determined via typeof(x)
        # annoyingly for some types h5type needs the instance
        h5type(f, writtenas, rconvert(ty, newstruct(writtenas)))
    else
        h5fieldtype(f, writtenas, ty, Val{false})
    end
end

_odr(writtenas::Type{T}, readas::Type{T}, odr) where {T<:Array} = odr
_odr(writtenas::Type{T}, readas::DataType, odr) where {T<:Array} =
    CustomSerialization{writtenas,RelOffset}

function constructrr(::JLDFile, ::Type{T}, dt::BasicDatatype,
                     attrs::Vector{ReadAttribute}) where T<:Array
    dt.class == DT_REFERENCE || throw(UnsupportedFeatureException())
    (ReadRepresentation{Array, RelOffset}(), true)
end

## SimpleVectors

writeas(::Type{Core.SimpleVector}) = Vector{Any}
wconvert(::Type{Vector{Any}}, x::Core.SimpleVector) = collect(Any, x)
rconvert(::Type{Core.SimpleVector}, x::Vector{Any}) = Core.svec(x...)

## Dicts

writeas(::Type{Dict{K,V}}) where {K,V} = Vector{Pair{K,V}}
writeas(::Type{IdDict}) = Vector{Pair{Any,Any}}
writeas(::Type{Base.ImmutableDict{K,V}}) where {K,V} = Vector{Pair{K,V}}
wconvert(::Type{Vector{Pair{K,V}}}, x::AbstractDict{K,V}) where {K,V} = collect(x)
function rconvert(::Type{T}, x::Vector{Pair{K,V}}) where {T<:AbstractDict,K,V}
    d = T()
    isa(d, Dict) && sizehint!(d::Dict, length(x))
    for (k,v) in x
        d[k] = v
    end
    d
end

function rconvert(::Type{<:Base.ImmutableDict}, x::Vector{Pair{K,V}}) where {K,V}
    @assert !isempty(x)
    d = Base.ImmutableDict(x[1])
    for p in (@view x[2:end])
        d = Base.ImmutableDict(d, p)
    end
    d
end

## NTuples
# Immutable objects are stored as HDF5 structs and inlined into
# parent structures. HDF5 only allows typemax(Unt16) bytes
# for struct description. NTuples are the most common offender for
# exploding struct size. (e.g. in the form of large StaticArrays)
# The definitions below prevent inlining of large NTuples and
# convert to an array instead.
const NTUPLE_INLINE_THRESHOLD = 10

function writeas(NT::Type{NTuple{N,T}}) where {N,T}
    if N > NTUPLE_INLINE_THRESHOLD
        return Vector{T}
    else
        return NT
    end
end

wconvert(::Type{Vector{T}}, x::NTuple{N,T}) where {N,T} = collect(x)
rconvert(::Type{NTuple{N,T}}, x::Vector{T}) where {N,T} = NTuple{N,T}(x)


text:
using Colors, Gadfly, RDatasets

set_default_plot_size(5inch,4inch)
iris = dataset("datasets","iris")

p = plot(iris, x=:SepalLength, y=:PetalLength, color=:Species, Geom.point,
     layer(Stat.smooth(method=:lm, levels=[0.90, 0.99]), Geom.line, Geom.ribbon),
    Theme(alphas=[0.6], key_position=:inside)
)


text:
module DaylightPlots

using RecipesBase
using Reexport, RecipesBase
@reexport using Base.Dates, TimeZones

include("sunrise.jl")
include("plots.jl")

end # module


text:
module _EnofEn
export EnofEn
using StatsBase: countmap, Histogram, fit
    """
        EoE, AvEn, S2 = EnofEn(Sig)

    Returns the entropy of entropy (`EoE`), the average Shannon entropy (`AvEn`), and
    the number of levels (`S2`) across all windows estimated from the data sequence (`Sig`) using
    the default parameters:
    window length (samples) = 10, slices = 10, logarithm = natural,
    heartbeat interval range (xmin, xmax) = (min(Sig), max(Sig))

        EoE, AvEn, S2 = EnofEn(Sig::AbstractArray{T,1} where T<:Real; tau::Int=10, S::Int=10, Xrange::Tuple{Real,REal}, Logx::Real=exp(1))

    Returns the entropy of entropy (`EoE`) estimated from the data sequence
    (`Sig`)  using the specified 'keyword' arguments:

    # Arguments:
    `tau`    - Window length, an integer > 1 \n
    `S`      - Number of slices (s1,s2), a two-element tuple of integers > 2 \n
    `Xrange` - The min and max heartbeat interval, a two-element tuple where X[1] <= X[2]\n
    `Logx`   - Logarithm base, a positive scalar  \n

    # See also `SampEn`, `MSEn`, `ApEn`

    # References:
        [1] Chang Francis Hsu, et al.,
            "Entropy of entropy: Measurement of dynamical complexity for
            biological systems."
            Entropy
            19.10 (2017): 550.

    """
    function EnofEn(Sig::AbstractArray{T,1} where T<:Real; tau::Int=10, S::Int=10,
                        Xrange::Tuple{Real,Real}=(minimum(Sig),maximum(Sig)), Logx::Real=exp(1))

    N = size(Sig,1)
    (N > 10) ? nothing :  error("Sig:   must be a numeric vector")
    (tau > 1 && tau < length(Sig)) ? nothing : error("tau:   must be an integer > 1")
    (S > 1) ? nothing : error("S:  must be an integer > 1")
    (length(Xrange)==2 && (Xrange[1]<=Xrange[2])) ? nothing : error("Xrange:   must be a two-element numeric tuple where Xrange[1]<Xrange[2]")
    (Logx>0) ? nothing :  error("Logx:     must be a positive number > 0")

    Wn = Int(floor(N/tau))
    Wj = transpose(reshape(Sig[1:Wn*tau],tau,Wn))
    Yj = zeros(Wn)
    #Edges = collect(range(minimum(Sig),maximum(Sig),length=(S[1]+1)))
    Edges = collect(range(Xrange[1],Xrange[2],length=(S+1)))
    Edges[1] -= .1;  Edges[end] += .1
    for n = 1:Wn
        Temp = fit(Histogram,Wj[n,:],Edges).weights/tau
        Temp = Temp[Temp.>0]
        Yj[n] = -sum(Temp.*log.(Logx, Temp))
    end

    AvEn = sum(Yj)/Wn
    #Edges = collect(range(minimum(Yj),maximum(Yj),length=(S[2]+1)))
    #Edges[1] -= .1;  Edges[end] += .1
    #Pjl = fit(Histogram,Yj,Edges).weights/Wn
    #Pjl = Pjl[Pjl.>0]
    Pjl = collect(values(countmap(round.(Yj,digits=12))))./Wn
    S2 = length(Pjl)
    if round(sum(Pjl),digits=5) != 1
        @warn("Possible error estimating probabilities")
    end

    EoE = -sum(Pjl.*log.(Logx, Pjl))
    return EoE, AvEn, S2

    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


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

"""
    @opaque ([type, ]args...) -> body

Marks a given closure as "opaque". Opaque closures capture the
world age of their creation (as opposed to their invocation).
This allows for more aggressive optimization of the capture
list, but trades off against the ability to inline opaque
closures at the call site, if their creation is not statically
visible.

An argument tuple type (`type`) may optionally be specified, to
specify allowed argument types in a more flexible way. In particular,
the argument type may be fixed length even if the function is variadic.

!!! warning
    This interface is experimental and subject to change or removal without notice.
"""
macro opaque(ex)
    esc(Expr(:opaque_closure, ex))
end

macro opaque(ty, ex)
    esc(Expr(:opaque_closure, ty, ex))
end

# OpaqueClosure construction from pre-inferred CodeInfo/IRCode
using Core.Compiler: IRCode
using Core: CodeInfo

function compute_ir_rettype(ir::IRCode)
    rt = Union{}
    for i = 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        if isa(stmt, Core.Compiler.ReturnNode) && isdefined(stmt, :val)
            rt = Core.Compiler.tmerge(Core.Compiler.argextype(stmt.val, ir), rt)
        end
    end
    return Core.Compiler.widenconst(rt)
end

function Core.OpaqueClosure(ir::IRCode, env...;
        nargs::Int = length(ir.argtypes)-1,
        isva::Bool = false,
        rt = compute_ir_rettype(ir))
    if (isva && nargs > length(ir.argtypes)) || (!isva && nargs != length(ir.argtypes)-1)
        throw(ArgumentError("invalid argument count"))
    end
    src = ccall(:jl_new_code_info_uninit, Ref{CodeInfo}, ())
    src.slotflags = UInt8[]
    src.slotnames = fill(:none, nargs+1)
    src.slottypes = copy(ir.argtypes)
    Core.Compiler.replace_code_newstyle!(src, ir, nargs+1)
    Core.Compiler.widen_all_consts!(src)
    src.inferred = true
    # NOTE: we need ir.argtypes[1] == typeof(env)

    ccall(:jl_new_opaque_closure_from_code_info, Any, (Any, Any, Any, Any, Any, Cint, Any, Cint, Cint, Any),
          Tuple{ir.argtypes[2:end]...}, Union{}, rt, @__MODULE__, src, 0, nothing, nargs, isva, env)
end

function Core.OpaqueClosure(src::CodeInfo, env...)
    M = src.parent.def
    sig = Base.tuple_type_tail(src.parent.specTypes)

    ccall(:jl_new_opaque_closure_from_code_info, Any, (Any, Any, Any, Any, Any, Cint, Any, Cint, Cint, Any),
          sig, Union{}, src.rettype, @__MODULE__, src, 0, nothing, M.nargs - 1, M.isva, env)
end

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
module VisUtilModule

using FinEtools
using PlotlyJS
using JSON
using LinearAlgebra: norm, I

const BACKGROUNDCOLOR = "rgb(230, 230,230)"

function removepair(nt, key)
    return (; delete!(Dict(pairs(nt)), key)...)
end

function addpair(nt, key, val)
    d = Dict()
    for p in pairs(nt)
        d[p[1]] = p[2]
    end
    d[key] = val
    return (; d...)
end

function local_frame!(F0, x0, x1x2_vector)
    # This is the element frame in the configuration t=0
    F0[:,1] = (x0[2,:]-x0[1,:]);
    L0 = norm(@view F0[:,1]);
    F0[:,1] /= L0;
    #     F0(:,3)=skewmat(F0(:,1))*x1x2_vector;
    F0[:,3] .= (-F0[3,1]*x1x2_vector[2]+F0[2,1]*x1x2_vector[3],
                 F0[3,1]*x1x2_vector[1]-F0[1,1]*x1x2_vector[3],
                -F0[2,1]*x1x2_vector[1]+F0[1,1]*x1x2_vector[2]);
    F0[:,3] /= norm(@view F0[:,3]);
    #     F0(:,2)=skewmat(F0(:,3))*F0(:,1);
    F0[:,2] .= (-F0[3,3]*F0[2,1]+F0[2,3]*F0[3,1],
                 F0[3,3]*F0[1,1]-F0[1,3]*F0[3,1],
                -F0[2,3]*F0[1,1]+F0[1,3]*F0[2,1]);
    return L0, F0
end

"""
    space_aspectratio(xyz; margin = 0.0)

Compute the aspect ratio of the space given by the array of locations
"""
function space_aspectratio(xyz; margin = 0.0)
    xr = maximum(xyz[:, 1]) - minimum(xyz[:, 1])
    yr = maximum(xyz[:, 2]) - minimum(xyz[:, 2])
    zr = maximum(xyz[:, 3]) - minimum(xyz[:, 3])
    if margin == 0.0
        margin = (xr + yr + zr) / 10000
    end
    xr += margin
    yr += margin
    zr += margin
    return [1.0, yr / xr, zr / xr]
end

"""
    plot_space_box(xyz; color = BACKGROUNDCOLOR)

Plot points to set up the global 3D box
"""
function plot_space_box(xyz; color = BACKGROUNDCOLOR)
    return [scatter3d(;x=xyz[:, 1],y=xyz[:, 2], z=xyz[:, 3], mode="markers",
        marker=attr(color=color, size=1, symbol="circle",
            line=attr(color=color, width=1)),
        line=attr(color=color, width=1))]
end

"""
    plot_points(xyz; color = "#1f77b4", msize = 4)

Plot points as markers
"""
function plot_points(xyz; color = "#1f77b4", msize = 4)
    return [scatter3d(;x=xyz[:, 1],y=xyz[:, 2], z=xyz[:, 3], mode="markers",
        marker=attr(color=color, size=msize, symbol="circle",
            line=attr(color=color, width=3)),
        line=attr(color="#1f77b4", width=4))]
end

"""
    plot_nodes(fens; color = "#1f77b4")

Plot nodes as markers
"""
function plot_nodes(fens; color = "#1f77b4")
    return plot_points(fens.xyz; color = color)
end

"""
    plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4)

Plot the midline
"""
function plot_midline(fens, fes; kwargs...)
    x = deepcopy(fens.xyz)
    if :x in keys(kwargs)
        x = kwargs[:x]; kwargs = removepair(kwargs, :x)
    end
    u = fill(0.0, size(x))
    if :u in keys(kwargs)
        u = kwargs[:u]; kwargs = removepair(kwargs, :u)
    end
    lwidth = 4
    if :lwidth in keys(kwargs)
        lwidth = kwargs[:lwidth]; kwargs = removepair(kwargs, :lwidth)
    end
    color = "rgb(155, 155, 255)"
    if :color in keys(kwargs)
        color = kwargs[:color]; kwargs = removepair(kwargs, :color)
    end
    xyz = fill(0.0, 2, 3)
    t = PlotlyBase.GenericTrace[]
    for c in fes.conn
        xyz[1, :] .= x[c[1], :] .+ u[c[1], :]
        xyz[2, :] .= x[c[2], :] .+ u[c[2], :]
        push!(t, scatter3d(;x=xyz[:, 1], y=xyz[:, 2], z=xyz[:, 3], mode="lines", line=attr(color=color, width=lwidth)))
    end
    return t
end

function _circular(dimensions, buffers, x1x2_vector, R1I, R1J; kwargs...)
    radius = dimensions[1]
    F0, x0, xt, c, s, xc, faces, facecolors = buffers
    L0, F0 = local_frame!(F0, x0, x1x2_vector)
    FtI=R1I*F0;
    FtJ=R1J*F0;
    nseg = length(c)-1
    for j in 1:nseg+1
        xc[j,:] = xt[1,:]+radius*c[j]*FtI[:,2]+radius*s[j]*FtI[:,3];
        xc[j+(nseg+1),:] = xt[2,:]+radius*c[j]*FtJ[:,2]+radius*s[j]*FtJ[:,3];
    end
    return mesh3d(;x=xc[:, 1],y=xc[:, 2], z=xc[:, 3], i=faces[:, 1].-1, j=faces[:, 2].-1, k=faces[:, 3].-1, facecolor=facecolors, kwargs...)
end

function _rectangular(dimensions, buffers, x1x2_vector, R1I, R1J; kwargs...)
    b, h = dimensions[1:2]./2 # departures from the midline by half the dimension
    F0, x0, xt, c, s, xc, faces, facecolors = buffers
    L0, F0 = local_frame!(F0, x0, x1x2_vector)
    FtI=R1I*F0;
    FtJ=R1J*F0;
    for j in 1:4+1
        xc[j,:] = xt[1,:]+c[j]*b*FtI[:,2]+s[j]*h*FtI[:,3];
        xc[j+(4+1),:] = xt[2,:]+c[j]*b*FtJ[:,2]+s[j]*h*FtJ[:,3];
    end
    return mesh3d(;x=xc[:, 1],y=xc[:, 2], z=xc[:, 3], i=faces[:, 1].-1, j=faces[:, 2].-1, k=faces[:, 3].-1, facecolor=facecolors, kwargs...)
end

"""
    plot_solid(fens, fes; kwargs...)

Plot beam structure as solid
"""
function plot_solid(fens, fes; kwargs...)
    x = deepcopy(fens.xyz)
    if :x in keys(kwargs)
        x = kwargs[:x]; kwargs = removepair(kwargs, :x)
    end
    u = fill(0.0, size(x))
    if :u in keys(kwargs)
        u = kwargs[:u]; kwargs = removepair(kwargs, :u)
    end
    R = fill(0.0, size(x, 1), 9)
    for j in 1:size(R, 1)
        R[j, :] .= Matrix(1.0 * I, 3, 3)[:]
    end
    if :R in keys(kwargs)
        R = kwargs[:R]; kwargs = removepair(kwargs, :R)
    end
    facecolor = "rgb(255, 155, 55)"
    if :facecolor in keys(kwargs)
        facecolor = kwargs[:facecolor]; kwargs = removepair(kwargs, :facecolor)
    end
    if !(:flatshading in keys(kwargs))
        kwargs = addpair(kwargs, :flatshading, true);
    end
    F0 = fill(0.0, 3, 3)
    x0 = fill(0.0, 2, 3)
    xt = fill(0.0, size(x0))
    R1I = Matrix(1.0 * I, 3, 3)
    R1J = Matrix(1.0 * I, 3, 3)
    if fes.crosssection.shape == "circle"
        nseg = 13
        v = vec((collect(1:nseg+1).-1)./nseg)
        c = cos.(2*pi*v);
        s = sin.(2*pi*v);
        faces = fill(0, 2*nseg, 3)
        for i in 1:nseg
            faces[i, :] .= (i+1, i, i+nseg+1)
            faces[i+nseg, :] .= (i+nseg+1, i+nseg+2, i+1)
        end
        xc = zeros(2*(nseg+1),3);
        facecolors = fill(facecolor, 2*nseg)
        buffers = (F0, x0, xt, c, s, xc, faces, facecolors)
        _draw = _circular
    elseif fes.crosssection.shape == "rectangle"
        faces = fill(0, 2*4, 3)
        for i in 1:4
            faces[i, :] .= (i+1, i, i+4+1)
            faces[i+4, :] .= (i+4+1, i+4+2, i+1)
        end
        c = [1, -1, -1, 1, 1]
        s = [1, 1, -1, -1, 1]
        xc = zeros(2*(4+1),3);
        facecolors = fill(facecolor, 2*4)
        buffers = (F0, x0, xt, c, s, xc, faces, facecolors)
        _draw = _rectangular
    end

    t = PlotlyBase.GenericTrace[]
    for i in 1:count(fes)
        c = fes.conn[i]
        x0[1, :] .= x[c[1], :]
        x0[2, :] .= x[c[2], :]
        xt[1, :] .= x[c[1], :] .+ u[c[1], :]
        xt[2, :] .= x[c[2], :] .+ u[c[2], :]
        R1I[:] .= R[c[1], :]
        R1J[:] .= R[c[2], :]
        push!(t, _draw(fes.dimensions[i], buffers, fes.x1x2_vector[i], R1I, R1J; kwargs...))
    end
    return t
end

"""
    render(traces; kwargs...)

Render a plot

Key word arguments:
- `config = PlotConfig()`
"""
function render(traces; kwargs...)
    layout = default_layout_3d(; kwargs...)
    if :layout in keys(kwargs)
        layout = kwargs[:layout]; kwargs = removepair(kwargs, :layout)
    end
    # Default options: show the "save to chart studio" button
    config = PlotConfig(
        showLink=true,
        plotlyServerURL="https://chart-studio.plotly.com"
        )
    # Should we override config because they were supplied as argument?
    if :config in keys(kwargs)
        config = kwargs[:config]; kwargs = removepair(kwargs, :config)
    end
    p = plot(traces, layout; config = config)
    display(p)
    return p
end

"""
    default_layout_3d(; kwargs...)

Set up the default layout for three-dimensional plots
"""
function default_layout_3d(; kwargs...)
    aspectratio = "data" # plotly.js tries to be smart by scaling back the scene to a cube when x/y/z are too dissimilar.
    # By setting scene.aspectratio: 'data', the x/y/z range are always honored.
    if :aspectratio in keys(kwargs)
        aspectratio = kwargs[:aspectratio]; kwargs = removepair(kwargs, :aspectratio)
    end
    title = ""
    if :title in keys(kwargs)
        title = kwargs[:title]; kwargs = removepair(kwargs, :title)
    end
    autosize = true
    if :autosize in keys(kwargs)
        autosize = kwargs[:autosize]; kwargs = removepair(kwargs, :autosize)
    end
    width = 360
    if :width in keys(kwargs)
        width = kwargs[:width]; kwargs = removepair(kwargs, :width)
    end
    height = 240
    if :height in keys(kwargs)
        height = kwargs[:height]; kwargs = removepair(kwargs, :height)
    end
    scene=attr(
        xaxis = attr(title="X", gridcolor="rgb(255, 255, 255)",
          zerolinecolor="rgb(255, 255, 255)",
          showbackground=true,
          backgroundcolor=BACKGROUNDCOLOR),
        yaxis = attr(title="Y", gridcolor="rgb(255, 255, 255)",
           zerolinecolor="rgb(255, 255, 255)",
           showbackground=true,
           backgroundcolor=BACKGROUNDCOLOR),
        zaxis = attr(title="Z", gridcolor="rgb(255, 255, 255)",
           zerolinecolor="rgb(255, 255, 255)",
           showbackground=true,
           backgroundcolor=BACKGROUNDCOLOR),
        camera = attr(
            up=attr(x=0, y=0, z=1),
            center=attr(x=0, y=0, z=0),
            eye=attr(x=1.25, y=1.25, z=1.25),
            projection = attr(type = "orthographic") # alternative: projection = attr(type = "perspective")
            ),
        aspectratio = aspectratio,
        aspectmode = "manual")
    if autosize
        layout = Layout(autosize=autosize, title=title, scene=scene, showlegend=false)
    else
        layout = Layout(width=width, height=height, title=title, scene=scene, showlegend=false,)
    end

    return layout
end

"""
    save_to_json(pl, fn)

Save plot to a JSON file
"""
function save_to_json(pl, fn)
    savejson(pl, fn)
end

"""
    plot_from_json(fn)

Retrieve plot from a JSON file and display it
"""
function plot_from_json(fn)
    pl = open(fn, "r") do f;
        plot(JSON.parse(Plot, String(read(f))))
    end
    return pl
end

end # FinEtoolsBeamsVis


text:
SetupCommBuffers(edgeNodes::Int, dom::AbstractDomain) = nothing
BuildMesh(nx::Int, edgeNodes::Int, edgeElems::Int, domNodes::Int, padded_domElems::Int, x_h::Vector{FT}, y_h::Vector{FT}, z_h::Vector{FT}, nodelist_h::Vector{Int}, dom::AbstractDomain) where {FT} = nothing

getMyRank(comm::MPI.Comm) = MPI.Comm_rank(comm)
getMyRank(::Nothing) = 0

getNumRanks(comm::MPI.Comm) = MPI.Comm_size(comm)
getNumRanks(::Nothing) = 1

getWtime(::MPI.Comm) = MPI.Wtime()
getWtime(::Nothing) = time()

comm_max(data::Float64, comm::MPI.Comm) = MPI.Allreduce(data, MPI.MAX, comm)
comm_max(data::Float64, ::Nothing) = data


text:
# Implements a Multilinear Expression data structure

export
   Monomial,
   MultilinearExpression

# DATA TYPE DEFINITIONS
"""
    Monomial

Specifies a monomial as a coefficient of numeric type T and a sorted list of variable indices.
We assume the index 0 represents the constant (variable-free) monomial.
"""
struct Monomial #{T}
    #coef::T # coefficient
    vars::Tuple{Vararg{Int}} # list of variable indices
    # function Monomial{T}(coef::T,vars) where T <: Number
    #     @assert Base.issorted(vars)
    #     new(coef,vars)
    # end
    function Monomial(vars)
        @assert Base.issorted(vars)
        new(vars)
    end
end
"""
    MultilinearExpression

Implements a multilinear expression data type as a dictionary of monoids to coefficients.

### Examples
```julia
# Creates a constant expression
@show c = MultilinearExpression(7.4)
# Creates expressions with single term
@show e1 = MultilinearExpression(4.0,1)
@show e2 = MultilinearExpression(3.0,1,3)
````
"""
struct MultilinearExpression
    terms::Dict{Monomial,<:Number}
    MultilinearExpression(d::Dict{Monomial,<:Number}) = new(d)
end
# CONSTRUCTORS
"Creates an expression with a single constant c"
MultilinearExpression(c::Number) = MultilinearExpression(Dict(Monomial(tuple()) => c))
#Monomial(coef::T, vars::Tuple{Vararg{Int}}) where T <: Number = Monomial{T}(coef,vars)
"Creates an expression with a single term (given as a pair monomial => coefficient"
MultilinearExpression(t::Pair{Monomial,<:Number}) = MultilinearExpression(Dict(t))
# function MultilinearExpression(e::MultilinearExpression,t::Pair{Monomial,<:Number})
#     terms = copy(e.terms)
#     terms[t[1]] = t[2]
#     MultilinearExpression(terms)
# end
"Creates an expression with given coefficient and variables (indices)."
MultilinearExpression(c::Number,v1::Int,v...) = MultilinearExpression(Dict(Monomial((v1,v...)) => c))
#monomial(coef,v1,v...) = MultilinearExpression( Monomial((v1,v...)) => coef )
#monomial(coef,v1,v...) = MultilinearExpression( (Monomial(coef,(v1,v...)),) )
#monomial(c::Number) = MultilinearExpression(Monomial(tuple()) => c)
"Returns the aditivity identity expression."
Base.zero(::Type{MultilinearExpression}) = MultilinearExpression(0) # falls back to most basic identity constant
Base.zero(e::MultilinearExpression) = MultilinearExpression(zero(first(e.terms)[2]))
"Returns the multiplicative identity expression."
Base.one(::Type{MultilinearExpression}) = MultilinearExpression(1) # falls back to most basic numeric identity constant
Base.one(e::MultilinearExpression) = MultilinearExpression(one(first(e.terms)[2]))
# # GETTERS
Base.getindex(e::MultilinearExpression,m::Monomial) = get(e.terms,m,zero(first(e.terms)[2]))
Base.setindex!(e::MultilinearExpression,c::Number,m::Monomial) = setindex!(e.terms,c,m)
@inline terms(e::MultilinearExpression) = e.terms
# coeffficient(m::Monomial) = m.coef
@inline variables(m::Monomial) = m.vars
# # HELPERS
# "Verify if monomials are ordered in an expression"
# issorted(e::MultilinearExpression) = length(e) < 2 || ((variables(e[1]) < variables(e[2])) && issorted(Base.tail(e)))
"Iterate over the pairs monomials => coeeficient in the expression."
Base.iterate(e::MultilinearExpression) = Base.iterate(terms(e))
Base.iterate(e::MultilinearExpression,state) = Base.iterate(terms(e),state)
"Returns the number of monomials in the expression."
Base.length(e::MultilinearExpression) = length(terms(e))
# # INPUT/OUTPUT
Base.show(io::IO, m::Monomial) = print(io, "χ", join(variables(m),"χ"))
function Base.show(io::IO, t::Pair{Monomial,<:Number})
    if length(variables(t[1])) == 0
        print(io, t[2])
    else
        print(io, t[2], t[1])
    end
    nothing
end
#Base.show(io::IO, e::MultilinearExpression) = print(io, join(map(string,collect(terms(e)))," + "))
Base.show(io::IO, e::MultilinearExpression) = print(io, join(terms(e), " + "))
# # OPERATIONS
# """
#     Adds two expressions. Assumes the terms are ordered in each expression.
# """
function Base.:+(e1::MultilinearExpression,e2::MultilinearExpression)
    e = MultilinearExpression(copy(terms(e1)))
    for (m,c) in e2
        e[m] += c
    end
    e
#     @assert issorted(e1) && issorted(e2)
#     i, j = 1, 1
#     e = Vector{Monomial}()
#     while i <= length(e1) && j <= length(e2)
#         c1, c2 = coeffficient(e1[i]), coeffficient(e2[j])
#         v1, v2 = variables(e1[i]), variables(e2[j])
#         if v1 == v2
#             if c1 + c2 != zero(c1)
#                 push!(e,Monomial(c1 + c2, v1))
#             end
#             i += 1
#             j += 1
#         elseif v1 < v2
#             if c1 != zero(c1)
#                 push!(e,e1[i])
#             end
#             i += 1
#         else
#             if c2 != zero(c2)
#                 push!(e,e2[j])
#             end
#             j += 1
#         end
#     end
#     for k = i:length(e1)
#         push!(e,e1[k])
#     end
#     for k = j:length(e2)
#         push!(e,e2[k])
#     end
#     Tuple(e)
end
## Base.:+(e::MultilinearExpression,m::Monomial) = Base.:+(e,(m,))
# # Base.:+(m::Monomial,e::MultilinearExpression) = Base.:+((m,),e)
Base.:+(c::Number,e::MultilinearExpression) = MultilinearExpression(c) + e
Base.:+(e::MultilinearExpression,c::Number) = c + e
# Base.:*(c::T,m::Monomial{T}) where T <: Number = Monomial(c*coeffficient(m),variables(m))
Base.:*(c::Number,m::Monomial) = MultilinearExpression(m => c)
Base.:*(m::Monomial,c::Number) = c * m
function Base.:*(m1::Monomial,m2::Monomial)
    vars = union(variables(m1),variables(m2))
    @assert length(vars) == length(variables(m1)) + length(variables(m2)) "Multiplication of monomial with common variable is not multilinear"
    Monomial(Tuple(sort!(vars)))
end
function Base.:*(e1::MultilinearExpression,e2::MultilinearExpression)
    MultilinearExpression(Dict(
        m1 * m2 => c1 * c2
        for (m1,c1) in e1, (m2,c2) in e2
    ))
    # for (m1,c1) in e1
    #     for (m2,c2) in e2
    #         e[ m1 * m2 ] = c1 * c2
    #     end
    # end
    # e
    # @assert issorted(e1) && issorted(e2)
    # e = Dict{Tuple{Vararg{Int}},Monomial}()
    # Tuple( map(x -> x[2], sort!( collect(e), by = x -> x[1] ) ) )
end
Base.:*(c::Number,e::MultilinearExpression) = MultilinearExpression(c) * e
Base.:*(e::MultilinearExpression,c::Number) = c * e
## COMPARISONS
Base.:(==)(m1::Monomial,m2::Monomial) = variables(m1) == variables(m2)
function Base.:(==)(e1::MultilinearExpression,e2::MultilinearExpression)
    if length(e1) != length(e2)
        return false
    end
    for (m1,c1) in e1
        if e2[m1] != c1
            return false
        end
    end
    true
end
Base.hash(m::Monomial, h::UInt) = hash(variables(m), h)
Base.hash(t::Pair{Monomial,<:Number}, h::UInt) = hash(t[1],hash(t[2],h))
Base.hash(e::MultilinearExpression, h::UInt) = foldl( (t1,t2) -> hash(t1,hash(t2,h)), e ; init=hash(0))
## PROMOTION AND CONVERSION
Base.convert(::Type{MultilinearExpression},c::Number) = MultilinearExpression(c)
Base.convert(::Type{MultilinearExpression},p::Pair{Monomial,Number}) = MultilinearExpression(p)
Base.promote_rule(::Type{MultilinearExpression}, ::Type{T}) where T <: Number = MultilinearExpression


text:
function tm_score(P, Q, matches; metric=SqEuclidean())
    matchedP, matchedQ = matched_points(P, Q, matches)
    LN = size(Q,2)
    tform = kabsch(matchedP, matchedQ)
    return 1/LN * sum(1 ./ (1 .+ colwise(metric, tform(matchedP), matchedQ) ./ d0(LN)^2))
end