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

"""

    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.
* `upscale=1`: upscale 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, upscale=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, upscale=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, upscale=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
    if alpha!=1 && mode==0; error("Gradient of pool(alpha!=1,mode=0) broken in CUDNN"); end
    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...)

"""

Unpooling; `reverse` of pooling.

    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, upscale=1, mode=0)
        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(upscale,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(upscale,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(upscale,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)

function cdims(w,x; padding=0, stride=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)
            1 + div(size(x,i) - size(w,i) + 2*pi, si)
        elseif i == N-1
            size(w,N)
        else # i == N
            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 dcdims(w,x; padding=0, stride=1, o...)
    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)
            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 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)

### CPU convolution using im2col from Mocha.jl

# w=Ww,Hw,Cx,Cy
# x=Wx,Hx,Cx,Nx
# y=Wy,Hy,Cy,Nx
# if we apply im2col to a single image:
# w2=(Ww*Hw*Cx),Cy  ;; simple reshape
# x2=(Wy*Hy),(Ww*Hw*Cx)
# y2=(Wy*Hy),Cy     ;; simple reshape after y2=x2*w2

function conv4(w::AbstractArray{T,4}, x::AbstractArray{T,4};
                  padding=0, stride=1, upscale=1, mode=0, alpha=1,
                  o...) where {T} # Ignoring handle, algo, workSpace, workSpaceSizeInBytes
    if upscale != 1; throw(ArgumentError("CPU conv4 only supports upscale=1.")); end
    if mode != 0 && mode != 1; throw(ArgumentError("conv4 only supports mode=0 or 1.")); end
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    if Cx!=C1; throw(DimensionMismatch()); end
    Wy,Hy,Cy,Ny = cdims(w,x;padding=padding,stride=stride)
    # @assert Cy==C2 && Ny==Nx
    y = similar(x, (Wy,Hy,Cy,Ny))
    x2dims = im2col_dims(w,x,y)
    x2 = similar(x, x2dims)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    M,N,K,Y = Wy*Hy,Cy,Ww*Hw*Cx,Wy*Hy*Cy
    alpha,beta,yidx = T(alpha),T(0),1
    @inbounds for n in 1:Nx
        im2col!(w, x, x2, n, p1, p2, s1, s2, mode)
        gemm!('N','N',M,N,K,alpha,pointer(x2),pointer(w),beta,pointer(y,yidx))
        yidx += Y
    end
    return y
end

function conv4w(w::AbstractArray{T,4},x::AbstractArray{T,4},dy::AbstractArray{T,4};
                   padding=0, stride=1, upscale=1, mode=0, alpha=1,
                   o...) where {T} # Ignoring handle, algo, workSpace, workSpaceSizeInBytes
    # dw = x'*dy
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    Wy,Hy,Cy,Ny = size(dy)
    # if upscale != 1; throw(ArgumentError("CPU conv4 only supports upscale=1.")); end
    # if mode != 0 && mode != 1; throw(ArgumentError("conv4 only supports mode=0 or 1.")); end
    # @assert Cx==C1 && Cy==C2 && Ny==Nx
    dw = zero(w)
    x2dims = im2col_dims(w,x,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Cy,Ww*Hw*Cx,Cy,Wy*Hy
    alpha,beta = T(alpha),T(1)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        im2col!(w, x, x2, n, p1, p2, s1, s2, mode)
        gemm!('T','N',M,N,K,alpha,pointer(x2),pointer(dy,dyi),beta,pointer(dw))
        dyi += Y
    end
    return dw
end

function conv4x(w::AbstractArray{T,4},x::AbstractArray{T,4},dy::AbstractArray{T,4};
                   padding=0, stride=1, upscale=1, mode=0, alpha=1,
                   o...) where {T} # Ignoring handle, algo, workSpace, workSpaceSizeInBytes
    # dx = dy*w'
    Wx,Hx,Cx,Nx = size(x)
    Ww,Hw,C1,C2 = size(w)
    Wy,Hy,Cy,Ny = size(dy)
    # if upscale != 1; throw(ArgumentError("CPU conv4 only supports upscale=1.")); end
    # if mode != 0 && mode != 1; throw(ArgumentError("conv4 only supports mode=0 or 1.")); end
    @assert Cx==C1 && Cy==C2 && Ny==Nx
    dx = similar(x)
    x2dims = im2col_dims(w,x,dy)
    x2 = similar(x, x2dims)
    # op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
    Y,M,N,K = Wy*Hy*Cy,Wy*Hy,Ww*Hw*Cx,Cy
    alpha,beta = T(alpha),T(0)
    (p1,p2) = psize(padding,x)
    (s1,s2) = psize(stride,x)
    dyi = 1
    @inbounds for n in 1:Nx
        gemm!('N','T',M,N,K,alpha,pointer(dy,dyi),pointer(w),beta,pointer(x2))
        col2im!(w,dx,x2,n,p1,p2,s1,s2,mode)
        dyi += Y
    end
    return dx
end

im2col_dims(w,x,y)=(size(y,1)*size(y,2), size(w,1)*size(w,2)*size(w,3))

# Functions from conv.cpp:

for (T,S) in ((Float32,32), (Float64,64)); @eval begin

    function im2col!(w::AbstractArray{$T,4}, x::AbstractArray{$T,4}, x2::AbstractArray{$T,2},
                     n::Int, p1::Int, p2::Int, s1::Int, s2::Int, mode::Int)
        Wx,Hx,Cx,Nx = size(x)
        Ww,Hw,C1,C2 = size(w)
        xn = pointer(x, Wx*Hx*Cx*(n-1)+1)
        @knet8($("im2col$S"),
               (Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
               xn,x2,Wx,Hx,Cx,Ww,Hw,p1,p2,s1,s2,mode)
        return x2
    end

    function col2im!(w::AbstractArray{$T,4}, x::AbstractArray{$T,4}, x2::AbstractArray{$T,2},
                     n::Int, p1::Int, p2::Int, s1::Int, s2::Int, mode::Int)
        Wx,Hx,Cx,Nx = size(x)
        Ww,Hw,C1,C2 = size(w)
        xn = pointer(x, Wx*Hx*Cx*(n-1)+1)
        @knet8($("col2im$S"),
               (Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
               x2,xn,Wx,Hx,Cx,Ww,Hw,p1,p2,s1,s2,mode)
        return x
    end

    ### CPU pooling from Mocha.jl

    function pool(x::AbstractArray{$T,4}; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0, alpha=1, handle=nothing)
        if maxpoolingNanOpt!=0; throw(ArgumentError("CPU pool only supports maxpoolingNanOpt=0")); end
        Wx,Hx,Cx,Nx = size(x);
        Wy,Hy,Cy,Ny = pdims(x;window=window,padding=padding,stride=stride)
        y = similar(x, (Wy,Hy,Cy,Ny))
        (w1,w2) = psize(window, x)
        (p1,p2) = psize(padding, x)
        (s1,s2) = psize(stride, x)
        if mode == 0
            @knet8($("max_pooling_fwd$S"),
                   (Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                   x,y,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
        elseif mode == 1 || (mode == 2 && p1==p2==0)
            @knet8($("mean_pooling_fwd$S"),
                   (Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                   x,y,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
        else
            throw(ArgumentError("mode $mode not supported by cpu pool"))
        end
        if alpha != 1; lmul!(alpha,y); end
        return y
    end

    function poolx(x::AbstractArray{$T,4}, y::AbstractArray{$T,4}, dy::AbstractArray{$T,4};
                   window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0, alpha=1, handle=nothing)
        if maxpoolingNanOpt!=0; throw(ArgumentError("CPU pool only supports maxpoolingNanOpt=0")); end
        Wx,Hx,Cx,Nx = size(x);
        Wy,Hy,Cy,Ny = size(y);
        dx = similar(x)
        (w1,w2) = psize(window, x)
        (p1,p2) = psize(padding, x)
        (s1,s2) = psize(stride, x)
        if mode == 0
            if alpha != 1; y = y ./ alpha; end
            @knet8($("max_pooling_bwd$S"),
                   (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                   x,y,dy,dx,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
        elseif mode == 1 || (mode == 2 && p1==p2==0)
            @knet8($("mean_pooling_bwd$S"),
                   (Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                   dx,dy,Wx,Hx,Cx,Nx,Wy,Hy,w1,w2,p1,p2,s1,s2)
        else
            throw(ArgumentError("mode $mode not supported by cpu pool"))
        end
        if alpha != 1; lmul!(alpha,dx); end
        return dx
    end
end;end

# 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('*')
        @cudnn(cudnnFindConvolutionForwardAlgorithm,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr),
              handle,TD(x),FD(w),CD(w,x;o...),TD(y),requestedAlgoCount,returnedAlgoCount,perfResults)
        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('*')
        @cudnn(cudnnFindConvolutionBackwardFilterAlgorithm,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr),
              handle,TD(x),TD(dy),CD(w,x;o...),FD(dw),requestedAlgoCount,returnedAlgoCount,perfResults)
        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('*')
        @cudnn(cudnnFindConvolutionBackwardDataAlgorithm,
              (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Cint},Cptr),
              handle,FD(w),TD(dy),CD(w,x;o...),TD(dx),requestedAlgoCount,returnedAlgoCount,perfResults)
        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 = div(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

# TODO: this assumes one workspace per gpu. What about streams?
const CUDNN_WORKSPACE = []
function cudnnWorkSpace(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


text:
"""
The dimension of the space of Lagrange shape functions of degree d over a simplex of
dimension n is binom(n+d,d) == binom(n+d,n)
"""
function lagdimension end

# D: degree
# C: continuity
# M: mesh type
# T: field type
# NF: number of local shape functions
mutable struct LagrangeBasis{D,C,M,T,NF,P} <: Space{T}
  geo::M
  fns::Vector{Vector{Shape{T}}}
  pos::Vector{P}
end

# Constructor that automatically deduces MeshType and ScalarType but requires specification
# of the Degree, Cont, and NumFns type parameters
function LagrangeBasis{D,C,N}(mesh::M, fns::Vector{Vector{Shape{T}}}, pos::Vector{P}) where {P,T,M,N,C,D}
    LagrangeBasis{D,C,M,T,N,P}(mesh, fns, pos)
end

refspace(space::LagrangeBasis{D,C,M,T,NF}) where {D,C,M,T,NF} = LagrangeRefSpace{T,D,dimension(geometry(space))+1,NF}()
subset(s::S,I) where {S <: Space} = S(s.geo, s.fns[I], s.pos[I])

function lagrangecxd0(mesh)

    U = universedimension(mesh)
    D1 = dimension(mesh)+1

    geometry = mesh
    num_cells = numcells(mesh)

    # create the local shapes
    fns = Vector{Vector{Shape{Float64}}}(undef,num_cells)
    pos = Vector{vertextype(mesh)}(undef,num_cells)
    for (i,cell) in enumerate(cells(mesh))
        fns[i] = [Shape(i, 1, 1.0)]
        pos[i] = cartesian(center(chart(mesh, cell)))
  end

  NF = 1
  LagrangeBasis{0,-1,NF}(geometry, fns, pos)
end

"""
    lagrangec0d1(mesh[, bnd])

Construct the basis of continuous, piecewise linear basis functions subordinate to mesh `mesh`. Basis functions are constructed at vertices in the interionr of the mesh and on the closure of 'bnd'. In particular, leaving out the second argument creates a finite element space subject to homogeneous Dirichlet boundary conditions.
"""
function lagrangec0d1_dirichlet(mesh)

    verts = skeleton(mesh, 0)
    detached = trues(numvertices(mesh))
    for v in cells(verts)
        detached[v] = false
    end

    bnd = boundary(mesh)
    bndverts = skeleton(bnd, 0)
    notonbnd = trues(numvertices(mesh))
    for v in cells(bndverts)
        notonbnd[v] = false
    end

    vertexlist = findall(notonbnd .& .!detached)
    lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1})
end

function interior_and_junction_vertices(mesh, jct)
    verts = skeleton(mesh, 0)
    detached = trues(numvertices(mesh))
    for v in cells(verts)
        detached[v] = false
    end

    bndfaces = boundary(mesh)
    bndverts = skeleton(bndfaces, 0)
    notonbnd = trues(numvertices(mesh))
    for v in cells(bndverts)
        notonbnd[v] = false
    end

    onjct = broadcast(!, notonbnd)
    overlap_with_junction = overlap_gpredicate(jct)
    for indices in cells(bndfaces)
        bndface = simplex(vertices(bndfaces, indices))
        if overlap_with_junction(bndface)
            continue
        else
            for v in indices
                onjct[v] = false
            end
        end
    end

    vertexlist = findall(broadcast(|, onjct, notonbnd) .& broadcast(!,detached))
end

"""
    duallagrangecxd0(mesh, jct) -> basis

Build dual Lagrange piecewise constant elements. Boundary nodes are only considered if they are in the interior of `jct`.
"""
function duallagrangecxd0(mesh, jct=CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1))
    vertexlist = interior_and_junction_vertices(mesh, jct)
    duallagrangecxd0(mesh, vertexlist)
end

function duallagrangecxd0(mesh, vertexlist::Vector{Int})

    T = coordtype(mesh)

    fns = Vector{Vector{Shape{T}}}(undef,length(vertexlist))
    pos = Vector{vertextype(mesh)}()

    fine = barycentric_refinement(mesh)
    vtoc, vton = vertextocellmap(fine)
    verts = vertices(mesh)
    for (k,v) in enumerate(vertexlist)
        n = vton[v]
        F = vtoc[v,1:n]
        fns[k] = singleduallagd0(fine, F, v)
        push!(pos, verts[v])
    end

    NF = 1
    LagrangeBasis{0,-1,NF}(fine, fns, pos)
end

function duallagrangecxd0(mesh, vertices::CompScienceMeshes.AbstractMesh{U,1}) where {U}
    vertexlist = Int[v[1] for v in vertices]
    return duallagrangecxd0(mesh, vertexlist)
end

"""
    singleduallagd0(fine, F, v)

Build a single dual constant Lagrange element a mesh `fine`. `F` contains the indices to cells in the support and v is the index in the vertex list of the defining vertex.
"""
function singleduallagd0(fine, F, v)

    T = coordtype(fine)
    fn = Shape{T}[]
    for cellid in F
        cell = cells(fine)[cellid]
        ptch = chart(fine, cell)
        coeff = 1 / volume(ptch) / length(F)
        refid = 1
        push!(fn, Shape(cellid, refid, coeff))
    end

    return fn
end

"""
    lagrangec0d1(mesh; dirichlet=[true|false]) -> basis

Build lagrangec0d1 elements, including (dirichlet=false) or excluding (dirichlet=true) those attached to boundary vertices.
"""
function lagrangec0d1(mesh; dirichlet::Bool=true)
    if dirichlet == false
        return lagrangec0d1(mesh, boundary(mesh))
    else
        return lagrangec0d1_dirichlet(mesh)
    end
end

function lagrangec0d1(mesh, jct)
    vertexlist = interior_and_junction_vertices(mesh, jct)
    lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1})
end

# build continuous linear Lagrange elements on a 2D manifold
function lagrangec0d1(mesh, vertexlist::Vector, ::Type{Val{3}})

    T = coordtype(mesh)
    U = universedimension(mesh)

    cellids, ncells = vertextocellmap(mesh)

    Cells = cells(mesh)
    Verts = vertices(mesh)

    # create the local shapes
    fns = Vector{Shape{Float64}}[]
    pos = Vector{vertextype(mesh)}()

    sizehint!(fns, length(vertexlist))
    sizehint!(pos, length(vertexlist))
    for v in vertexlist

        numshapes = ncells[v]
        numshapes == 0 && continue

        shapes = Vector{Shape{Float64}}(undef,numshapes)
        for s in 1: numshapes
            c = cellids[v,s]
            # cell = mesh.faces[c]
            cell = Cells[c]

            localid = something(findfirst(isequal(v), cell),0)
            @assert localid != 0

            shapes[s] = Shape(c, localid, 1.0)
        end

        push!(fns, shapes)
        push!(pos, Verts[v])
    end

    NF = 3
    LagrangeBasis{1,0,NF}(mesh, fns, pos)
end

# for manifolds of dimension 1
function lagrangec0d1(mesh, vertexlist, ::Type{Val{2}})

    T = coordtype(mesh)
    U = universedimension(mesh)
    P = vertextype(mesh)

    geometry = mesh

    cellids, ncells = vertextocellmap(mesh)

    # create the local shapes
    numverts = numvertices(mesh)

    fns = Vector{Vector{Shape{T}}}()
    pos = Vector{P}()

    sizehint!(fns, length(vertexlist))
    sizehint!(pos, length(vertexlist))
    for v in vertexlist

        numshapes = ncells[v]
        numshapes == 0 && continue # skip detached vertices

        shapes = Vector{Shape{Float64}}(undef,numshapes)
        for s in 1: numshapes
            c = cellids[v,s]
            cell = mesh.faces[c]
            if cell[1] == v
                shapes[s] = Shape(c, 1, 1.0)
            elseif cell[2] == v
                shapes[s] = Shape(c, 2, 1.0)
            else
                error("Junctions not supported")
            end
        end

        push!(fns, shapes)
        push!(pos, mesh.vertices[v])
    end

    NF = 2
    LagrangeBasis{1,0,NF}(geometry, fns, pos)
end

function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where {U})

    Conn = connectivity(nodes, mesh, abs)
    rows = rowvals(Conn)
    vals = nonzeros(Conn)

    T = coordtype(mesh)
    P = vertextype(mesh)
    S = Shape{T}

    fns = Vector{Vector{S}}()
    pos = Vector{P}()
    for (i,node) in enumerate(nodes)
        fn = Vector{S}()
        for k in nzrange(Conn,i)
            cellid = rows[k]
            refid  = vals[k]
            push!(fn, Shape(cellid, refid, 1.0))
        end
        push!(fns,fn)
        push!(pos,cartesian(center(chart(nodes,node))))
    end

    LagrangeBasis{1,0,3}(mesh, fns, pos)
end

duallagrangec0d1(mesh) = duallagrangec0d1(mesh, barycentric_refinement(mesh), x->false, Val{dimension(mesh)+1})
function duallagrangec0d1(mesh, jct)
    jct_pred = inclosure_gpredicate(jct)
    refined = barycentric_refinement(mesh)
    duallagrangec0d1(mesh, refined, jct_pred, Val{dimension(mesh)+1})
end

function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})

    T = coordtype(mesh)
    num_faces = dimension(mesh)+1

    fns = Vector{Vector{Shape{T}}}(undef,numcells(mesh))
    pos = Vector{vertextype(mesh)}()

    # store the fine mesh's vertices in an octree for fast retrieval
    fine_vertices = Octree(vertices(refined))
    uv_ctr = ones(dimension(mesh))/(dimension(mesh)+1)

    vtoc, vton = vertextocellmap(refined)
    for (i,coarse_idcs) in enumerate(cells(mesh))
        fns[i] = Vector{Shape{T}}()
        push!(pos, cartesian(center(chart(mesh, coarse_idcs))))

        # It is assumed the vertices of this cell have the same index
        # mesh and its refinement.
        coarse_cell = chart(mesh, coarse_idcs)

        # get the index in fine.vertices of the centroid of coarse_cell
        centroid = barytocart(coarse_cell, uv_ctr)
        I = CollisionDetection.find(fine_vertices, centroid)
        @assert length(I) == 1
        centroid_id = I[1]

        # get the indx in fine.vertices of the centroid of the faces of coarse_cell
        face_center_ids = Vector{Int}(undef,num_faces)
        for f in 1:num_faces

            # prepare the barycentric coordinate
            uv_face_ctr = ones(dimension(mesh)+1)/(dimension(mesh))
            uv_face_ctr[f] = 0
            uv_face_ctr = uv_face_ctr[1:end-1]

            face_ctr = barytocart(coarse_cell, uv_face_ctr)
            I = CollisionDetection.find(fine_vertices, face_ctr)
            @assert length(I) == 1
            face_center_ids[f] = I[1]
        end

        n = vton[centroid_id]
        for c in vtoc[centroid_id,1:n]
            fine_idcs = cells(refined)[c]
            local_id = something(findfirst(isequal(centroid_id), fine_idcs), 0)
            @assert local_id != 0
            shape = Shape(c, local_id, 1.0)
            push!(fns[i], shape)
        end

        for f in 1:num_faces
            v = face_center_ids[f]
            jct_pred(vertices(refined)[v]) && continue
            n = vton[v]
            for c in vtoc[v,1:n]
                fine_idcs = cells(refined)[c]
                local_id = something(findfirst(isequal(v), fine_idcs),0)
                @assert local_id != 0
                shape = Shape(c, local_id, 1/n/2)
                push!(fns[i], shape)
            end
        end

        for f in 1:length(coarse_idcs)
            v = coarse_idcs[f]
            jct_pred(vertices(refined)[v]) && continue
            n = vton[v]
            for c in vtoc[v,1:n]
                fine_idcs = cells(refined)[c]
                local_id = something(findfirst(isequal(v), fine_idcs),0)
                @assert local_id != 0
                shape = Shape(c, local_id, 1/n/2)
                push!(fns[i], shape)
            end
        end
    end

    @assert length(fns) == length(pos)

    NF = 3
    return LagrangeBasis{1,0,NF}(refined, fns, pos)
end

"""
    duallagrangec0d1(originalmesh, refinedmesh)

It is the user responsibility to provide two meshes representing the same object.
The second mesh needs to be obtained using "barycentric_refinement(originalmesh)".
This basis function creats the dual Lagrange basis function and return an object that contains array of shapes [fns]
It also return a gemoetry containing the refined mesh.
"""
function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
  T = Float64
  U = universedimension(mesh)
  # get the information about number of vertices, number of faces , and the maping between vertices and faces for the original mesh
  numverts1 = numvertices(mesh)
  num_cells1 = numcells(mesh)
  cellids1, ncells1=vertextocellmap(mesh)
  # obtain the refined mesh from the original mesh
  #mesh2 = barycentric_refinement(mesh)
  # get the information about number of vertices, number of faces , and the maping between vertices and faces for the refined mesh
  num_cells2 = numcells(mesh2)
  numverts2 = numvertices(mesh2)
  geometry = mesh2
  cellids2, ncells2 = vertextocellmap(mesh2)

  fns = Vector{Vector{Shape{T}}}(undef,num_cells1)
  pos = Vector{vertextype(mesh)}()
  # We will iterate over the coarse mesh segments to assign all the functions to it.
  for segment_coarse in 1 : num_cells1
    # For the dual Lagrange there is a 6 shapes per segment
    numshapes = (ncells1[segment_coarse]*4) -2
    shapes = Vector{Shape{T}}(undef,numshapes)
    # Now we will get all the smaller faces within the coarse segment
    #i.e The coose segment will have two points, and these tow points are connected to two segmesnts in the finer mesh
    # This will give us a 4 smaller faces per Dual lagrange basis, we store them first in all_faces
    all_faces= Array{SVector{2,Int}}(undef,4)                      # faces in the original segment (4)
    # the follwoing code get the verteciec for each coarse segment
    # then it looks for the two faces connected to each point in the finer mesh
    # if for example the segment would connect to more than two faces we will have
    # mesh.faces[segment_coarse][1],n] and iterate over how many segments [n] are connected
    all_faces[1]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],1]]
    all_faces[2]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],2]]
    all_faces[3]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],1]]
    all_faces[4]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],2]]

    # now we now the first point of the corse segment will have the left hand side basis
    # and we know that that point is connected to faces 1,2 in the array all_faces
    # but now we want to which one of them is inner segment and which one is at the edge
    # as both of these two faces are in left hand side
    # The inner face will have the left hand side corase point in its first place for example
    # we have original point 2,3 , and the coarse segment (2,3)
    # after refinment we will have for example  points 20,2,21,3,22 with four faces (20,2),(2,21),(21,3),(3,22)
    # so the faces to the left (20,2),(2,21) we can decide the inner one if the original point(2) is at first index
    # same for right hand side faces (21,3),(3,22) if the original point(3) is in the second index
    # then the inner faces are (2,21) and (21,3)
    # the following code does the same checking process and assign the shapes for the dual lagragne right away

    # For the Left hand side faces
    if(all_faces[1][1] == mesh.faces[segment_coarse][1])
        shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,0.5)
        shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],1],2,0.5)
        shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,1.0)
    elseif(all_faces[2][1] == mesh.faces[segment_coarse][1])
        shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,1.0)
        shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],2],2,0.5)
        shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,0.5)
    end
    # For the Right hand side faces
    if(all_faces[3][2] == mesh.faces[segment_coarse][2])
        shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,1.0)
        shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],1],1,0.5)
        shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,0.5)
    elseif(all_faces[4][2] == mesh.faces[segment_coarse][2])
        shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,1.0)
        shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],2],1,0.5)
        shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,0.5)
    end
    # Now assign all of these shapes to the relevent segment in the coarse mesh
    fns[segment_coarse]=shapes
    push!(pos, cartesian(center(chart(mesh, mesh.faces[segment_coarse]))))
  end

  NF = 2
  LagrangeBasis{1,0,NF}(geometry, fns, pos)
end

gradient(space::LagrangeBasis{1,0}, geo, fns) = NDLCCBasis(geo, fns)
curl(space::LagrangeBasis{1,0}, geo, fns) = RTBasis(geo, fns)

#
# Sclar trace for Laggrange element based spaces
#
function strace(X::LagrangeBasis{1,0}, geo, fns::Vector)
    # dimension of the boundary
    n = dimension(geo)
    # degree of the space
    d = 1
    # number of local shape functions
    NF = binomial(n+d,n)

    trpos = deepcopy(positions(X))
    LagrangeBasis{1,0,NF}(geo, fns, trpos)
end

ext: jl


text:
module ODBC

using DataStreams, DataFrames, NullableArrays, CategoricalArrays, WeakRefStrings

if VERSION < v"0.5.0-dev+4267"
    if OS_NAME == :Windows
        const KERNEL = :NT
    else
        const KERNEL = OS_NAME
    end

    @eval is_apple()   = $(KERNEL == :Darwin)
    @eval is_linux()   = $(KERNEL == :Linux)
    @eval is_bsd()     = $(KERNEL in (:FreeBSD, :OpenBSD, :NetBSD, :Darwin, :Apple))
    @eval is_unix()    = $(is_linux() || is_bsd())
    @eval is_windows() = $(KERNEL == :NT)
else
    const KERNEL = Sys.KERNEL
end

if is_unix()
    using DecFP
end

if !isdefined(Core, :String)
    typealias String UTF8String
end

if !isdefined(Base, :unsafe_wrap)
    unsafe_wrap{A<:Array}(::Type{A}, ptr, len, own) = pointer_to_array(ptr, len, own)
end

if !isdefined(Base, :transcode)
    transcode(::Type{UInt8}, dat) = Base.encode_to_utf8(eltype(dat), dat, length(dat))
end

if !isdefined(Base, Symbol("@static"))
     macro static(ex)
        if isa(ex, Expr)
            if ex.head === :if
                cond = eval(current_module(), ex.args[1])
                if cond
                    return esc(ex.args[2])
                elseif length(ex.args) == 3
                    return esc(ex.args[3])
                else
                    return nothing
                end
            end
        end
        throw(ArgumentError("invalid @static macro"))
    end
end

export Data, DataFrame

include("API.jl")

"just a block of memory; T is the element type, `len` is total # of **bytes** pointed to, and `elsize` is size of each element"
type Block{T}
    ptr::Ptr{T}    # pointer to a block of memory
    len::Int       # total # of bytes in block
    elsize::Int    # size between elements in bytes
end

"""
Block allocator:
    -Takes an element type, and number of elements to allocate in a linear block
    -Optionally specify an extra dimension of elements that make up each element (i.e. container types)
"""
function Block{T}(::Type{T}, elements::Int, extradim::Integer=1)
    len = sizeof(T) * elements * extradim
    block = Block{T}(convert(Ptr{T}, Libc.malloc(len)), len, sizeof(T) * extradim)
    finalizer(block, x->Libc.free(x.ptr))
    return block
end

# used for getting messages back from ODBC driver manager; SQLDrivers, SQLError, etc.
Base.string(block::Block, len::Integer) = String(transcode(UInt8, unsafe_wrap(Array, block.ptr, len, false)))

type ODBCError <: Exception
    msg::String
end

const BUFLEN = 1024

function ODBCError(handle::Ptr{Void}, handletype::Int16)
    i = Int16(1)
    state = ODBC.Block(ODBC.API.SQLWCHAR, 6)
    native = Ref{ODBC.API.SQLINTEGER}()
    error_msg = ODBC.Block(ODBC.API.SQLWCHAR, BUFLEN)
    msg_length = Ref{ODBC.API.SQLSMALLINT}()
    while ODBC.API.SQLGetDiagRec(handletype, handle, i, state.ptr, native, error_msg.ptr, BUFLEN, msg_length) == ODBC.API.SQL_SUCCESS
        st  = string(state, 5)
        msg = string(error_msg, msg_length[])
        println("[ODBC] $st: $msg")
        i += 1
    end
    return true
end

#Macros to to check if a function returned a success value or not
macro CHECK(handle, handletype, func)
    str = string(func)
    quote
        ret = $func
        ret != ODBC.API.SQL_SUCCESS && ret != ODBC.API.SQL_SUCCESS_WITH_INFO && ODBCError($handle, $handletype) &&
            throw(ODBCError("$($str) failed; return code: $ret => $(ODBC.API.RETURN_VALUES[ret])"))
        nothing
    end
end

"List ODBC drivers that have been installed and registered"
function listdrivers()
    descriptions = String[]
    attributes   = String[]
    driver_desc = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    driver_attr = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDrivers(ENV, dir, driver_desc.ptr, BUFLEN, desc_length, driver_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(driver_desc, desc_length[]))
        push!(attributes,   string(driver_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"List ODBC DSNs, both user and system, that have been previously defined"
function listdsns()
    descriptions = String[]
    attributes   = String[]
    dsn_desc    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC.API.SQLSMALLINT}()
    dsn_attr    = Block(ODBC.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC.API.SQLSMALLINT}()
    dir = ODBC.API.SQL_FETCH_FIRST
    while ODBC.API.SQLDataSources(ENV, dir, dsn_desc.ptr, BUFLEN, desc_length, dsn_attr.ptr, BUFLEN, attr_length) == ODBC.API.SQL_SUCCESS
        push!(descriptions, string(dsn_desc, desc_length[]))
        push!(attributes,   string(dsn_attr, attr_length[]))
        dir = ODBC.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"""
A DSN represents an established ODBC connection.
It is passed to most other ODBC methods as a first argument
"""
type DSN
    dsn::String
    dbc_ptr::Ptr{Void}
    stmt_ptr::Ptr{Void}
    stmt_ptr2::Ptr{Void}
end

Base.show(io::IO,conn::DSN) = print(io, "ODBC.DSN($(conn.dsn))")

"""
Construct a `DSN` type by connecting to a valid ODBC DSN or by specifying a valid connection string.
Takes optional 2nd and 3rd arguments for `username` and `password`, respectively.
1st argument `dsn` can be either the name of a pre-defined ODBC DSN or a valid connection string.
A great resource for building valid connection strings is [http://www.connectionstrings.com/](http://www.connectionstrings.com/).
"""
function DSN(dsn::AbstractString, username::AbstractString=String(""), password::AbstractString=String(""); driver_prompt::Integer=ODBC.API.SQL_DRIVER_NOPROMPT)
    dbc = ODBC.ODBCAllocHandle(ODBC.API.SQL_HANDLE_DBC, ODBC.ENV)
    dsns = ODBC.listdsns()
    found = false
    for d in dsns[:,1]
        dsn == d && (found = true)
    end
    if found
        @CHECK dbc ODBC.API.SQL_HANDLE_DBC ODBC.API.SQLConnect(dbc, dsn, username, password)
    else
        dsn = ODBCDriverConnect!(dbc, dsn, driver_prompt % UInt16)
    end
    stmt = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    stmt2 = ODBCAllocHandle(ODBC.API.SQL_HANDLE_STMT, dbc)
    conn = DSN(dsn, dbc, stmt, stmt2)
    return conn
end

"disconnect a connected `DSN`"
function disconnect!(conn::DSN)
    ODBCFreeStmt!(conn.stmt_ptr)
    ODBCFreeStmt!(conn.stmt_ptr2)
    ODBC.API.SQLDisconnect(conn.dbc_ptr)
    return nothing
end

type Statement
    dsn::DSN
    stmt::Ptr{Void}
    query::String
    task::Task
end

"An `ODBC.Source` type executes a `query` string upon construction and prepares data for streaming to an appropriate `Data.Sink`"
type Source <: Data.Source
    schema::Data.Schema
    dsn::DSN
    query::String
    columns::Vector{Any}
    status::Int
    rowsfetched::Ref{ODBC.API.SQLLEN}
    rowoffset::Int
    boundcols::Vector{Any}
    indcols::Vector{Vector{ODBC.API.SQLLEN}}
    sizes::Vector{ODBC.API.SQLULEN}
    ctypes::Vector{ODBC.API.SQLSMALLINT}
    jltypes::Vector{DataType}
end

Base.show(io::IO, source::Source) = print(io, "ODBC.Source:\n\tDSN: $(source.dsn)\n\tstatus: $(source.status)\n\tschema: $(source.schema)")

include("Source.jl")
include("Sink.jl")

function __init__()
    global const ENV = ODBC.ODBCAllocHandle(ODBC.API.SQL_HANDLE_ENV, ODBC.API.SQL_NULL_HANDLE)
end

# used to 'clear' a statement of bound columns, resultsets,
# and other bound parameters in preparation for a subsequent query
function ODBCFreeStmt!(stmt)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_CLOSE)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_UNBIND)
    ODBC.API.SQLFreeStmt(stmt,ODBC.API.SQL_RESET_PARAMS)
end

end #ODBC module


text:
module X

# TODO: add Solution type that keeps the stuff we are currently returning as
# well as other info such as execution time, versioninfo() stats, date/time of
# computation ect.

using JLD
using CompEcon

const DATA_FILE = "bfz_output.jld"

if !isfile(DATA_FILE)
    f = jldopen(DATA_FILE, "w")

    # create model map
    mm = Dict{NTuple{12,Float64}, Int}()
    write(f, "model_map", mm)

    # close the file
    close(f)
end

immutable EZAgent
    ρ::Float64
    α::Float64
    β::Float64
end

_unpack(a::EZAgent) = (a.ρ, a.α, a.β)

immutable Exog
    A::Float64
    B::Float64
    τ::Float64
    φᵥ::Float64
    vbar::Float64
    lgbar::Float64
    ϵp::Matrix{Float64}
    Π::Vector{Float64}
    xp::Matrix{Float64}
    vp::Matrix{Float64}
    gp::Matrix{Float64}
end

function step(A::Float64, B::Float64, φᵥ::Float64, vbar::Float64, τ::Float64, x,
              v, ϵ1, ϵ2)
    xp = A*x .+ B*sqrt(v).*ϵ1
    vp = (1-φᵥ)*vbar .+ φᵥ*v .+ τ.*ϵ2
    xp, vp
end

step(ex::Exog, x, v, ϵ1=ex.ϵp[:, 1], ϵ2=ex.ϵp[:, 2]) =
    step(ex.A, ex.B, ex.φᵥ, ex.vbar, ex.τ, x, v, ϵ1, ϵ2)

function Exog(A::Float64, B::Float64, τ::Float64, φᵥ::Float64, vbar::Float64,
              lgbar::Float64, nϵ1::Int, nϵ2::Int, grid::Matrix{Float64})

    # extract x_t and v_t from grid
    x = grid[:, 2]
    v = grid[:, 3]

    # construct exogenous process
    ϵp, Π = qnwnorm([nϵ1, nϵ2], [0.0, 0.0], eye(2))
    xp, vp = step(A, B, φᵥ, vbar, τ, x, v, ϵp[:, 1]', ϵp[:, 2]')
    gp = exp(lgbar + xp)

    # package -- note, transpose (⋅)p variables so 1 state = 1 column
    Exog(A, B, τ, φᵥ, vbar, lgbar, ϵp, Π, xp', vp', gp')
end

function simulate(ex::Exog; T::Int=25_000, v0::Float64=0.001, x0::Float64=0.0)
    A, B, φᵥ, τ, vbar = ex.A, ex.B, ex.φᵥ, ex.τ, ex.vbar
    v = Array(Float64, T)
    x = Array(Float64, T)

    # set initial conditions
    v[1] = v0
    x[1] = x0

    # constant in law of motion for v
    c = (1 - φᵥ) * vbar

    @inbounds for t=1:T-1
        x[t+1], v[t+1] = step(ex, x[t], v[t], rand(), rand())
    end

    x, v
end

# ----- #
# Model #
# ----- #

immutable BFZ
    agent::EZAgent
    η::Float64
    ν::Float64
    δ::Float64
    grid::Matrix{Float64}
    grid_transpose::Matrix{Float64}  # for iterating state by state over cols
    basis::Basis{3}
    basis_struc::BasisStructure
    exog::Exog
end

function BFZ(;ρ::Real=-1, α::Real=-9, β::Real=0.9995,
              η::Real=1/3, ν::Float64=0.1, δ::Real=0.025,
              lgbar::Real=0.004, e::Real=1.0, A::Real=0., B::Real=1.,
              vbar::Real=0.01^2, φv::Real=0.975, τ::Real=(vbar/3)*sqrt(1-φv^2),
              nk::Int=30, nx::Int=8, nv::Int=8, nϵ1::Int=3, nϵ2::Int=3)

    basis = Basis(SplineParams(collect(linspace(8.0, 32.0, nk)), 0, 3),
                  SplineParams(collect(linspace(-.04, .04, nx)), 0, 1),
                  SplineParams(collect(linspace(1e-4, 4e-4, nv)), 0, 1))

    grid = nodes(basis)[1]

    exog = Exog(A, B, τ, φv, vbar, lgbar, nϵ1, nϵ2, grid)

    basis_struc = BasisStructure(basis, Direct())

    # Create the agent
    agent = EZAgent(ρ, α, β)

    # special case for ν. Make sure it isn't negative and then clamp it to be no
    # smaller than 1e-6. ν=0 ⟶ production is cobb-Douglass. Numerically we can't
    # let ν=0, but if we set it very close to zero then we approximate
    # Cobb-Douglass really well. With some testing I found 1e-8 was optimal
    ν < 0 && throw(ArgumentError("ν cannot be negative"))
    ν = clamp(ν, 1e-8, Inf)

    # package everything up and return
    BFZ(agent, η, ν, δ, grid, grid', basis, basis_struc, exog)
end

function BFZ(params::NTuple{12,Float64})
    ρ, α, β, ν, η, δ, lgbar, A, B, vbar, φᵥ, τ = params
    BFZ(ρ=ρ, α=α, β=β, ν=ν, η=η, δ=δ, lgbar=lgbar, A=A, B=B, vbar=vbar, φv=φᵥ)
end

_params(m::BFZ) = (m.agent.ρ, m.agent.α, m.agent.β, m.ν, m.η, m.δ, m.exog.lgbar,
                   m.exog.A, m.exog.B, m.exog.vbar, m.exog.φᵥ, m.exog.τ)

const DEFAULT_PARAMS = _params(BFZ())

function Base.writemime(io::IO, ::MIME"text/plain", m::BFZ)
    m = """\
    BFZ model with the following parameters:

    ρ: $(m.agent.ρ)
    α: $(m.agent.α)
    β: $(m.agent.β)
    ν: $(m.ν)
    ω(η): $(m.η)
    δ: $(m.δ)
    lgbar: $(m.exog.lgbar)
    A: $(m.exog.A)
    B: $(m.exog.B)
    vbar: $(m.exog.vbar)
    φᵥ: $(m.exog.φᵥ)
    τ: $(m.exog.τ)
    """
    println(m)
end

_Y(m, k) = (m.η*k.^(m.ν) + (1-m.η)).^(1/m.ν)
_VF(m, c, μ) = ((1-m.agent.β)*c^m.agent.ρ + m.agent.β*μ^m.agent.ρ)^(1/m.agent.ρ)
_CE(m, jp, gp) = dot((gp.*jp).^(m.agent.α), m.exog.Π).^(1/m.agent.α)

function solve_this_state(m::BFZ, i::Int, xv_mat::AbstractMatrix, coefs::Vector,
                          y::Real, j::Real, jk::Real, gp::Vector)
    k, x, v = m.grid_transpose[:, i]
    ρ, α, β = _unpack(m.agent)
    ν, η, δ = m.ν, m.η, m.δ

    N = length(gp)

    # solve for c and k'
    c = (jk/(j^(1-ρ) * (1-β) * (y^(1-ν)*η*k^(ν-1) + 1 - δ))).^(1/(ρ - 1))
    kp = (y + (1-δ)*k - c) ./ gp

    # compute j' to evaluate μ
    kp_basis_mat = CompEcon.evalbase(m.basis.params[1], kp, 0)[1]
    jp = CompEcon.row_kron(xv_mat, kp_basis_mat) * coefs

    μ = _CE(m, jp, gp)
    jm = _VF(m, c, μ)

    return jm, c, kp
end

function basis_mat_xv(m::BFZ, x::Vector, v::Vector)
    mat_v = CompEcon.evalbase(m.basis.params[3], v, 0)[1]
    mat_x = CompEcon.evalbase(m.basis.params[2], x, 0)[1]
    CompEcon.row_kron(mat_v, mat_x)::typeof(mat_v)
end

function solve(m::BFZ, re_compute::Bool=false; kwargs...)
    jldopen(DATA_FILE, "r+") do f
        model_map = read(f, "model_map")
        the_params = _params(m)

        if haskey(model_map, the_params) && !re_compute
            grp_key = model_map[the_params]
            grp = f[string(grp_key)]
            coefs = read(grp, "coefs")
            j = read(grp, "j")
            c = read(grp, "c")
            kp = read(grp, "kp")
            println("Returning cached solution")
        else
            if haskey(model_map, the_params)
                # delete the group so we can re-compute
                delete!(f, "$(model_map[the_params])")
                grp_key = model_map[the_params]
            else
                grp_key = length(model_map) > 0 ? maximum(values(model_map)) + 1 : 1
            end

            coefs, j, c, kp = solve_ecm(m; kwargs...)

            for (nm, val) in zip([:j, :coefs, :c, :kp], (j, coefs, c, kp))
                write(f, "$grp_key/$nm", val)
            end

            # update the model map
            model_map[the_params] = grp_key
            delete!(f, "model_map")
            write(f, "model_map", model_map)
        end

        coefs, j, c, kp, grp_key
    end
end

function solve_ecm(m::BFZ; ξ::Float64=0.4, tol::Float64=1e-10, maxiter::Int=5000)
    # initialize J
    ρ, α, β = _unpack(m.agent)
    ν, η, δ = m.ν, m.η, m.δ
    M = size(m.grid, 1)
    N = size(m.exog.xp, 1)

    jm_old = (1-β)^(1/ρ)*(_Y(m, m.grid[:, 1]) + (1 - δ)*m.grid[:, 1]) ./ 1000.0
    coefs = funfitxy(m.basis, m.basis_struc, jm_old)[1]

    J_k_basis_struc = BasisStructure(m.basis, Direct(),
                                     nodes(m.basis)[1], [1 0 0])

    xv_mats = [basis_mat_xv(m, m.exog.xp[:, i], m.exog.vp[:, i]) for i=1:M]
    err = 1.0
    it = 0

    jm = similar(jm_old)
    y_all = _Y(m, m.grid[:, 1])

    local out  # declare out local so we can return it after the loop
    while err > tol && it < maxiter
        it += 1
        # compute J at each state
        J = funeval(coefs, m.basis_struc, [0 0 0])
        J_k = funeval(coefs, J_k_basis_struc, [1 0 0])

        # solve at each state in parallel
        ssi(i) = solve_this_state(m, i, xv_mats[i], coefs, y_all[i], J[i],
                                  J_k[i], m.exog.gp[:, i])
        out = map(ssi, 1:M)

        # make sure we have jm::Vector{Float64}
        jm = Float64[i[1] for i in out]

        # now let's get ourselves a new coefficient vector
        coef_hat = funfitxy(m.basis, m.basis_struc, jm)[1]
        coefs = (1-ξ)*coef_hat + ξ*coefs
        err = maxabs(jm - jm_old)
        copy!(jm_old, jm)

        @show it, err

    end

    # extract solution
    j = Float64[i[1] for i in out]
    c = Float64[i[2] for i in out]
    kp = hcat(Vector{Float64}[i[3] for i in out]...)

    return coefs, j, c, kp
end

function simulate_solution(m::BFZ, coefs::Vector, c::Vector; T::Int=75_000,
                           burn::Int=floor(Int, T/10))
    # simulate exog process
    x, v = simulate(m.exog; T=T)

    # allocate memory for simulated c and k
    c_sim = Array(Float64, T-1)
    k_sim = Array(Float64, T)
    y_sim = Array(Float64, T-1)

    # "arbitrary" initial condition for k. Mean of k on the grid
    k_sim[1] = mean(m.grid[:, 1])

    # get coefficients for approximation of c
    c_coefs = get_coefs(m.basis, m.basis_struc, c)

    gbar = exp(m.exog.lgbar)

    # simulate c, j forward
    @inbounds for t=1:T-1
        c_sim[t] = funeval(c_coefs, m.basis, [k_sim[t] x[t] v[t]])[1]
        gp = gbar*exp(x[t+1])
        y_sim[t] = _Y(m, k_sim[t])
        k_sim[t+1] = (y_sim[t] + (1-m.δ)*k_sim[t] - c_sim[t]) / gp
    end
    x, v, c_sim, k_sim, y_sim = (x[burn:end], v[burn:end], c_sim[burn:end],
                                 k_sim[burn:end], y_sim[burn:end])

    c_coefs, x, v, c_sim, k_sim, y_sim
end

function euler_errors_sim(m::BFZ, coefs::Vector, c_coefs::Vector, x_sim, v_sim,
                          c_sim, k_sim, y_sim)
    ρ, α, β = _unpack(m.agent)
    ν, δ, η = m.ν, m.δ, m.η

    # array to hold the errors, which are computed by solving
    # (c*(1-err))^{ρ-1} = RHS for err ⟶ err = 1- RHS^{1/(ρ-1)}/c
    err_sim = similar(c_sim)

    gbar = exp(m.exog.lgbar)

    for t=1:length(c_sim)-1
        # extract stuff for this state
        ct, kt, yt, xt, vt = c_sim[t], k_sim[t], y_sim[t], x_sim[t], v_sim[t]

        # push exog forward using quadrature nodes
        xp, vp = step(m.exog, xt, vt)
        gp = gbar*exp(xp)

        # get kp; then yp, cp, and jp; then μ
        kp = (yt + (1-m.δ)*kt - ct) ./ gp
        yp = _Y(m, kp)
        bs = BasisStructure(m.basis, Direct(), [kp xp vp], [0 0 0])
        cp = funeval(c_coefs, bs, [0 0 0])
        jp = funeval(coefs, bs, [0 0 0])
        μ = _CE(m, jp, gp)

        # now I can eval RHS
        to_E = gp.^(α-1) .* jp.^(α-ρ) .* cp.^(ρ-1) .* (η*yp.^(1-ν).*kp.^(ν-1)+1-δ)
        RHS = β * μ^(ρ-α) * dot(to_E, m.exog.Π)

        # now fill in err_sim
        err_sim[t] = 1 - RHS^(1/(ρ-1))/ct

    end
    err_sim
end

function main()
    m = BFZ()
    solve(m)
end

end  # module

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
@testset "Function tests" begin

@testset "Logs"  begin

    var_states = Dict{Int64,Vector{Int64}}()
    var_changes = Dict{Int64,Vector{Tuple{Symbol, Int64, Int64, Int64}}}()
    children = CS.TreeLogNode[]
    l1 = CS.TreeLogNode(0,:Open,0,0,0,0,var_states, var_changes, children)
    l2 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children)
    @test !CS.same_logs(l1,l2)

    # different children order
    children = CS.TreeLogNode[]
    tln1 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln2 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, CS.TreeLogNode[])
    push!(children, tln1)
    push!(children, tln2)
    l1 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children[[1,2]])
    l2 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children[[2,1]])
    @test !CS.same_logs(l1,l2)

    # missing one child
    children = CS.TreeLogNode[]
    tln1 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln2 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, CS.TreeLogNode[])
    push!(children, tln1)
    push!(children, tln2)
    l1 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children[[1,2]])
    l2 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children[[1]])
    @test !CS.same_logs(l1,l2)

    # 2 layer children different struture
    children1 = CS.TreeLogNode[]
    tln13 = CS.TreeLogNode(0,:Open,0,3,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln11 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, [tln13])
    tln12 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, CS.TreeLogNode[])
    push!(children1, tln11)
    push!(children1, tln12)
    children2 = CS.TreeLogNode[]
    tln23 = CS.TreeLogNode(0,:Open,0,3,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln21 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln22 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, [tln23])
    push!(children2, tln21)
    push!(children2, tln22)
    l1 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children1)
    l2 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children2)
    @test !CS.same_logs(l1,l2)

    # 2 layer children different position
    children1 = CS.TreeLogNode[]
    tln13 = CS.TreeLogNode(0,:Open,0,3,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln14 = CS.TreeLogNode(0,:Open,0,4,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln11 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, [tln13,tln14])
    tln12 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, CS.TreeLogNode[])
    push!(children1, tln11)
    push!(children1, tln12)
    children2 = CS.TreeLogNode[]
    tln23 = CS.TreeLogNode(0,:Open,0,3,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln24 = CS.TreeLogNode(0,:Open,0,4,0,0,var_states, var_changes, CS.TreeLogNode[])
    tln21 = CS.TreeLogNode(0,:Open,0,1,0,0,var_states, var_changes, [tln24, tln23])
    tln22 = CS.TreeLogNode(0,:Open,0,2,0,0,var_states, var_changes, CS.TreeLogNode[])
    push!(children2, tln21)
    push!(children2, tln22)
    l1 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children1)
    l2 = CS.TreeLogNode(0,:Closed,0,0,0,0,var_states, var_changes, children2)
    @test !CS.same_logs(l1,l2)
end

end


text:
using Plots
using JSON

jstxt = open("result.json") do file
  read(file, String)
end

data = JSON.parse(jstxt)

function to_sec(x)
  strm,strs = split(x, "m")
  m=parse(Float64,strm)
  s=parse(Float64,strs[1:length(strs)-1])
  m*60+s
end

envs = [ "2.6.2", "2.6.2/--jit", "jruby-9.2.6.0", "truffleruby-1.0.0-rc14" ]
k = sort(collect(keys(data[envs[1]])))
println(k)
y = map( env->begin
  values = map( x->data[env][x], k )
  return map( to_sec, values )
end, envs )

plot( y,
  legend = :topleft,
  label = envs,
  xlabel = "Argument of fibo",
  ylabel = "result of time (real)",
  xticks = (1:length(k), k),
  yscale=:log10
)

savefig("plot.pdf")

max_stars_repo_name: nabetani/qiita-scribbles

max_issues_repo_name: nabetani/qiita-scribbles

max_forks_repo_name: nabetani/qiita-scribbles


text:
import RowEchelon

#=
Homework 5
Section 2.2
Problem 5
=#

# Adapt system into a matrix
A = [8 6 ; 5 4]
b = [2 ;-1]

# Find Ax = b with x = invA(b)
invA = inv(A)
x = invA*b

println("The solution is:")
print("x = ")
println(x)


text:
using DataTypesBasic

# Monoid support for Exceptions
TypeClasses.neutral(::Type{<:Exception}) = MultipleExceptions()
TypeClasses.combine(e1::Exception, e2::Exception) = MultipleExceptions(e1, e2)


text:
import Pkg
Pkg.rm(["Grids", "GridGames"])
Pkg.develop(["Grids", "GridGames"])

using MuZero
using GridGames
using Base.Threads: @spawn
cd(@__DIR__)
include("game.jl")
include("muzero_config.jl")
config = getconfig(Tac)
include("inference_tests.jl")
pv_tests = makepolicyvaluetests(Tac)
pv_tests.period = 100
learner = Learner(
    config;
    uuid="435c256d-1c60-4e07-ac03-a2430ae7b7bd",
    dir="~/.julia/dev/MuZeroModels/tictactoe/",
    tests=[pv_tests],
    schedule=tacschedule!
)
getpvtestresult(learner, 3335)
getpvtestresult(learner, 3335) |> svg
length(learner.replay_buffers)
mzgame = learner.replay_buffers[1]
mzgame.game.states[end]
target = getstate(Tac([1,4,8,2])) #6,3,6
found = findmuzerogames(learner, target)

length(learner.replay_buffers)
learner.stats.policy_test_fails
#40 fails: 6,8,25,31,47,51,57,78,116,138,139,142,143,180,185,186,187,201,202,218,225,259,284,300,323,334,418,512,514,538,628,725,778,1098,1114,1188,1413,2120,2760,3335
showresultsfull(learner; mean=100)
get_steps(learner)
acc, fails, value_mse = runtests(pv_tests, training_network(learner))

using Flux
test_net = new_test_network(learner.config)
testmode!.(get_nets(test_net), true)
updatetestparams!(learner.train, test_net)
testmode!.(get_nets(test_net), true)
acc2, fails2, value_mse2 = runtests(pv_tests, test_net)

#run all tests one by one! (instead of in batch mode)
fails = []
for pvtest in pv_tests
    out = initial_inference(test_net, pvtest.state) |> NetworkOutput
    append!(fails, getpolicyfails(pvtest, out))
end

fit!(learner; steps=500)
get_steps(learner)
learner.stats.policy_test_fails
# 45 fails: 6,25,31,47,51,57,78,116,121,128,131,138,139,142,172,180,185,187,201,218,225,259,284,323,334,418,512,514,538,628,725,761,778,839,1098,1114,1188,1318,1413,1606,2066,2120,2760,3254,3335

getpvtestresult(learner, 3254) |> svg
target2 = getstate(Tac([1,5,4,6,2,9]))
found = findmuzerogames(
    learner,
    getpolicyvaluetests(learner)[3254].state
)

fit!(learner; steps=90)
# 47 fails: 6,25,31,47,51,57,78,116,121,128,131,139,142,167,172,180,185,187,201,218,225,238,259,284,285,323,334,418,512,514,538,628,778,834,839,1098,1114,1188,1318,1413,1606,1861,2120,2760,3254,3286,3335
found = findmuzerogames(
    learner,
    getpolicyvaluetests(learner)[3286].state
)
getpolicy(
    learner,
    getpolicyvaluetests(learner)[3254].state
)
getpvtestresult(learner, 3254) |> svg


text:
"
 startup julia script

 * create link
 ln -nfs ~/dotfiles/startup.jl ~/.julia/config/startup.jl
"

ENV["JULIA_EDITOR"] = "vim"
try
    ENV["PYTHON"]=Sys.which("python3")
catch
    ENV["PYTHON"]=""
end

atreplinit() do repl
    # for OhMyREPL
    try
        @eval using OhMyREPL
    catch e
        @warn "error while importing OhMyREPL" e
    end

    repl.options.iocontext[:compact] = false
end

# For Revise.jl
try
    using Revise
catch e
    @warn(e.msg)
end

function recursive_includet(filename)
    already_included = copy(Revise.included_files)
    includet(filename)
    newly_included = setdiff(Revise.included_files, already_included)
    for (mod, file) in newly_included
        Revise.track(mod, file)
    end
end

# For BenchmarkTools
atreplinit() do repl
    try
        @eval using BenchmarkTools
    catch e
        @warn(e.msg)
    end
end

# For PackageCompiler.jl
try
    using PackageCompiler
catch e
    @warn(e.msg)
end

@eval using Pkg

DEFAULT_PKGS = [:OhMyREPL, :Revise, :BenchmarkTools, :PyPlot]

function install_default_pkgs()
    println("Installing default packages to latest ones...")
    for pkg in DEFAULT_PKGS
        Pkg.add(String(pkg))
    end
end

function generate_default_sysimage()

    println("Updating default packages to latest ones...")
    for pkg in DEFAULT_PKGS
        println("Updating $(String(pkg))")
        Pkg.update(String(pkg))
    end

    println("Restore sysimage once...")
    try
        restore_default_sysimage()
    catch e
        println(e.msg)
    end

    create_sysimage(DEFAULT_PKGS;replace_default=true)

    println("Done!!. If you want to restore default sysimage, run `restore_default_sysimage()`")

end

max_stars_repo_name: AtsushiSakai/dotfiles

max_issues_repo_name: AtsushiSakai/dotfiles

max_forks_repo_name: AtsushiSakai/dotfiles


text:
using JSON

get_dataset_path(;dataset_name="dataset_path") = JSON.parsefile("paths.json")[dataset_name]

norm(arr, min, max) = (arr .- min) ./ (max - min)
denorm(arr, min, max) = arr .* (max - min) + min
norm_to_range(arr, min_x, max_x, range_a, range_b) = (range_b - range_a) .* ((arr .- min_x) ./ (max_x - min_x)) .+ range_a

has_inf(arr) = sum(isinf.(arr)) > 0
has_nan(arr) = sum(isnan.(arr)) > 0

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Kolakoski

using Nemo

export ModuleKolakoski
export KolakoskiList, C000002, I000002, L000002

"""

* KolakoskiList, C000002, I000002, L000002
"""
const ModuleKolakoski = ""

"""

Generate the Kolakoski sequence which is the unique sequence over the alphabet ``{1, 2}`` starting with ``1`` and having the sequence of run lengths identical with itself.
"""
C000002() =
    Channel(csize = 10) do c
        x = y = Int(-1)

        while true
            put!(c, [2, 1][(x&1)+1])
            f = y & ~(y + 1)
            x = xor(x, f)
            y = (y + 1) | (f & (x >> 1))
        end
    end

struct KolakoskiSeq
    count::Int
    ch::Channel
    KolakoskiSeq(count) = new(count, C000002())
end

function Base.iterate(I::KolakoskiSeq)
    if I.count == 0
        close(I.ch)
        return nothing
    end
    (take!(I.ch), (0))
end

function Base.iterate(I::KolakoskiSeq, S)
    j = S[1] + 1
    if I.count == j
        close(I.ch)
        return nothing
    end
    (take!(I.ch), (j))
end

Base.length(I::KolakoskiSeq) = I.count
Base.eltype(I::KolakoskiSeq) = Int

"""

Iterate over the first ``n`` Kolakoski numbers.
"""
I000002(n::Int) = KolakoskiSeq(n)

"""

Return the list of the first ``n`` terms of the Kolakoski sequence.
"""
function KolakoskiList(len::Int)
    len ≤ 0 && return []
    generator = C000002()
    L = [take!(generator) for _ ∈ 1:len]
    close(generator)
    return L
end

"""

Return the list of the first ``n`` terms of the Kolakoski sequence.
"""
L000002(n::Int) = KolakoskiList(n)

#START-TEST-########################################################

using Test

function test()

    @testset "Kolakoski" begin
        K = KolakoskiList(100)
        @test K[1] == 1
        @test K[33] == 2
        @test K[72] == 2

        generator = C000002()
        for n ∈ [1, 33, 72]
            k = take!(generator)
            @test K[n] == k
        end
        close(generator)
    end
end

function demo()
    println(KolakoskiList(20))

    generator = C000002()
    o = e = 0
    for n ∈ 1:80
        take!(generator) == 1 ? o += 1 : e += 1
        print(o - e, " ")
    end
    println()
    close(generator)

    for f ∈ I000002(20)
        print(f, ", ")
    end
    println()
    print(L000002(20))
    println()
end

"""

I000002(100000) ::
    0.000035 seconds (31 allocations: 2.969 KiB)
KolakoskiList(10000) ::
    0.086202 seconds (120.03 k allocations: 10.226 MiB, 17.61% gc time)
"""
function perf()
    @time I000002(100000)
    @time KolakoskiList(100000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=

[1, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1]

1 0 -1 0 1 0 1 0 -1 0 -1 -2 -1 0 -1 0 1 0 -1 0 -1 0 1 0 1 0 -1 0
1 0 1 2 1 2 1 0 1 0-1 0 1 0 1 0 -1 0 -1 0 1 0 1 2 1 0 1 0 -1 0 1
0 1 0 -1 0 -1 -2 -1 0 -1 0 1 0 1 0 -1 0 -1 0 1 0

1, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1,

[1, 2, 2, 1, 1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1]

  0.000035 seconds (31 allocations: 2.969 KiB)
  0.086202 seconds (120.03 k allocations: 10.226 MiB, 17.61% gc time)

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
module GeometryRandomFields

export EΓ,Γ,grf,GaussianCovarianceFunction,MatérnCovarianceFunction

include("covariancefunctions.jl")
include("generators.jl")
include("spectralrepresentation.jl")
include("characteristic.jl")

end # module


text:
using Documenter, MakieLayout, Makie, Animations, ColorSchemes, Colors

# don't open windows while generating animations
Makie.AbstractPlotting.inline!(true)

makedocs(
    sitename="MakieLayout.jl",
    pages = [
        "index.md",
        "GridLayout" => "grids.md",
        "LAxis" => "laxis.md",
        "LLegend" => "llegend.md",
        "Other Layoutables" => "layoutables.md",
        "How layouting works" => "layouting.md",
        "Frequently Asked Questions" => "faq.md",
    ],
    format = Documenter.HTML(
            prettyurls = get(ENV, "CI", nothing) == "true"
        )
    )

# struct Local <: Documenter.DeployConfig end
#
# function Documenter.deploy_folder(cfg::Local;
#         repo, devbranch, push_preview, devurl, kwargs...)
#
#     folder = if get(ENV, "PUSH_LOCAL_BUILD", false) == "true"
#         @warn("Setting ENV[\"PUSH_LOCAL_BUILD\"] = \"false\", remember to set it to true for the next push.")
#         ENV["PUSH_LOCAL_BUILD"] = "false"
#         devurl
#     else
#         @warn("Set ENV[\"PUSH_LOCAL_BUILD\"] = \"true\" if you want your local build to be pushed to Github Pages.")
#         nothing
#     end
# end
#
# Documenter.authentication_method(::Local) = Documenter.SSH
#
# function Documenter.documenter_key(::Local)
#     open(readline, expanduser("~/.ssh/documenter_makielayout"))
# end

struct Gitlab <: Documenter.DeployConfig
    commit_branch::String
    pull_request_iid::String
    repo_slug::String
    commit_tag::String
    pipeline_source::String
end

function Gitlab()
    commit_branch = get(ENV, "CI_COMMIT_BRANCH", "")
    pull_request_iid = get(ENV, "CI_EXTERNAL_PULL_REQUEST_IID", "")
    repo_slug = get(ENV, "CI_PROJECT_PATH_SLUG", "")
    commit_tag = get(ENV, "CI_COMMIT_TAG", "")
    pipeline_source = get(ENV, "CI_PIPELINE_SOURCE", "")
    Gitlab(
        commit_branch,
        pull_request_iid,
        repo_slug,
        commit_tag,
        pipeline_source,
    )
end

function Documenter.deploy_folder(cfg::Gitlab; repo, devbranch, push_preview, devurl, kwargs...)

    io = IOBuffer()
    all_ok = true

    println(io, "Gitlab config:\n", cfg)

    subfolder = if cfg.commit_tag != ""
        tag_ok = occursin(Base.VERSION_REGEX, cfg.commit_tag)
        println("tag_ok: ", tag_ok)
        all_ok &= tag_ok

        cfg.commit_tag
    else
        devurl
    end

    key_ok = haskey(ENV, "DOCUMENTER_KEY")
    println(io, "key_ok: ", key_ok)
    all_ok &= key_ok

    @info String(take!(io))

    return all_ok ? subfolder : nothing
end

Documenter.authentication_method(::Gitlab) = Documenter.SSH

function Documenter.documenter_key(::Gitlab)
    return ENV["DOCUMENTER_KEY"]
end

deploydocs(
    repo = "github.com/jkrumbiegel/MakieLayout.jl.git",
    deploy_config = Gitlab(),
)


text:
"""Julia program to implement Bubble Sort
In Bubble Sort we compare the adjacent elements and swap them such that the greater
element among them comes in the right-side. After each iteration, one smallest element
goes towars the end of the array.
"""

function bubble_sort(arr, n)
    # Iterate the entire array
    for i in 1:n
        # Iterate till the concerned part only.
        for j in 1:(n-i)
            # Swap adjacent elements if they are in wrong order
            if(arr[j] > arr[j+1])
                arr[j], arr[j+1] = arr[j+1], arr[j]
            end
        end
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int64, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
bubble_sort(arr, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers do you want to sort? 5
Enter the numbers: 4 2 5 1 3
The numbers in sorted order is [1, 2, 3, 4, 5]

"""


text:
module EuclideanDistanceMatrices

using LinearAlgebra, Statistics
import Pkg
using StatsBase
using TotalLeastSquares
import Convex
using SCS
using MonteCarloMeasurements
using Turing, Distributions
using Turing2MonteCarloMeasurements

export complete_distmat, rankcomplete_distmat, reconstruct_pointset, denoise_distmat, lowrankapprox, procrustes, posterior, align_to_mean

"""
    D̃, S = complete_distmat(D, W, λ = sqrt(count(W .== 0)))

Takes an incomplete squared Euclidean distance matrix `D` and fills in the missing entries indicated by the mask `W`. `W` is a `BitArray` or array of {0,1} with 0 denoting a missing value. Returns the completed matrix and an SVD object that allows reconstruction of the generating point set `X`.

*NOTE* This function is only available after `using Convex, SCS`, and it will run out of memory for large matrices. See `rankcomplete_distmat` for a less accuraet function that works for larger matrices.

# Arguments:
- `D`: The incomplete matrix
- `W`: The mask
- `λ`: Regularization parameter. A higher value enforces better data fitting, which might be required if the number of entries in `D` is very small.

# Example:
```julia
using EuclideanDistanceMatrices, Distances
P = randn(2,40)
D = pairwise(SqEuclidean(), P)
W = rand(size(D)...) .> 0.3 # Create a random mask
W = (W + W') .> 0           # It makes sense for the mask to be symmetric
W[diagind(W)] .= true
D0 = W .* D                 # Remove missing entries

D2, S = complete_distmat(D0, W)

@show (norm(D-D2)/norm(D))
@show (norm(W .* (D-D2))/norm(D))
```

The set of points that created `D` can be reconstructed up to an arbitrary rotation and translation, `X` contains the reconstruction in the `d` first rows, where `d` is the dimension of the point coordinates. To reconstruct `X` using `S`, do
```julia
X  = reconstruct_pointset(S, 2)

# Verify that reconstructed `X` is correct up to rotation and translation
R,t = procrustes(X,P)
P2 = R*X .+ t
norm(P-P2)/norm(P) # Should be small
```

Ref: Algorithm 5 from "Euclidean Distance Matrices: Essential Theory, Algorithms and Applications"
Ivan Dokmanic, Reza Parhizkar, Juri Ranieri and Martin Vetterli https://arxiv.org/pdf/1502.07541.pdf
"""
function complete_distmat(D, W, λ=sqrt(count(W .== 0)))
    @assert all(==(1), diag(W)) "The diagonal is always observed and equal to 0. Make sure the diagonal of W is true"
    @assert all(iszero, diag(D)) "The diagonal of D is always 0"
    n = size(D, 1)
    x = -1/(n + sqrt(n))
    y = -1/sqrt(n)
    V = [fill(y, 1, n-1); fill(x, n-1,n-1) + I(n-1)]
    e = ones(n)
    G = Convex.Variable((n-1, n-1))
    B = V*G*V'
    E = diag(B)*e' + e*diag(B)' - 2*B
    problem = Convex.maximize(tr(G)- λ * norm(vec(W .* (E - D))), [G ∈ :SDP])
    Convex.solve!(problem, SCS.Optimizer)
    Int(problem.status) == 1 || @error problem.status
    B  = Convex.evaluate(B)
    D2 = diag(B)*e' + e*diag(B)' - 2*B
    @info "Data fidelity (norm(W .* (D-D̃))/norm(D))", (norm(W .* (D-D2))/norm(D))
    s  = svd(B)
    D2, s
end

"""
    rankcomplete_distmat(D, W, dim; μ=mean(D[W]), iters=100, tol=1e-6, verbose=true)

Takes an incomplete squared Euclidean distance matrix `D` and fills in the missing entries indicated by the mask `W::BitArray` with 0 denoting a missing value. Returns the completed matrix and an `Eigen` object that allows reconstruction of the generating point set `X`.

This function works for larger matrices than `complete_distmat`, but is much less accurate. See `complete_distmat` for examples.

# Arguments:
- `D`: The incomplete matrix
- `W`: The mask
- `dim`: The dimension of the points that created `D`
- `μ`: Initial value for unobserved entries.
- `iters`: Maximum number of iterations
- `tol`: Tolerance for the change in `D` between iterations.
- `verbose`: print status
"""
function rankcomplete_distmat(D, W, dim; μ=mean(D[W]), iters=100, tol=1e-6, verbose=true)
    @assert all(==(1), diag(W)) "The diagonal is always observed and equal to 0. Make sure the diagonal of W is true"
    @assert all(iszero, diag(D)) "The diagonal of D is always 0"
    n    = size(D, 1)
    D2   = copy(D)
    Dold = copy(D2)
    D2[.!W] .= μ
    local E
    for iter = 1:iters
        D2, E = eigval_th(Symmetric(D2), dim+2)
        D2[W] .= D[W]
        D2[diagind(D2)] .= 0
        D2 .= max.(D2, 0)
        err = sqrt(sum(abs2(d1-d2) for (d1,d2) in zip(Dold, D2)))/norm(D2)
        verbose && @info "Iter $iter Change: $err"
        err < tol && break
        Dold .= D2
    end
    D2, E
end

function eigval_th(D, r)
    E = eigen(D)
    E.values[1:end-r+1] .= 0
    E.vectors*Diagonal(E.values)*E.vectors', E
end

"""
    reconstruct_pointset(D, dim)

Takes a squared distance matrix or the SVD of one and reconstructs the set of points embedded in dimension `dim` that generated `D; up to a translation and rotation/reflection. See `procrustes` for help with aligning the result to a collection of anchors.

The algorithm used if `D` is a distance matrix is called Multidimensional Scaling (MDS) https://en.wikipedia.org/wiki/Multidimensional_scaling
"""
function reconstruct_pointset(D::AbstractMatrix, dim)
    n = size(D,1)
    J = I - fill(1/n, n, n)
    G = -1/2 * J*D*J
    E = eigen(Symmetric(G))
    reconstruct_pointset(E, dim)
end

function reconstruct_pointset(E::Eigen, dim)
    Diagonal(sqrt.(max.(E.values[end-dim+1:end], 0)))*E.vectors[:, end-dim+1:end]'
end

function reconstruct_pointset(S::SVD,dim)
    X  = Diagonal(sqrt.(S.S[1:dim]))*S.Vt[1:dim, :]
end

"""
    denoise_distmat(D, dim, p = 2)

Takes a noisy squared distance matrix and returns a denoised version. `p` denotes the "norm" used in measuring the error. `p=2` assumes that the error is Gaussian, whereas `p=1` assumes that the error is large but sparse.

# Arguments:
- `dim`: The dimension of the points that generated `D`
"""
function denoise_distmat(D, dim, p=2)
    if p == 2
        s = svd(D)
        return lowrankapprox(s, dim+2)
    elseif p == 1
        A,E,s,sv = rpca(D, nonnegA=true)
        return lowrankapprox(s, dim+2)
    else
        throw(ArgumentError("p must be 1 or 2"))
    end
end

function lowrankapprox(s::SVD, r)
    @views s.U[:,1:r] * Diagonal(s.S[1:r]) * s.Vt[1:r,:]
end
lowrankapprox(D, r) = D = lowrankapprox(svd(D), r)

"""
    R,t = procrustes(X, Y)

Find rotation matrix `R` and translation vector `t` such that `R*X .+ t ≈ Y`
"""
function procrustes(X,Y)
    mX = mean(X, dims=2)
    mY = mean(Y, dims=2)
    Xb,Yb = X .- mX, Y .- mY
    s = svd(Xb*Yb')
    R = s.V*s.U'
    R, mY - R*mX
end

"""
    part, result = posterior(locations::AbstractMatrix, distances, args...; nsamples = 3000, sampler = NUTS(), σL = 0.3, σD = 0.3)

Compute the full Bayesian posterior over locations given noisy measurements of both locations and distances.
Missing values can be indicated by using either the `missing` element, or by providing a rough estimate and setting a very large standard deviation for the guessed element.

# Arguments:
- `locations`: A matrix with one point in each column.
- `distances`: A vector of tuples `(i,j,d)` where `i,j` denotes the indices between the distance measured is `d` (not squared).
- `args`: If any provided, they are sent to `optimize`. Can be used to supply a choice of optimizer of `MAP` estimation.
- `nsamples`: How many samples to draw for MCMC sampling.
- `sampler`: Any of the samplers or `MAP/MLE` supported by Turing
- `σL`: The noise std in the locations.
- `σD`: The noise std in the distances.
- `tdoa`: Indicates whether or not `distances` contains TDOA measurements. See note below.

`part` contains `part.P` and `part.d` with the posterior over locations and distances.

## TDOA
When using `tdoa=true`, `distances` is expected to contain tuples `(i,j,tdoa)`, and the last column in `locations` is expected to contain an estimate of the source location. Note that the standard deviation for locations can be a vector, aiding in situations when this location is known with more or less accuracy than the sensors.

# Example
```julia
N = 10

σL = 0.1
σD = 0.01

P = randn(2,N)
Pn = P + σL*randn(size(P))
D = pairwise(Euclidean(), P, dims=2)
Dn = D + σD*randn(size(D))
Dn[diagind(Dn)] .= 0

distances = []
noisy_distances = []
p = 0.5
for i = 1:N
    for j = i+1:N
        if rand() < p
            push!(distances, (i,j,D[i,j]))
            push!(noisy_distances, (i,j,Dn[i,j]))
        end
    end
end
@show length(distances)
@show expected = p*((N^2-N)÷2)

part, chain = posterior(
    Pn,
    noisy_distances;
    nsamples = 2000,
    sampler = NUTS(),
    σL = σL,
    σD = σD
)

norm(mean.(part.P) - P) < norm(Pn - P)

scatter(part.P[1,:], part.P[2,:], markersize=6, layout=2, sp=1)
scatter!(P[1,:], P[2,:], lab="True positions", sp=1)
scatter!(Pn[1,:], Pn[2,:], lab="Measured positions", sp=1)
bar!(getindex.(distances, 3), sp=2, lab="True dist")
bar!(getindex.(noisy_distances, 3), sp=2, lab="Measured dist", alpha=0.7)
scatter!(part.d, sp=2, seriestype=:scatter) |> display
```
"""
function posterior(
    locations::AbstractMatrix{S},
    distances,
    args...;
    tdoa = false,
    nsamples = 3000,
    sampler = NUTS(),
    σL = 0.3,
    σD = 0.3,
) where S
    dim, N = size(locations)
    Nd = length(distances)

    if tdoa
        d = abs.(getindex.(distances, 3))
        Turing.@model model(locations, distances, d, ::Type{T} = Float64) where {T} = begin
            P0 ~ MvNormal(vec(locations), σL) # These denote the true locations
            P = reshape(P0, dim, N)
            dh = Vector{T}(undef, Nd) # These are the predicted distance measurements
            for ind in eachindex(distances)
                (i,j,_) = distances[ind]
                di = norm(P[:,i] - P[:,end]) # Distance from source to i
                dj = norm(P[:,j] - P[:,end]) # Distance from source to j
                dh[ind] = abs(di-dj) # This is the predicted TDOA given the posterior locations
            end
            d ~ MvNormal(dh, σD) # Observe TDOAs
        end
    else
        d = (getindex.(distances, 3))
        Turing.@model model(locations, distances, d, ::Type{T} = Float64) where {T} = begin
            P0 ~ MvNormal(vec(locations), σL) # These denote the true locations
            P = reshape(P0, dim, N)
            dh = Vector{T}(undef, Nd) # These are the predicted distance measurements
            for ind in eachindex(distances)
                (i,j,_) = distances[ind]
                dh[ind] = norm(P[:,i] - P[:,j]) # This is the predicted Euclidean distance given the posterior location
            end
            d ~ MvNormal(dh, σD) # Observe distances
        end
    end

    m = model(locations, distances, d)

    if sampler isa Turing.Inference.InferenceAlgorithm
        @info "Starting sampling (this might take a while)"
        @time chain = sample(m, sampler, nsamples, drop_warmup=true)
        nt = Particles(chain)
        P = reshape(nt.P0, dim, N)
        nt = (nt..., P=P, d=[norm(P[:,i] - P[:,j]) for i in 1:N, j in 1:N])
        @info "Done"
        return nt, chain
    elseif typeof(sampler) ∈ (Turing.MAP, Turing.MLE)
        res = optimize(m, sampler, args...)
        c = StatsBase.coef(res)
        C = StatsBase.vcov(res)
        names = StatsBase.params(res)

        Pinds = findfirst.([==("P0[$i]") for i in eachindex(locations)], Ref(names))
        # dinds = findfirst.([==("d[$i]") for i in eachindex(distances)], Ref(names))

        Pde = Particles(MvNormal(Vector(c), Symmetric(Matrix(C) + 0.1*min(σL,σD)^2*I)))
        P = reshape(Pde[Pinds], dim, N)
        # de = Pde[dinds]
        (P=P, d=[norm(P[:,i] - P[:,j]) for i in 1:N, j in 1:N]), res
    end
    # chain
end

"""
    align_to_mean(P)

Takes an array of particles `P` and aligns all samples to the mean array using `procrustes`. This takes `P` from an estimate of absolute quantities to an estimate of relative quantities, which, in general, will exhibit lower variance.
"""
function align_to_mean(P)
    m = mean.(P)
    function apply_procrustes(X)
        R, t = procrustes(X, m)
        R * X .+ t
    end
    aligned_P = ℝⁿ2ℝⁿ_function(apply_procrustes, P)
end

end


text:
sparse_options = (; conic_solver = sparse_hypatia)

insts = OrderedDict()

insts["test"] = [
    ((1, 2, 2, true),),
    ((1, 2, 2, false),),
    ((2, 3, 3, true), sparse_options),
    ((2, 3, 3, false), sparse_options),
    ((6, 6, 3, true), sparse_options),
    ((6, 6, 3, false), sparse_options),
    ((20, 20, 3, true), sparse_options),
]

function polyfacilitylocation_insts(use_nat::Bool, n_max::Int)
    return [((n, 2n, 3, use_nat), sparse_options) for n in vcat(5, 5:5:n_max)]
end

insts["nat"] = polyfacilitylocation_insts(true, 50)
insts["ext"] = polyfacilitylocation_insts(false, 25)

return (PolyFacilityLocation, insts)


text:
include("../emi.jl")
using EMI
using EMI.Draw, EMI.Gmsh

canvas = Canvas()
for N in 3:10
    ngon = NGon(N, 1) + (N-3)*Point(2, 0)  # Fixed size shit
    canvas = canvas + ngon
end

set_bbox!(canvas, 0.2, 0.3)
gmsh_script(canvas, 0.2)

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
"""
    KNN(k=5,
        metric=Euclidean,
        algorithm=:kdtree,
        leafsize=10,
        reorder=true,
        reduction=:maximum)

Calculate the anomaly score of an instance based on the distance to its k-nearest neighbors.

Parameters
----------
$_k_param

$_knn_shared

    reduction::Symbol
One of `(:maximum, :median, :mean)`. (`reduction=:maximum`) was proposed by [1]. Angiulli et al. [2] proposed sum to
reduce the distances, but mean has been implemented for numerical stability.

Examples
--------
$(_score_unsupervised("KNN"))

References
----------
[1] Ramaswamy, Sridhar; Rastogi, Rajeev; Shim, Kyuseok (2000): Efficient Algorithms for Mining Outliers from Large Data
Sets.

[2] Angiulli, Fabrizio; Pizzuti, Clara (2002): Fast Outlier Detection in High Dimensional Spaces.
"""
MMI.@mlj_model mutable struct KNN <: UnsupervisedDetector
    k::Integer = 5::(_ > 0)
    metric::DI.Metric = DI.Euclidean()
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :balltree))
    leafsize::Integer = 10::(_ ≥ 0)
    reorder::Bool = true
    parallel::Bool = false
    reduction::Symbol = :maximum::(_ in (:maximum, :median, :mean))
end

struct KNNModel <: Model
    tree::NN.NNTree
end

function fit(detector::KNN, X::Data)::Fit
    # create the specified tree
    tree = buildTree(X, detector.metric, detector.algorithm, detector.leafsize, detector.reorder)

    # use tree to calculate distances
    _, dists = knn_others(tree, X, detector.k)

    # reduce distances to outlier score
    scores = _knn(dists, detector.reduction)
    Fit(KNNModel(tree), scores)
end

@score function score(detector::KNN, model::Fit, X::Data)::Result
    if detector.parallel
        idxs, dists = knn_parallel(model.tree, X, detector.k)
        return _knn(dists, detector.reduction)
    else
        idxs, dists = NN.knn(model.tree, X, detector.k)
        return _knn(dists, detector.reduction)
    end
end

@inline function _knn(distances::AbstractVector{<:AbstractVector}, reduction::Symbol)::Scores
    # Helper function to reduce `k` distances to a single distance.
    if reduction == :maximum
        return maximum.(distances)
    elseif reduction == :median
        return median.(distances)
    elseif reduction == :mean
        return mean.(distances)
    end
end


text:
module Heading

import Genie, Stipple

export heading

function heading(title::String = "";
                  class::String = "",
                  img::Union{String,Nothing} = nothing,
                  h = Genie.Renderer.Html.h1(title, class="st-header__title text-h3"),
                  kwargs...)

  Genie.Renderer.Html.header(; class="$(class) st-header q-pa-sm", kwargs...) do; [
    img
    h
  ]end
end

end