text:
### High level interface to CUDNN:
# 0. Init
# 1. Descriptors
# 2. Convolution
# 3. Pooling
# 4. LRNCrossChannel
# 5. DivisiveNormalization
# 6. Batch Normalization
# 7. Other tensor functions
# 8. Helper functions

### 0. INIT ##################################################

isdefined(Base, :__precompile__) && __precompile__()
module CUDNN
using Compat
using CUDArt

cudnnHandle = nothing

function __init__()
    # Setup default cudnn handle
    global cudnnHandle
    cudnnHandlePtr = cudnnHandle_t[0]
    cudnnCreate(cudnnHandlePtr)
    cudnnHandle = cudnnHandlePtr[1]
    # destroy cudnn handle at julia exit
    atexit(()->cudnnDestroy(cudnnHandle))
    global CUDNN_VERSION = convert(Int, ccall((:cudnnGetVersion,libcudnn),Csize_t,()))
end

const libcudnn = Libdl.find_library(["libcudnn"])
isempty(libcudnn) && error("CUDNN library cannot be found")

export CUDNN_VERSION
export cudnnTransformTensor, cudnnAddTensor, cudnnSetTensor, cudnnScaleTensor
export CUDNN_ADD_IMAGE, CUDNN_ADD_SAME_HW, CUDNN_ADD_FEATURE_MAP, CUDNN_ADD_SAME_CHW, CUDNN_ADD_SAME_C, CUDNN_ADD_FULL_TENSOR
export cudnnActivationForward, cudnnActivationBackward
export CUDNN_ACTIVATION_SIGMOID, CUDNN_ACTIVATION_RELU, CUDNN_ACTIVATION_TANH
export cudnnSoftmaxForward, cudnnSoftmaxBackward
export cudnnPoolingForward, cudnnPoolingBackward
export cudnnConvolutionForward
export CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION
export cudnnConvolutionBackwardBias, cudnnConvolutionBackwardFilter, cudnnConvolutionBackwardData
export CUDNN_POOLING_MAX, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING, CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING
export cudnnGetConvolutionNdForwardOutputDim, cudnnGetPoolingNdForwardOutputDim, cudnnGetConvolutionForwardWorkspaceSize
export cudnnLRNCrossChannelForward, cudnnLRNCrossChannelBackward, cudnnDivisiveNormalizationForward, cudnnDivisiveNormalizationBackward

import Base: unsafe_convert, conv2, strides
import CUDArt: free

# This is missing from CUDArt:
type cudaStream; end
typealias cudaStream_t Ptr{cudaStream}

# This is missing from cudnn.h:
const CUDNN_DIM_MAX = 8

# error handling function
function cudnnCheck(status)
    if status == CUDNN_STATUS_SUCCESS
        return nothing
    end
    # Because try/finally may disguise the source of the problem,
    # let's show a backtrace here
    warn("CUDNN error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(cudnnGetErrorString(status))
end

# These are semi-automatically generated using Clang with wrap_cudnn.jl:
include("types.jl")
include("libcudnn.jl")

# Conversion between Julia and CUDNN datatypes
cudnnDataType(::AbstractCudaArray{Float16})=CUDNN_DATA_HALF
cudnnDataType(::AbstractCudaArray{Float32})=CUDNN_DATA_FLOAT
cudnnDataType(::AbstractCudaArray{Float64})=CUDNN_DATA_DOUBLE
cudnnDataType(::Type{Float16})=CUDNN_DATA_HALF
cudnnDataType(::Type{Float32})=CUDNN_DATA_FLOAT
cudnnDataType(::Type{Float64})=CUDNN_DATA_DOUBLE
juliaDataType(a)=(a==CUDNN_DATA_HALF ? Float16 :
                  a==CUDNN_DATA_FLOAT ? Float32 :
                  a==CUDNN_DATA_DOUBLE ? Float64 : error())

### 1. DESCRIPTORS ##################################################

type TD; ptr; end
type FD; ptr; end
type CD; ptr; end
type PD; ptr; end
type LD; ptr; end
type AD; ptr; end

free(td::TD)=cudnnDestroyTensorDescriptor(td.ptr)
free(fd::FD)=cudnnDestroyFilterDescriptor(fd.ptr)
free(cd::CD)=cudnnDestroyConvolutionDescriptor(cd.ptr)
free(pd::PD)=cudnnDestroyPoolingDescriptor(pd.ptr)
free(ld::LD)=cudnnDestroyLRNDescriptor(ld.ptr)
free(ad::AD)=cudnnDestroyActivationDescriptor(ad.ptr)

unsafe_convert(::Type{cudnnTensorDescriptor_t}, td::TD)=td.ptr
unsafe_convert(::Type{cudnnFilterDescriptor_t}, fd::FD)=fd.ptr
unsafe_convert(::Type{cudnnConvolutionDescriptor_t}, cd::CD)=cd.ptr
unsafe_convert(::Type{cudnnPoolingDescriptor_t}, pd::PD)=pd.ptr
unsafe_convert(::Type{cudnnLRNDescriptor_t}, ld::LD)=ld.ptr
unsafe_convert(::Type{cudnnActivationDescriptor_t}, ad::AD)=ad.ptr

function TD(a::AbstractCudaArray, dims=ndims(a))
    (sz, st) = tensorsize(a, dims)
    d = cudnnTensorDescriptor_t[0]
    cudnnCreateTensorDescriptor(d)
    cudnnSetTensorNdDescriptor(d[1], cudnnDataType(a), length(sz), sz, st)
    this = TD(d[1])
    finalizer(this, free)
    return this
end

function FD(a::AbstractCudaArray, dims=ndims(a), format=CUDNN_TENSOR_NCHW)
    # The only difference of a FilterDescriptor is no strides.
    (sz, st) = tensorsize(a, dims)
    d = cudnnFilterDescriptor_t[0]
    cudnnCreateFilterDescriptor(d)
    CUDNN_VERSION >= 4000 ?
    cudnnSetFilterNdDescriptor_v4(d[1], cudnnDataType(a), format, length(sz), sz) :
    cudnnSetFilterNdDescriptor(d[1], cudnnDataType(a), length(sz), sz)
    this = FD(d[1])
    finalizer(this, free)
    return this
end

function CD(nd, padding, stride, upscale, mode, xtype)
    cd = cudnnConvolutionDescriptor_t[0]
    cudnnCreateConvolutionDescriptor(cd)
    CUDNN_VERSION >= 4000 ? cudnnSetConvolutionNdDescriptor(cd[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,cudnnDataType(xtype)) :
    CUDNN_VERSION >= 3000 ? cudnnSetConvolutionNdDescriptor_v3(cd[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,cudnnDataType(xtype)) :
    cudnnSetConvolutionNdDescriptor(cd[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode)
    this = CD(cd[1])
    finalizer(this, free)
    return this
end

CD(; ndims=2, padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION, eltype=Float32)=CD(ndims, padding, stride, upscale, mode, eltype)
CD(a::AbstractCudaArray; o...)=CD(; o..., ndims=ndims(a)-2, eltype=eltype(a))

function PD(nd, window, padding, stride, mode, maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    pd = cudnnPoolingDescriptor_t[0]
    cudnnCreatePoolingDescriptor(pd)
    CUDNN_VERSION >= 4000 ?
    cudnnSetPoolingNdDescriptor_v4(pd[1],mode,maxpoolingNanOpt,nd,pdsize(window,nd),pdsize(padding,nd),pdsize(stride,nd)) :
    cudnnSetPoolingNdDescriptor(pd[1],mode,nd,pdsize(window,nd),pdsize(padding,nd),pdsize(stride,nd))
    this = PD(pd[1])
    finalizer(this, free)
    return this
end

PD(; ndims=2, window=2, padding=0, stride=window, mode=CUDNN_POOLING_MAX, maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN)=PD(ndims,window,padding,stride,mode,maxpoolingNanOpt)
PD(a::AbstractCudaArray; o...)=PD(; o..., ndims=ndims(a)-2)

function LD(lrnN, lrnAlpha, lrnBeta, lrnK)
    ld = cudnnLRNDescriptor_t[0]
    cudnnCreateLRNDescriptor(ld)
    cudnnSetLRNDescriptor(ld[1], lrnN, lrnAlpha, lrnBeta, lrnK)
    this = LD(ld[1])
    finalizer(this, free)
    return this
end

# lrnN: Normalization window width in elements. LRN layer uses a window [center-lookBehind, center+lookAhead], where lookBehind = floor( (lrnN-1)/2 ), lookAhead = lrnN-lookBehind-1. So for lrnN=10, the window is [k-4...k...k+5] with a total of 10 samples. For DivisiveNormalization layer the window has the same extents as above in all 'spatial' dimensions (dimA[2], dimA[3], dimA[4]). By default lrnN is set to 5 in cudnnCreateLRNDescriptor.
# lrnAlpha: Value of the alpha variance scaling parameter in the normalization formula. Inside the library code this value is divided by the window width for LRN and by (window width)^#spatialDimensions for DivisiveNormalization. By default this value is set to 1e-4 in cudnnCreateLRNDescriptor.
# lrnBeta: Value of the beta power parameter in the normalization formula. By default this value is set to 0.75 in cudnnCreateLRNDescriptor.
# lrnK: Value of the k parameter in normalization formula. By default this value is set to 2.0.
LD(; lrnN=5, lrnAlpha=1e-4, lrnBeta=0.75, lrnK=2.0)=LD(lrnN, lrnAlpha, lrnBeta, lrnK)

# mode: Enumerant to specify the activation mode.
# reluNanOpt: Enumerant to specify the Nan propagation mode.
# reluCeiling: floating point number to specify the clipping threashod when the activation mode is set to CUDNN_ACTIVATION_CLIPPED_RELU.

function AD(mode, reluNanOpt, reluCeiling)
    ad = cudnnActivationDescriptor_t[0]
    cudnnCreateActivationDescriptor(ad)
    cudnnSetActivationDescriptor(ad[1], mode, reluNanOpt, reluCeiling)
    this = AD(ad[1])
    finalizer(this, free)
    return this
end

AD(; mode=CUDNN_ACTIVATION_RELU, reluNanOpt=CUDNN_NOT_PROPAGATE_NAN, reluCeiling=1.0)=AD(mode, reluNanOpt, reluCeiling)

# This is missing from CUDNN although mentioned in the documentation
const CUDNN_MAX_DIM = 8

function cudnnGetTensorNdDescriptor(td::TD)
    nd=CUDNN_MAX_DIM
    dataType = cudnnDataType_t[0]
    np = Array(Cint, 1)
    dimA = Array(Cint, nd)
    strideA = Array(Cint, nd)
    cudnnGetTensorNdDescriptor(td,nd,dataType,np,dimA,strideA)
    n = np[1]
    return (dataType[1], n, inttuple(dimA[1:n]), inttuple(strideA[1:n]))
end

function cudnnGetFilterNdDescriptor(fd::FD)
    nd=CUDNN_MAX_DIM
    dataType = cudnnDataType_t[0]
    np = Array(Cint, 1)
    dimA = Array(Cint, nd)
    format = cudnnTensorFormat_t[CUDNN_TENSOR_NCHW]
    CUDNN_VERSION >= 4000 ?
    cudnnGetFilterNdDescriptor_v4(fd,nd,dataType,format,np,dimA) :
    cudnnGetFilterNdDescriptor(fd,nd,dataType,np,dimA)
    n = np[1]
    return (dataType[1], format[1], n, inttuple(dimA[1:n]))
end

function cudnnGetConvolutionNdDescriptor(cd::CD)
    nd = CUDNN_DIM_MAX - 2
    np = Cint[0]
    p = Array(Cint, nd)
    s = Array(Cint, nd)
    u = Array(Cint, nd)
    m = cudnnConvolutionMode_t[0]
    t = cudnnDataType_t[0]
    CUDNN_VERSION >= 4000 ? cudnnGetConvolutionNdDescriptor(cd, nd, np, p, s, u, m, t) :
    CUDNN_VERSION >= 3000 ? cudnnGetConvolutionNdDescriptor_v3(cd, nd, np, p, s, u, m, t) :
    cudnnGetConvolutionNdDescriptor(cd, 2, np, p, s, u, m)
    n = np[1]
    (n, inttuple(p[1:n]), inttuple(s[1:n]), inttuple(u[1:n]), m[1], juliaDataType(t[1]))
end

function cudnnGetPoolingNdDescriptor(pd::PD)
    nd = CUDNN_MAX_DIM - 2
    np = Cint[0]
    m = cudnnPoolingMode_t[0]
    s = Array(Cint, nd)
    p = Array(Cint, nd)
    t = Array(Cint, nd)
    maxpoolingNanOpt = cudnnNanPropagation_t[CUDNN_NOT_PROPAGATE_NAN]
    CUDNN_VERSION >= 4000 ?
    cudnnGetPoolingNdDescriptor_v4(pd, nd, m, maxpoolingNanOpt, np, s, p, t) :
    cudnnGetPoolingNdDescriptor(pd, nd, m, np, s, p, t)
    n = np[1]
    (m[1], maxpoolingNanOpt[1], n, inttuple(s[1:n]), inttuple(p[1:n]), inttuple(t[1:n]))
end

function cudnnGetLRNDescriptor(ld::LD)
    lrnN = Cuint[0]; lrnAlpha=Cdouble[0]; lrnBeta=Cdouble[0]; lrnK=Cdouble[0]
    cudnnGetLRNDescriptor(ld, lrnN, lrnAlpha, lrnBeta, lrnK)
    (lrnN[1], lrnAlpha[1], lrnBeta[1], lrnK[1])
end

### 2. CONVOLUTION ##################################################

# CUDNN/Caffe sizes for various arrays in column-major notation:
# conv x: (W,H,C,N): W,H=image size, C=channels, N=instances
# conv w: (X,Y,C,K): X,Y=filter size, C=input channels, K=output channels
# conv y: (W-X+1,H-Y+1,K,N)
# conv b: (1,1,K,1)
# mmul x: (1,1,C,N): C=features, N=instances
# mmul w: (1,1,K,C)
# mmul y: (1,1,K,N)
# mmul b: (1,1,K,1)

# MATLAB equivalents:

# matlab 1-D y=conv(x,w,'full')  "Full convolution (default)."
# y[k] = sum[j] x[j] w[k-j+1]
# k=1:x+w-1
# j=max(1,k+1-w):min(k,x)
# cudnn: ndims=2, padding=w-1, stride=1, upscale=1, mode=CUDNN_CONVOLUTION
# Note: for cudnn you need to reshape the vectors to (n,1,1,1) and get
# the central (only non-zero) column of the result.

# matlab 1-D y=conv(x,w,'same')  "Central part of the convolution of the same size as u."
# y[k] = sum[j] x[j] w[k-j+w/2+1]
# k=1:x
# j=max(1,k-w/2+1):min(k+w/2,x)
# cudnn: ndims=2, padding=w/2, stride=1, upscale=1, mode=CUDNN_CONVOLUTION
# Note: This works only if length(w) is odd.  There is no setting to get 'same' with w even.

# matlab 1-D y=conv(x,w,'valid')  "Only those parts of the convolution that are computed without the zero-padded edges."
# y[k] = sum[j] x[j] w[k-j+w]
# k=1:x-w+1
# j=k:k+w-1
# cudnn: ndims=2, padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION

# TODO: use the new doc system.

# This function returns the dimensions of the resulting n-D tensor of a nbDims-2-D
# convolution, given the convolution descriptor, the input tensor descriptor and the filter
# descriptor This function can help to setup the output tensor and allocate the proper
# amount of memory prior to launch the actual convolution.
# Each dimension of the (nbDims-2)-D images of the output tensor is computed as
# followed:
#  outputDim = 1 + (inputDim + 2*pad - filterDim)/convolutionStride;

function cudnnGetConvolutionNdForwardOutputDim(src::AbstractCudaArray, filter::AbstractCudaArray; cd=nothing, o...)
    cd1 = (cd == nothing ? CD(src; o...) : cd)
    nbDims = ndims(src)
    outputDim = Array(Cint, nbDims)
    cudnnGetConvolutionNdForwardOutputDim(cd1, TD(src), FD(filter), nbDims, outputDim)
    # cd1 === cd || free(cd1)
    tuple(Int[reverse(outputDim)...]...)
end

# These are related to convolution algorithm selection:
# cudnnGetConvolutionForwardAlgorithm
# cudnnGetConvolutionForwardWorkspaceSize
#
# From the v2 release notes it seems like IMPLICIT_PRECOMP_GEMM is a
# good default especially with our memory limitations.  In v3 support for
# 3D is only available for IMPLICIT_GEMM and the speed difference does not
# seem significant, so I switched the default to IMPLICIT_GEMM.
#
# Forward convolution is now implemented via several different algorithms, and
# the interface allows the application to choose one of these algorithms specifically
# or to specify a strategy (e.g., prefer fastest, use no additional working space) by
# which the library should select the best algorithm automatically. The four
# algorithms currently given are as follows:
# o IMPLICIT_GEMM corresponds to the sole algorithm that was provided in
# cuDNN Release 1; it is the only algorithm that supports all input sizes while
# using no additional working space.
# o IMPLICIT_PRECOMP_GEMM is a modification of this approach that uses a
# small amount of working space (specifically, C*R*S*sizeof(int) bytes,
# where C is the number of input feature maps and R and S are the filter height
# and width, respectively, for the case of 2D convolutions) to achieve
# significantly higher performance in most cases. This algorithm achieves its
# highest performance when zero-padding is not used.
# o GEMM is an “im2col”-based approach, explicitly expanding the input data
# in memory and then using an otherwise-pure matrix multiplication that
# obeys cuDNN’s input and output stridings, avoiding a separate transposition
# step on the input or output. Note that this algorithm requires significant
# working space, though it may be faster than either of the two “implicit
# GEMM” approaches in some cases. As such, the PREFER_FASTEST
# forward convolution algorithm preference may sometimes recommend this
# approach. When memory should be used more sparingly,
# SPECIFY_WORKSPACE_LIMIT can be used instead of PREFER_FASTEST
# to ensure that the algorithm recommended will not require more than a given
# amount of working space.
# o DIRECT is a placeholder for a future implementation of direct convolution.

function cudnnGetConvolutionForwardAlgorithm(src::AbstractCudaArray, filter::AbstractCudaArray, dest::AbstractCudaArray;
                                             handle=cudnnHandle,
                                             preference=CUDNN_CONVOLUTION_FWD_PREFER_FASTEST, memoryLimitInbytes=10^9,
                                             cd=nothing, o...)
    cd1 = (cd == nothing ? CD(src; o...) : cd)
    algo = cudnnConvolutionFwdAlgo_t[0]
    cudnnGetConvolutionForwardAlgorithm(handle,TD(src),FD(filter),cd1,TD(dest),preference,memoryLimitInbytes,algo)
    # cd1 === cd || free(cd1)
    return algo[1]
end

function cudnnGetConvolutionForwardWorkspaceSize(src::AbstractCudaArray, filter::AbstractCudaArray, dest::AbstractCudaArray;
                                                 handle=cudnnHandle, algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,
                                                 cd=nothing, o...)
    cd1 = (cd == nothing ? CD(src; o...) : cd)
    sizeInBytes = Cint[0]
    cudnnGetConvolutionForwardWorkspaceSize(handle,TD(src),FD(filter),cd1,TD(dest),algorithm,sizeInBytes)
    # cd1 === cd || free(cd1)
    return Int(sizeInBytes[1])
end

function cudnnConvolutionForward(src, filter, dest=nothing;
                                 handle=cudnnHandle, alpha=1.0, beta=0.0,
                                 algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,
                                 workSpace=C_NULL, workSpaceSizeInBytes=0,
                                 cd=nothing, o...)
    # The user can specify the convolution descriptor by providing cd, or giving padding, stride, etc. in o...
    cd1 = (cd == nothing ? CD(src; o...) : cd)
    @assert eltype(filter) == eltype(src)
    osize = cudnnGetConvolutionNdForwardOutputDim(src,filter; cd=cd1)
    (dest == nothing) && (dest = CudaArray(eltype(src), osize))
    @assert osize == size(dest)
    @assert eltype(dest) == eltype(src)
    wsize = cudnnGetConvolutionForwardWorkspaceSize(src, filter, dest; algorithm=algorithm, cd=cd1)
    if ((wsize > 0) && (workSpace == C_NULL || workSpaceSizeInBytes < wsize))
        workSpaceSizeInBytes = wsize
        ws = CudaArray(Int8, workSpaceSizeInBytes)
    else
        ws = workSpace
    end
    cudnnConvolutionForward(handle,
                            cptr(alpha,src),TD(src),src,
                            FD(filter),filter,
                            cd1,algorithm,ws,workSpaceSizeInBytes,
                            cptr(beta,dest),TD(dest),dest)
    ws === workSpace || free(ws)
    # cd1 === cd || free(cd1)
    return dest
end

# conv2(x,w) in Julia performs 2-D convolution with padding equal to
# one less than the w dimensions.

Base.conv2{T}(src::AbstractCudaArray{T}, filter::AbstractCudaArray{T}, dest=nothing)=
    cudnnConvolutionForward(src, filter, dest; cd=CD(src; padding = (size(filter,1)-1, size(filter,2)-1)))

# n=h=w=1 for dest and c same as input.  CUDNN seems to assume a
# single scalar bias per output channel, i.e. the same number is added
# to every image and every pixel on the same output channel.  Say our
# input was x, we convolved it going forward with w, added bias b, and
# got tensor y as a result.  i.e. y=w*x+b where * means convolution
# and + means broadcasting addition. Now we get dJ/dy=src, and we want
# dJ/db=dest.  Well, we simply need to add up all dy where a
# particular b was added to get db.  This means dest is just the sum
# of src entries for every channel.
function cudnnConvolutionBackwardBias(src, dest=CudaArray(eltype(src),(1,1,size(src,3),1));
                                      handle=cudnnHandle, alpha=1.0, beta=0.0)
    cudnnConvolutionBackwardBias(handle,cptr(alpha,src),TD(src,4),src,cptr(beta,dest),TD(dest,4),dest)
    return dest
end

# I am guessing if y=w*x+b going forward, the arguments below
# correspond to src=x, diff=dy, grad=dw.
function cudnnConvolutionBackwardFilter(src, diff, grad;
                                        handle=cudnnHandle, alpha=1.0, beta=0.0,
                                        algo=CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0,
                                        workSpace=C_NULL, workSpaceSizeInBytes=0,
                                        cd=nothing, o...)
    cd1 = (cd == nothing ? CD(src; o...) : cd)
    CUDNN_VERSION >= 4000 ?
    cudnnConvolutionBackwardFilter(handle,cptr(alpha,src),TD(src),src,
                                   TD(diff),diff,cd1,
                                   algo, workSpace, workSpaceSizeInBytes,
                                   cptr(beta,grad),FD(grad),grad) :
    CUDNN_VERSION >= 3000 ?
    cudnnConvolutionBackwardFilter_v3(handle,
                                      cptr(alpha,src),TD(src),src,
                                      TD(diff),diff,cd1,
                                      algo, workSpace, workSpaceSizeInBytes,
                                      cptr(beta,grad),FD(grad),grad) :
    cudnnConvolutionBackwardFilter(handle,
                                   cptr(alpha,src),TD(src),src,
                                   TD(diff),diff,cd1,
                                   cptr(beta,grad),FD(grad),grad)
    return grad
end

# I am guessing if y=w*x+b going forward, the arguments below
# correspond to filter=w, diff=dy, grad=dx.
function cudnnConvolutionBackwardData(filter, diff, grad;
                                      handle=cudnnHandle, alpha=1.0, beta=0.0,
                                      algo=CUDNN_CONVOLUTION_BWD_DATA_ALGO_0,
                                      workSpace=C_NULL, workSpaceSizeInBytes=0,
                                      cd=nothing, o...)
    cd1 = (cd == nothing ? CD(filter; o...) : cd)
    CUDNN_VERSION >= 4000 ?
    cudnnConvolutionBackwardData(handle,cptr(alpha,diff),
                                 FD(filter),filter,
                                 TD(diff),diff,cd1,
                                 algo, workSpace, workSpaceSizeInBytes,
                                 cptr(beta,grad),TD(grad),grad) :
    CUDNN_VERSION >= 3000 ?
    cudnnConvolutionBackwardData_v3(handle,cptr(alpha,diff),
                                    FD(filter),filter,
                                    TD(diff),diff,cd1,
                                    algo, workSpace, workSpaceSizeInBytes,
                                    cptr(beta,grad),TD(grad),grad) :
    cudnnConvolutionBackwardData(handle,cptr(alpha,diff),
                                 FD(filter),filter,
                                 TD(diff),diff,cd1,
                                 cptr(beta,grad),TD(grad),grad)
    # cd1 === cd || free(cd1)
    return grad
end

### 3. POOLING ##################################################

# The cudnn documentation again does not tell anything about how the
# pooling parameters are used.  According to the caffe source
# (pooling_layer.cpp) here is what the parameters in the
# PoolingDescriptor mean: windowHeight and windowWidth give the size
# of the pooling area.  To compute the pooled output entry at position
# (ph,pw), we look at the input entries in the range:
#
# hstart = ph*verticalStride - verticalPadding
# wstart = pw*horizontalStride - horizontalPadding
# hend = hstart + windowHeight
# wend = wstart + windowWidth
#
# This is 0 based, the end points are exclusive, C style.
# The ranges are truncated if they are outside the input size.
# This makes the output size:
#
# pooledHeight = ceil((inputHeight + 2 * verticalPadding - windowHeight) / verticalStride) + 1
# pooledWidth = ceil((inputWidth + 2 * horizontalPadding - windowWidth) / horizontalStride) + 1
#
# If the output tensor is smaller, cudnn fills it with the upper left
# region of the pooled output.  If the output tensor is larger, it is
# filled up to the size of the input + 2 * padding, with the rest of
# the entries set to -inf (for max pooling) and NaN (for avg pooling).
# Actually there is a single row/col of NaNs and 0s outside.
#
# CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING is buggy when padding >
# 0 or else I have no idea what it is doing.

function cudnnPoolingForward(src, dest;
                             handle=cudnnHandle, alpha=1.0, beta=0.0,
                             pd = nothing, o...)
    pd1 = (pd == nothing ? PD(src; o...) : pd)
    cudnnPoolingForward(handle, pd1,
                        cptr(alpha,src), TD(src), src,
                        cptr(beta,dest), TD(dest), dest)
    # pd1 === pd || free(pd1)
    return dest
end

function cudnnPoolingBackward(src, srcDiff, dest, destDiff;
                              handle=cudnnHandle, alpha=1.0, beta=0.0,
                              pd = nothing, o...)
    pd1 = (pd == nothing ? PD(src; o...) : pd)
    cudnnPoolingBackward(handle, pd1,
                         cptr(alpha,src), TD(src), src,
                         TD(srcDiff), srcDiff,
                         TD(dest), dest,
                         cptr(beta,destDiff), TD(destDiff), destDiff)
    # pd1 === pd || free(pd1)
    return destDiff
end

# v2: This does not seem to exist in libcudnn.so even though it is declared in cudnn.h and the docs
# It is also in the static library but not the .so!  Weird.
# v3 update: now it is in the library, but it is buggy, just gives the dimensions of src back.
# v4 update: still buggy, I am amazed.
function cudnnGetPoolingNdForwardOutputDim_buggy(pd::PD, src::AbstractCudaArray)
    if CUDNN_VERSION >= 3000
        nbDims = ndims(src)
        outputTensorDimA = Array(Cint, nbDims)
        cudnnGetPoolingNdForwardOutputDim(pd.ptr, TD(src), nbDims, outputTensorDimA)
        tuple(Int[reverse(outputTensorDimA)...]...)
    end
end

# Here is the corrected version based on the specs
function cudnnGetPoolingNdForwardOutputDim(pd::PD, input::AbstractCudaArray)
    dims = [size(input)...]
    (mode, mpno, pdims, window, padding, stride) = cudnnGetPoolingNdDescriptor(pd)
    for i=1:length(dims)-2
        dims[i] = 1 + ceil((dims[i] + 2*padding[i] - window[i]) / stride[i])
    end
    tuple(dims...)
end

### 4. LRNCrossChannel

function cudnnLRNCrossChannelForward(src, dest;
                                     handle=cudnnHandle,
                                     mode=CUDNN_LRN_CROSS_CHANNEL_DIM1,
                                     alpha=1.0, beta=0.0,
                                     ld=nothing, o...)
    ld1 = (ld != nothing ? ld : LD(; o...))
    cudnnLRNCrossChannelForward(handle, ld1, mode,
                                cptr(alpha,src), TD(src), src,
                                cptr(beta,dest), TD(dest), dest)
    return dest
end

function cudnnLRNCrossChannelBackward(src, srcDiff, dest, destDiff;
                                      handle=cudnnHandle,
                                      mode=CUDNN_LRN_CROSS_CHANNEL_DIM1,
                                      alpha=1.0, beta=0.0,
                                      ld=nothing, o...)
    ld1 = (ld != nothing ? ld : LD(; o...))
    cudnnLRNCrossChannelBackward(handle, ld1, mode,
                                 cptr(alpha,src), TD(src), src,
                                 TD(srcDiff), srcDiff,
                                 TD(dest), dest,
                                 cptr(beta,destDiff), TD(destDiff), destDiff)
    return destDiff
end

### 5. DIVISIVE NORMALIZATION ##################################################

function cudnnDivisiveNormalizationForward(src, dest;
                                           means=nothing, temp1=nothing, temp2=nothing,
                                           handle=cudnnHandle,
                                           mode=CUDNN_DIVNORM_PRECOMPUTED_MEANS,
                                           alpha=1.0, beta=0.0,
                                           ld=nothing, o...)
    ld1 = (ld != nothing ? ld : LD(; o...))
    means == nothing && (means = C_NULL)
    p1 = (temp1 == nothing ? similar(src) : temp1)
    p2 = (temp2 == nothing ? similar(src) : temp2)
    cudnnDivisiveNormalizationForward(handle, ld1, mode,
                                      cptr(alpha,src), TD(src), src,
                                      means, p1, p2,
                                      cptr(beta,dest), TD(dest), dest)
    p1 === temp1 || free(p1)
    p2 === temp2 || free(p2)
    return dest
end

function cudnnDivisiveNormalizationBackward(src, srcDiff, destDiff;
                                            srcMeans=nothing, destMeansDiff=nothing,
                                            temp1=nothing, temp2=nothing,
                                            handle=cudnnHandle,
                                            mode=CUDNN_DIVNORM_PRECOMPUTED_MEANS,
                                            alpha=1.0, beta=0.0,
                                            ld=nothing, o...)
    ld1 = (ld != nothing ? ld : LD(; o...))
    srcMeans == nothing && (srcMeans = C_NULL)
    destMeansDiff == nothing && (destMeansDiff = C_NULL)
    p1 = (temp1 == nothing ? similar(src) : temp1)
    p2 = (temp2 == nothing ? similar(src) : temp2)
    cudnnDivisiveNormalizationBackward(handle, ld1, mode,
                                       cptr(alpha,src), TD(src), src,
                                       srcMeans, srcDiff, p1, p2,
                                       cptr(beta,destDiff), TD(destDiff), destDiff, destMeansDiff)
    p1 === temp1 || free(p1)
    p2 === temp2 || free(p2)
    return destDiff
end

### 6. BATCH NORMALIZATION ##################################################

# Modes:
# CUDNN_BATCHNORM_PER_ACTIVATION: Normalization is performed per-activation. This mode is intended to be used after non- convolutional network layers. In this mode bnBias and bnScale tensor dimensions are 1xCxHxW.
# CUDNN_BATCHNORM_SPATIAL: Normalization is performed over N+spatial dimensions. This mode is intended for use after convolutional layers (where spatial invariance is desired). In this mode bnBias, bnScale tensor dimensions are 1xCx1x1.

# This function performs the forward BatchNormalization layer
# computation for inference phase. This layer is based on the paper
# "Batch Normalization: Accelerating Deep Network Training by Reducing
# Internal Covariate Shift", S. Ioffe, C. Szegedy, 2015.  Only 4D and 5D
# tensors are supported.  The epsilon value has to be the same during
# training, backpropagation and inference.  For training phase use
# cudnnBatchNormalizationForwardTraining.  Much higher performance when
# HW-packed tensors are used for all of x, dy, dx.

# TODO: read the paper and cudnn doc to figure out a good interface.

### 7. OTHER TENSOR FUNCTIONS ##################################################

# alpha * src + beta * dest -> dest
# Both beta and dest optional, beta=0 if not specified, dest is allocated with ones if not specified.
# These defaults give the expected answers for (alpha * src) and (alpha * src + beta)
# The doc says no in-place, i.e. src and dest should not overlap.
# My experiments say otherwise but I will go with the doc just in case.

function cudnnTransformTensor(alpha::Number, src::AbstractCudaArray, beta::Number=0, dest::AbstractCudaArray=ones(src);
                              handle=cudnnHandle)
    cudnnTransformTensor(handle,
                         cptr(alpha,src), TD(src,4), src,
                         cptr(beta,dest), TD(dest,4), dest)
    return dest
end

# Refer to cudnn doc to see what different add modes do

function cudnnAddTensor(bias::AbstractCudaArray, src::AbstractCudaArray;
                        handle=cudnnHandle, alpha=1.0, beta=1.0, mode=CUDNN_ADD_SAME_C)
    CUDNN_VERSION >= 4000 ? cudnnAddTensor(handle, cptr(alpha,bias), TD(bias,4), bias, cptr(beta,src), TD(src,4), src) :
    CUDNN_VERSION >= 3000 ? cudnnAddTensor_v3(handle, cptr(alpha,bias), TD(bias,4), bias, cptr(beta,src), TD(src,4), src) :
    CUDNN_VERSION >= 2000 ? cudnnAddTensor(handle, mode, cptr(alpha,bias), TD(bias,4), bias, cptr(beta,src), TD(src,4), src) :
    error("Unsupported version $CUDNN_VERSION")
    return src
end

# src .= value

function cudnnSetTensor(src::AbstractCudaArray, value::Number; handle=cudnnHandle)
    cudnnSetTensor(handle, TD(src,4), src, cptr(value,src))
    return src
end

# src .*= alpha

function cudnnScaleTensor(src::AbstractCudaArray, alpha::Number; handle=cudnnHandle)
    cudnnScaleTensor(handle, TD(src,4), src, cptr(alpha,src))
    return src
end

# Apply activation fn to each element of src (dest optional, in-place by default)
# mode is one of CUDNN_ACTIVATION_{SIGMOID,RELU,TANH}

function cudnnActivationForward(src::AbstractCudaArray, dest::AbstractCudaArray=src; handle=cudnnHandle,
                                mode=CUDNN_ACTIVATION_RELU, reluNanOpt=CUDNN_NOT_PROPAGATE_NAN, reluCeiling=1.0,
                                alpha=1.0, beta=0.0)
    if CUDNN_VERSION >= 4000
        cudnnActivationForward_v4(handle, AD(mode, reluNanOpt, reluCeiling),
                                  cptr(alpha,src), TD(src,4), src,
                                  cptr(beta,dest), TD(dest,4), dest)
    else
        cudnnActivationForward(handle, mode,
                               cptr(alpha,src), TD(src,4), src,
                               cptr(beta,dest), TD(dest,4), dest)
    end
    return dest
end

# Compute activation fn gradient.  The naming is a bit confusing.  In
# KUnet terminology, if we have input x and output y going forward, we
# get dy, the gradient wrt output y, and compute dx, the gradient wrt
# input x (dx overwrites dy by default).  In cudnn, x=dest, y=src,
# dx=destDiff, dy=srcDiff.  The gradient calculation should not need
# x=dest.  cudnn seems to set dx=0 whereever x=0.

function cudnnActivationBackward(src::AbstractCudaArray, srcDiff::AbstractCudaArray, dest::AbstractCudaArray=src, destDiff::AbstractCudaArray=srcDiff;
                                 mode=CUDNN_ACTIVATION_RELU, reluNanOpt=CUDNN_NOT_PROPAGATE_NAN, reluCeiling=1.0,
                                 handle=cudnnHandle, alpha=1.0, beta=0.0)
    if CUDNN_VERSION >= 4000
        cudnnActivationBackward_v4(handle, AD(mode, reluNanOpt, reluCeiling),
                                   cptr(alpha,src), TD(src,4), src,
                                   TD(srcDiff,4), srcDiff,
                                   TD(dest,4), dest,
                                   cptr(beta,destDiff), TD(destDiff,4), destDiff)
    else
        cudnnActivationBackward(handle, mode,
                                cptr(alpha,src), TD(src,4), src,
                                TD(srcDiff,4), srcDiff,
                                TD(dest,4), dest,
                                cptr(beta,destDiff), TD(destDiff,4), destDiff)
    end
    return destDiff
end

# In KUnet I implement softmax as a normalization function on the
# final output which is typically a small number of units representing
# classes.  CUDNN softmax acts on Tensors?  What does that mean?  Doc says:
# mode=INSTANCE: The softmax operation is computed per image (N) across the dimensions C,H,W.
# mode=CHANNEL: The softmax operation is computed per spatial location (H,W) per image (N) across the dimension C.
# From caffe docs: image data blob has NCHW dims.
# Regular feature vectors (that are input to InnerProduct layers) have NC11 dims. (N vectors holding C dims each)
# Convolution parameter blob has NCHW where N output filters, C input channels, HxW images.
# Inner product weight matrix has 11HW (or is it CK11) for H output and W input dims.
# Caffe softmax takes pred and label blobs and computes a single number.
# cudnnSoftmaxForward seems to just exponentiate and normalize.
# This is more like an activation function rather than a loss function.
# It takes a NC11 tensor (N vectors, holding C unnormalized logp each) and returns an NC11 tensor with normalized p.
# Note that an NC11 tensor is constructed from a julia array of size (1,1,C,N).
# Can we do in-place?  It turns out yes (undocumented).

function cudnnSoftmaxForward(src::AbstractCudaArray, dest::AbstractCudaArray=src;
                             handle=cudnnHandle,
                             algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                             mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                             alpha=1.0, beta=0.0)
    cudnnSoftmaxForward(handle, algorithm, mode,
                        cptr(alpha, src), TD(src,4), src,
                        cptr(beta, dest), TD(dest,4), dest)
    return dest
end

# If SoftmaxForward is implemented as an activation fn rather than a
# loss fn, it stands to reason that its backward partner is similar.
# If the input to the forward pass was x, array of unnormalized logp,
# and forward output was y, array of normalized p, then I am going to
# guess below that src=y, srcDiff=dy, and destDiff=dx.  Nothing in
# docs about in-place ops but I will assume it is possible as the
# structure is similar to cudnnActivateBackward (except dest=x is not in
# the arglist).
#
# But what does it compute?  Given a one-of-k vector of correct
# answers ai we have: dJ/dxi = (yi - ai), so what is srcDiff?  Do they
# really mean dJ/dyi?  Who computes dJ/dyi?  What idiot wrote this
# documentation without a single mathematical expression?
#
# OK, if src is normalized probabilities output by the model, and we
# assume y1 is the probability of the correct answer we have
# J=-log(y1) and we should specify srcDiff as dJ/dy1 = -1/y1 and
# dJ/dyi=1/y1 (for i!=1).  If we do, we should get destDiff
# dJ/dx1=y1-1 and dJ/dxi=yi (for i!=1).  We get twice those answers.
# I have no idea where the factor of 2 comes from.  I recommend
# dividing dJ/dyi with 2 for srcDiff to get the correct derivatives.
# You should also divide all dJ/dyi by the number of instances (N) to
# make learningRate specify the same step size regardless of batch
# size.  (Or maybe one should divide the loss J instead).

function cudnnSoftmaxBackward(src::AbstractCudaArray, srcDiff::AbstractCudaArray, destDiff::AbstractCudaArray=srcDiff;
                              handle=cudnnHandle,
                              algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                              mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                              alpha=1.0, beta=0.0)
    cudnnSoftmaxBackward(handle, algorithm, mode,
                         cptr(alpha, src), TD(src,4), src,
                         TD(srcDiff,4), srcDiff,
                         cptr(beta, destDiff), TD(destDiff,4), destDiff)
    return destDiff
end

### 8. HELPER FUNCTIONS ##################################################

# This is missing from CUDArt:
Base.strides(a::AbstractCudaArray)=map(i->stride(a,i), tuple(1:ndims(a)...))

# For low level cudnn functions that require a pointer to a number
cptr(x,a::CudaArray{Float64})=Float64[x]
cptr(x,a::CudaArray{Float32})=Float32[x]
cptr(x,a::CudaArray{Float16})=Float32[x]

inttuple(x)=tuple(Int[x...]...)

# size and strides in Julia are in column-major format, e.g. (W,H,C,N)
# whereas cudnn uses row-major, e.g. NCHW
# most cudnn ops still work only on 4-D arrays
# so we may need to modify size and strides

function tensorsize(a, dims)
    sz = Cint[reverse(size(a))...]
    st = Cint[reverse(strides(a))...]
    if length(sz) == 1 < dims
        unshift!(sz, 1)
        unshift!(st, 1)
    end
    while length(sz) < dims
        push!(sz, 1)
        push!(st, 1)
    end
    while length(sz) > dims
        d = pop!(sz)
        sz[length(sz)] *= d
        pop!(st)
        st[length(st)] = 1
    end
    (sz, st)
end

function cdsize(w, nd)
    isa(w,Integer) ? Cint[fill(w,nd)...] :
    length(w)!=nd ? error("Dimension mismatch") :
    Cint[reverse(w)...]
end

pdsize(w, nd)=Cint[reverse(psize(w,nd))...]
psize(w, nd)=(isa(w,Integer)  ? fill(w,nd) : length(w) != nd ? error("Dimension mismatch") : w)

end # module

### DEAD CODE:

# I am no longer keeping default conv or pool descriptors.  These
# descriptors are created and destroyed by the operations that need
# them.  Thus the following functions have been modified:

# function cudnnConvolutionForward(src::AbstractCudaArray, filter::AbstractCudaArray, dest=nothing;
#                                  handle=cudnnHandle, alpha=1.0, beta=0.0,
#                                  convDesc=defaultConvolutionDescriptor,
#                                  algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM,
#                                  workSpace=C_NULL, workSpaceSizeInBytes=0)
#     @assert eltype(filter) == eltype(src)
#     osize = cudnnGetConvolutionNdForwardOutputDim(src,filter;convDesc=convDesc)
#     (dest == nothing) && (dest = CudaArray(eltype(src), osize))
#     @assert osize == size(dest)
#     @assert eltype(dest) == eltype(src)
#     wsize = cudnnGetConvolutionForwardWorkspaceSize(src, filter, dest; algorithm=algorithm)
#     if ((wsize > 0) && (workSpace == C_NULL || workSpaceSizeInBytes < wsize))
#         workSpaceSizeInBytes = wsize
#         workSpace = CudaArray(Int8, workSpaceSizeInBytes)
#     end
#     cudnnConvolutionForward(handle,
#                             cptr(alpha,src),TD(src),src,
#                             FD(filter),filter,
#                             convDesc,algorithm,workSpace,workSpaceSizeInBytes,
#                             cptr(beta,dest),TD(dest),dest)
#     return dest
# end

# type ConvolutionDescriptor; ptr; padding; stride; upscale; mode; datatype; end

# # TODO: accept single number args for padding etc.
# function ConvolutionDescriptor(; padding=(0,0),
#                                stride=map(x->1,padding),
#                                upscale=map(x->1,padding),
#                                mode=CUDNN_CONVOLUTION,
#                                datatype=Float64)
#     @assert in(mode, (CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION))
#     # @assert length(padding) == length(stride) == length(upscale)
#     cd = cudnnConvolutionDescriptor_t[0]
#     cudnnCreateConvolutionDescriptor(cd)
#     cudnnSetConvolutionNdDescriptor(cd[1],
#                                     length(padding),
#                                     Cint[reverse(padding)...],
#                                     Cint[reverse(stride)...],
#                                     Cint[reverse(upscale)...],
#                                     mode, cudnnDataType(datatype))
#     this = ConvolutionDescriptor(cd[1], padding, stride, upscale, mode, datatype)
#     finalizer(this, free)
#     this
# end

# free(cd::ConvolutionDescriptor)=cudnnDestroyConvolutionDescriptor(cd.ptr)
# unsafe_convert(::Type{cudnnConvolutionDescriptor_t}, cd::ConvolutionDescriptor)=cd.ptr

# const defaultConvolutionDescriptor = ConvolutionDescriptor()

# type PoolingDescriptor; ptr; dims; padding; stride; mode; end

# # TODO: allow initialization with a single number

# function PoolingDescriptor(dims; padding=map(x->0,dims), stride=dims, mode=CUDNN_POOLING_MAX)
#     @assert in(mode, (CUDNN_POOLING_MAX,
#                       CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING,
#                       CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING))
#     @assert length(dims) == length(padding) == length(stride)
#     pd = cudnnPoolingDescriptor_t[0]
#     cudnnCreatePoolingDescriptor(pd)
#     cudnnSetPoolingNdDescriptor(pd[1],mode,length(dims),Cint[reverse(dims)...],Cint[reverse(padding)...],Cint[reverse(stride)...])
#     this = PoolingDescriptor(pd[1], dims, padding, stride, mode)
#     finalizer(this, free)
#     this
# end

# free(pd::PoolingDescriptor)=cudnnDestroyPoolingDescriptor(pd.ptr)
# unsafe_convert(::Type{cudnnPoolingDescriptor_t}, pd::PoolingDescriptor)=pd.ptr

# function cudnnPoolingForward(pd::PoolingDescriptor, src::AbstractCudaArray, dest::AbstractCudaArray;
#                              handle=cudnnHandle, alpha=1.0, beta=0.0)
#     cudnnPoolingForward(handle, pd,
#                         cptr(alpha,src), TD(src), src,
#                         cptr(beta,dest), TD(dest), dest)
#     return dest
# end

# # Read info from gpu for debugging
# function cudnnGetPoolingNdDescriptor(pd::PoolingDescriptor)
#     nd = length(pd.dims)
#     m = cudnnPoolingMode_t[0]
#     n = Cint[0]
#     s = Array(Cint, nd)
#     p = Array(Cint, nd)
#     t = Array(Cint, nd)
#     cudnnGetPoolingNdDescriptor(pd, nd, m, n, s, p, t)
#     inttuple(x)=tuple(Int[x...]...)
#     (m[1], n[1], inttuple(s), inttuple(p), inttuple(t))
# end

# type TD; ptr; end
# type FD; ptr; end
# free(td::TD)=cudnnDestroyTensorDescriptor(td.ptr)
# free(fd::FD)=cudnnDestroyFilterDescriptor(fd.ptr)
# unsafe_convert(::Type{cudnnTensorDescriptor_t}, td::TD)=td.ptr
# unsafe_convert(::Type{cudnnFilterDescriptor_t}, fd::FD)=fd.ptr


text:
using ModelingToolkit, GalacticOptim, Optim

@variables t x(t) v(t) u(t)
@parameters p[1:2]
D = Differential(t)

loss = (4-x)^2 + 2v^2 + u^2
eqs = [
    D(x) ~ v - p[2]*x
    D(v) ~ p[1]*u^3 + v
]

sys = ControlSystem(loss,eqs,t,[x,v],[u],p)
dt = 0.1
tspan = (0.0,1.0)
sys = runge_kutta_discretize(sys,dt,tspan)

u0 = rand(length(states(sys))) # guess for the state values
prob = OptimizationProblem(sys,u0,[0.1,0.1],grad=true)
sol = solve(prob,BFGS())

# issue #819
@testset "Combined system name collisions" begin
    eqs_short = [
        D(x) ~ - p[2]*x
        D(v) ~ p[1]*u^3
    ]
    sys1 = ControlSystem(loss,eqs_short, t, [x, v], [u], p, name = :sys1)
    sys2 = ControlSystem(loss,eqs_short, t, [x, v], [u], p, name = :sys1)
    @test_throws ArgumentError ControlSystem(loss, [sys2.v ~ sys1.v], t, [], [], [], systems = [sys1, sys2])
end

ext: jl


text:
module Factories

using MacroTools: block
using MacroTools: splitarg,combinefield,combinearg
using MacroTools: splitstructdef,combinestructdef
using MacroTools: splitdef,combinedef
using Printf: @printf,@sprintf
using ..TypeTraits: efficientoperations

import MacroTools: rmlines

export FExpr,escape
export Escaped,UnEscaped,MixEscaped
export RawExpr,rawexpr
export Inference,@inference
export Argument,@argument
export Parameter,@parameter
export Field,@field
export Block,@block,@push!,rmlines!,@rmlines!,rmlines,@rmlines
export FunctionFactory,@functionfactory,addargs!,@addargs!,addkwargs!,@addkwargs!,addwhereparams!,@addwhereparams!,extendbody!,@extendbody!
export TypeFactory,@typefactory,addfields!,@addfields!,addconstructors!,@addconstructors!
export addparams!,@addparams!

@generated function maxfieldnamelength(::T) where T
    result=1
    for name in T|>fieldnames
        result=max(result,name|>string|>length)
    end
    return :($result)
end

"Factory expression types, which is defined as `Union{Symbol,Expr}`."
const FExpr=Union{Symbol,Expr}

"Abstract escape mechanism."
abstract type EscapeMechanism end

"""
    Escaped(names::Symbol...)

Indicate that symbols of a factory should be escaped if they are in `names`.
"""
struct Escaped{N} <: EscapeMechanism
    names::NTuple{N,Symbol}
    Escaped(names::Symbol...)=new{names|>length}(names)
end

"""
    UnEscaped(names::Symbol...)

IIndicate that symbols of a factory should be escaped if they are not in `names`.
"""
struct UnEscaped{N} <: EscapeMechanism
    names::NTuple{N,Symbol}
    UnEscaped(names::Symbol...)=new{names|>length}(names)
end

"""
    MixEscaped(escaped::Escaped)
    MixEscaped(unescaped::UnEscaped)
    MixEscaped(escaped::Escaped,unescaped::UnEscaped)
    MixEscaped(unescaped::UnEscaped,escaped::Escaped)

Indicate that some parts of a factory use the `Escaped` mechanism while other parts use the `UnEscaped` mechanism.
"""
struct MixEscaped{N,M} <: EscapeMechanism
    escaped::Escaped{N}
    unescaped::UnEscaped{M}
end
MixEscaped()=MixEscaped(Escaped(),UnEscaped())
MixEscaped(escaped::Escaped)=MixEscaped(escaped,UnEscaped())
MixEscaped(unescaped::UnEscaped)=MixEscaped(Escaped(),unescaped)
MixEscaped(unescaped::UnEscaped,escaped::Escaped)=MixEscaped(escaped,unescaped)

"Raw expression without any variable escaped."
struct RawExpr <: EscapeMechanism end
"""
    rawexpr

Indicate that no variable in a factory should be escaped.
"""
const rawexpr=RawExpr()

"""
    escape(expr,::RawExpr) -> Any
    escape(expr,::Escaped) -> Any
    escape(expr,::UnEscaped) -> Any
    escape(expr::Symbol,em::Escaped) -> FExpr
    escape(expr::Expr,em::Escaped) -> Expr
    escape(expr::Symbol,em::UnEscaped) -> FExpr
    escape(expr::Expr,em::UnEscaped) -> Expr

Escape the variables in the input expression.
"""
escape(expr,::RawExpr)=expr
escape(expr,::Escaped)=expr
escape(expr,::UnEscaped)=expr
escape(expr::Symbol,em::Escaped)=expr ∈ em.names ? :($(esc(expr))) : expr
escape(expr::Expr,em::Escaped)=(expr.head==:escape ? expr : Expr(expr.head,(escape(arg,em) for arg in expr.args)...))
escape(expr::Symbol,em::UnEscaped)=expr ∉ em.names ? :($(esc(expr))) : expr
function escape(expr::Expr,em::UnEscaped)
    if expr.head==:escape
        @assert length(expr.args)==1 "escape error: internal error..."
        escape(expr.args[1],em)
    else
        Expr(expr.head,(escape(arg,em) for arg in expr.args)...)
    end
end

"""
    names(expr::Symbol) -> Vector{Symbol}
    names(expr::Expr) -> Vector{Symbol}

Get all the symbols in an expression.
"""
Base.names(expr::Symbol)=[expr]
function Base.names(expr::Expr)
    result=Symbol[]
    for arg in expr.args
        (isa(arg,Symbol)|| isa(arg,Expr)) && append!(result,names(arg))
    end
    return result
end

"""
    AbstractFactory

Abstract type for all concrete factories.
"""
abstract type AbstractFactory end

"""
    ==(f1::F,f2::F) where F<:AbstractFactory -> Bool
    isequal(f1::F,f2::F) where F<:AbstractFactory -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(f1::F,f2::F) where F<:AbstractFactory = ==(efficientoperations,f1,f2)
Base.isequal(f1::F,f2::F) where F<:AbstractFactory=isequal(efficientoperations,f1,f2)

"""
    show(io::IO,f::AbstractFactory)

Show a concrete `AbstractFactory`.
"""
function Base.show(io::IO,f::AbstractFactory)
    @printf io "%s(\n" f|>typeof|>nameof
    for name in f|>typeof|>fieldnames
        nstr=" "^2*string(name)*":"*" "^((f|>maxfieldnamelength)+1-(name|>string|>length))
        value=getfield(f,name)
        if isa(value,Nothing)
            vstr="nothing"
        elseif isa(value,AbstractFactory)
            vstr=string(value(rawexpr))
        elseif isa(value,Vector{<:AbstractFactory})
            eltypename=value|>eltype|>nameof
            strings=Tuple(string(af(rawexpr)) for af in value)
            if (length(strings)>0 ? sum(map(length,strings)) : 0)<90-(f|>maxfieldnamelength)
                vstr=@sprintf "%s[%s]" eltypename join(strings,", ")
            else
                vstr=replace((@sprintf "%s[\n%s\n]" eltypename join(strings,"\n")),"\n"=>"\n  ")
            end
        else
            vstr=string(value)
        end
        '\n' in vstr && (vstr=replace(vstr,"\n"=>"\n"*" "^((f|>maxfieldnamelength)+4)))
        @printf(io,"%s%s\n",nstr,vstr)
    end
    @printf io ")"
end

"""
    replace(f::AbstractFactory;kwargs...) -> typeof(f)

Return a copy of a concrete `AbstractFactory` with some of the field values replaced by the keyword arguments.
"""
Base.replace(f::AbstractFactory;kwargs...)=replace(efficientoperations,f;kwargs...)

"""
    Inference(head::Union{Symbol,Nothing},name::Union{Symbol,Nothing},params::Union{Inference,Vector{Inference},Nothing})
    Inference(;
            head::Union{Symbol,Nothing}=nothing,
            name::Union{Symbol,Nothing}=nothing,
            params::Union{Inference,Vector{Inference},Nothing}=nothing,
            )
    Inference(expr::FExpr)

The struct to describe a type inference.
"""
mutable struct Inference <: AbstractFactory
    head::Union{Symbol,Nothing}
    name::Union{Symbol,Nothing}
    params::Union{Inference,Vector{Inference},Nothing}
    function Inference(head::Union{Symbol,Nothing},name::Union{Symbol,Nothing},params::Union{Inference,Vector{Inference},Nothing})
        @assert (head===nothing)==(params===nothing) "Inference error: `head` and `params` must be or not be `nothing` at the same time."
        new(head,name,params)
    end
end
Inference(;
        head::Union{Symbol,Nothing}=nothing,
        name::Union{Symbol,Nothing}=nothing,
        params::Union{Inference,Vector{Inference},Nothing}=nothing
        )=Inference(head,name,params)
function Inference(expr::FExpr)
    if isa(expr,Symbol)
        return Inference(nothing,expr,nothing)
    elseif expr.head==:(<:)
        @assert length(expr.args)==1 "Inference error: wrong input expr."
        return Inference(expr.head,nothing,Inference(expr.args[1]))
    elseif expr.head==:curly
        return Inference(expr.head,expr.args[1],Inference.(@views(expr.args[2:end])))
    else
        error("Inference error: wrong input expr.")
    end
end

"""
    @inference expr::FExpr

Construct an `Inference` directly from a type inference.
"""
macro inference(expr::FExpr) expr=[expr];:(Inference($expr...)) end

"""
    (i::Inference)(em::RawExpr) -> FExpr
    (i::Inference)(em::UnEscaped) -> FExpr
    (i::Inference)(em::MixEscaped) -> FExpr

Convert a `Inference` to the `Expr` representation of the type inference it describes.
"""
(i::Inference)(em::MixEscaped)=i(em.unescaped)
function (i::Inference)(em::Union{RawExpr,<:UnEscaped})
    if i.head===nothing
        return escape(i.name,em)
    elseif i.head==:(<:)
        return Expr(i.head,i.params(em))
    else
        return Expr(i.head,escape(i.name,em),[param(em) for param in i.params]...)
    end
end

"""
    Argument(name::Union{Symbol,Nothing},type::Inference,slurp::Bool,default::Any)
    Argument(;name::Union{Symbol,Nothing}=nothing,type::Inference=Inference(:Any),slurp::Bool=false,default::Any=nothing)
    Argument(expr::FExpr)

The struct to describe a argument of a `function`.
"""
mutable struct Argument <: AbstractFactory
    name::Union{Symbol,Nothing}
    type::Inference
    slurp::Bool
    default::Any
end
Argument(;name::Union{Symbol,Nothing}=nothing,type::Inference=Inference(:Any),slurp::Bool=false,default::Any=nothing)=Argument(name,type,slurp,default)
Argument(expr::FExpr)=(cache=splitarg(expr);Argument(cache[1],Inference(cache[2]),cache[3],cache[4]))

"""
    @argument expr::FExpr

Construct an `Argument` directly from an argument statement.
"""
macro argument(expr::FExpr) expr=[expr];:(Argument($expr...)) end

"""
    (a::Argument)(em::RawExpr) -> Expr
    (a::Argument)(em::MixEscaped) -> Expr

Convert an `Argument` to the `Expr` representation of the argument it describes.
"""
(a::Argument)(em::RawExpr)=combinearg(a.name,a.type(em),a.slurp,escape(a.default,em))
(a::Argument)(em::MixEscaped)=combinearg(a.name,a.type(em.unescaped),a.slurp,escape(a.default,em.escaped))

"""
    Parameter(head::Union{Symbol,Nothing},name::Union{Symbol,Nothing},type::Union{Inference,Nothing})
    Parameter(;head::Union{Symbol,Nothing}=nothing,name::Union{Symbol,Nothing}=nothing,type::Union{Inference,Nothing}=nothing)
    Parameter(expr::FExpr)

The struct to describe a parameter of a `function` or a `type`.
"""
mutable struct Parameter <: AbstractFactory
    head::Union{Symbol,Nothing}
    name::Union{Symbol,Nothing}
    type::Union{Inference,Nothing}
end
Parameter(;head::Union{Symbol,Nothing}=nothing,name::Union{Symbol,Nothing}=nothing,type::Union{Inference,Nothing}=nothing)=Parameter(head,name,type)
function Parameter(expr::FExpr)
    if isa(expr,Symbol)
        head,name,type=nothing,expr,nothing
    elseif expr.head==:(<:)
        head=:(<:)
        (name,type)=length(expr.args)==2 ? expr.args : (nothing,expr.args[1])
        type=Inference(type)
    elseif expr.head==:(::)
        @assert length(expr.args)==1 "Parameter error: wrong input expr."
        head,name,type=:(::),nothing,Inference(expr.args[1])
    else
        error("Parameter error: wrong input expr.")
    end
    Parameter(head,name,type)
end

"""
    @parameter expr::FExpr

Construct a `Parameter` directly from an parameter statement.
"""
macro parameter(expr::FExpr) expr=[expr];:(Parameter($expr...)) end

"""
    (p::Parameter)(em::RawExpr) -> FExpr
    (p::Parameter)(em::UnEscaped) -> FExpr
    (p::Parameter)(em::MixEscaped) -> FExpr

Convert a `Parameter` to the `Expr` representation of the parameter it describes.
"""
function (p::Parameter)(em::Union{RawExpr,<:UnEscaped,<:MixEscaped})
    p.name===nothing ? (p.head==:(::) ? p.type(em) : Expr(:(<:),p.type(em))) : p.type===nothing ? p.name : Expr(:(<:),p.name,p.type(em))
end

"""
    Field(name::Symbol,type::Inference)
    Field(;name::Symbol,type::FExpr=Inference(:Any))
    Field(expr::FExpr)

The struct to describe a field of a `struct`.
"""
mutable struct Field <: AbstractFactory
    name::Symbol
    type::Inference
end
Field(;name::Symbol,type::Inference=Inference(:Any))=Field(name,type)
Field(expr::FExpr)=(cache=splitarg(expr);Field(cache[1],Inference(cache[2])))

"""
    @field expr::FExpr

Construct a `Field` directly from a field statement.
"""
macro field(expr::FExpr) expr=[expr];:(Field($expr...)) end

"""
    (f::Field)(em::RawExpr) -> Expr
    (f::Field)(em::UnEscaped) -> Expr
    (f::Field)(em::MixEscaped) -> Expr

Convert a `Field` to the `Expr` representation of the field it describes.
"""
(f::Field)(em::Union{RawExpr,<:UnEscaped,<:MixEscaped})=combinefield((f.name,f.type(em)))

"""
    Block(parts::FExpr...)

The struct to describe a `begin ... end` block.
"""
mutable struct Block <: AbstractFactory
    body::Vector{Any}
    Block(parts::FExpr...)=new(vcat((block(part).args for part in parts)...))
end

"""
    @block parts::FExpr...

Construct a `Block` directly from a `begin ... end` block definition.
"""
macro block(parts::FExpr...) :(Block($parts...)) end

"""
    (b::Block)(em::RawExpr) -> Expr
    (b::Block)(em::Escaped) -> Expr
    (b::Block)(em::MixEscaped) -> Expr

Convert a `Block` to the `Expr` representation of the `begin ... end` block it describes.
"""
(b::Block)(em::MixEscaped)=b(em.escaped)
(b::Block)(em::Union{RawExpr,<:Escaped})=Expr(:block,(escape(part,em) for part in b.body)...)

"""
    push!(b::Block,parts::FExpr...) -> Block
    push!(b::Block,parts::Block...) -> Block

Push other parts into the body of a block.
"""
Base.push!(b::Block)=b
Base.push!(b::Block,parts::FExpr...)=push!(b,Block(parts...))
function Base.push!(b::Block,parts::Block...)
    for part in parts
        append!(b.body,part.body)
    end
    return b
end

"""
    @push! b parts::FExpr...

Push other parts into the body of a block.
"""
macro push!(b,parts::FExpr...) :(push!($(esc(b)),$parts...)) end

"""
    rmlines!(b::Block) -> Block

Remove line number nodes in the body of a block.
"""
rmlines!(b::Block)=(filter!(part->!isa(part,LineNumberNode),b.body);b)

"""
    @rmlines! b::Expr

Remove line number nodes in the body of a block.
"""
macro rmlines!(b::Expr) :(rmlines!($(esc(b)))) end

"""
    rmlines(b::Block) -> Block

Return a copy of a block with the line number nodes removed.
"""
rmlines(b::Block)=Block(filter(part->!isa(part,LineNumberNode),b.body)...)

"""
    @rmlines b::Expr

Return a copy of a block with the line number nodes removed.
"""
macro rmlines(b::Expr) :(rmlines($(esc(b)))) end

"""
    FunctionFactory(name::FExpr,params::Vector{Inference},args::Vector{Argument},kwargs::Vector{Argument},rtype::Inference,whereparams::Vector{Parameter},body::Block)
    FunctionFactory(    ;name::FExpr,
                        params::Vector{Inference}=Inference[],
                        args::Vector{Argument}=Argument[],
                        kwargs::Vector{Argument}=Argument[],
                        rtype::Inference=Inference(:Any),
                        whereparams::Vector{Parameter}=Parameter[],
                        body::Block=Block()
                        )
    FunctionFactory(expr::Expr)

The struct to describe a `function`.
"""
mutable struct FunctionFactory <: AbstractFactory
    name::FExpr
    params::Vector{Inference}
    args::Vector{Argument}
    kwargs::Vector{Argument}
    rtype::Inference
    whereparams::Vector{Parameter}
    body::Block
end
function FunctionFactory(   ;name::FExpr,
                            params::Vector{Inference}=Inference[],
                            args::Vector{Argument}=Argument[],
                            kwargs::Vector{Argument}=Argument[],
                            rtype::Inference=Inference(:Any),
                            whereparams::Vector{Parameter}=Parameter[],
                            body::Block=Block(),
                            )
    FunctionFactory(name,params,args,kwargs,rtype,whereparams,body)
end
function FunctionFactory(expr::Expr)
    dict=splitdef(expr)
    FunctionFactory(    dict[:name],
                        Inference.(get(dict,:params,[])),
                        Argument.(dict[:args]),
                        Argument.(dict[:kwargs]),
                        Inference(get(dict,:rtype,:Any)),
                        Parameter.(collect(FExpr,dict[:whereparams])),
                        Block(dict[:body]),
                        )
end

"""
    @functionfactory expr::FExpr

Construct a `FunctionFactory` directly from a function definition.
"""
macro functionfactory(expr::Expr) expr=[expr];:(FunctionFactory($expr...)) end

"""
    (ff::FunctionFactory)(em::RawExpr) -> Expr
    (ff::FunctionFactory)(em::MixEscaped) -> Expr

Convert a `FunctionFactory` to the `Expr` representation of the `function` it describes.
"""
function (ff::FunctionFactory)(em::Union{RawExpr,<:MixEscaped})
    combinedef(Dict(
        :name           =>      escape(ff.name,isa(em,RawExpr) ? em : em.escaped),
        :params         =>      [param(em) for param in ff.params],
        :args           =>      [arg(em) for arg in ff.args],
        :kwargs         =>      [kwarg(em) for kwarg in ff.kwargs],
        :rtype          =>      ff.rtype(em),
        :whereparams    =>      [whereparam(em) for whereparam in ff.whereparams],
        :body           =>      ff.body(em)
    ))
end

"""
    addargs!(ff::FunctionFactory,args::Argument...) -> FunctionFactory
    addargs!(ff::FunctionFactory,args::FExpr...) -> FunctionFactory

Add a couple of positional arguments to a function factory.
"""
addargs!(ff::FunctionFactory)=ff
addargs!(ff::FunctionFactory,args::Argument...)=(push!(ff.args,args...);ff)
addargs!(ff::FunctionFactory,args::FExpr...)=addargs!(ff,Argument.(args)...)

"""
    @addargs! ff args::FExpr...

Add a couple of positional arguments to a function factory.
"""
macro addargs!(ff,args::FExpr...) :(addargs!($(esc(ff)),$args...)) end

"""
    addkwargs!(ff::FunctionFactory,kwargs::Argument...) -> FunctionFactory
    addkwargs!(ff::FunctionFactory,kwargs::FExpr...) -> FunctionFactory

Add a couple of keyword arguments to a function factory.
"""
addkwargs!(ff::FunctionFactory)=ff
addkwargs!(ff::FunctionFactory,kwargs::Argument...)=(push!(ff.kwargs,kwargs...);ff)
addkwargs!(ff::FunctionFactory,kwargs::FExpr...)=addkwargs!(ff,Argument.(kwargs)...)

"""
    @addkwargs! ff kwargs::FExpr...

Add a couple of keyword arguments to a function factory.
"""
macro addkwargs!(ff,kwargs::FExpr...) :(addkwargs!($(esc(ff)),$kwargs...)) end

"""
    extendbody!(ff::FunctionFactory,parts::FExpr...) -> FunctionFactory
    extendbody!(ff::FunctionFactory,parts::Block...) -> FunctionFactory

Extend the body of a function factory.
"""
extendbody!(ff::FunctionFactory)=ff
extendbody!(ff::FunctionFactory,parts::FExpr...)=extendbody!(ff,Block(parts...))
extendbody!(ff::FunctionFactory,parts::Block...)=(push!(ff.body,parts...);ff)

"""
    @extendbody! ff parts::FExpr...

Extend the body of a function factory.
"""
macro extendbody!(ff,parts::FExpr...) :(extendbody!($(esc(ff)),$parts...)) end

"""
    addwhereparams!(ff::FunctionFactory,whereparams::Parameter...) -> FunctionFactory
    addwhereparams!(ff::FunctionFactory,whereparams::FExpr...) -> FunctionFactory

Add a couple of method where parameters to a function factory or a type factory.
"""
addwhereparams!(ff::FunctionFactory)=ff
addwhereparams!(ff::FunctionFactory,whereparams::Parameter...)=(push!(ff.whereparams,whereparams...);ff)
addwhereparams!(ff::FunctionFactory,whereparams::FExpr...)=addwhereparams!(ff,Parameter.(whereparams)...)

"""
    @addwhereparams! ff whereparams::FExpr...

Add a couple of method parameters to a function factory or a type factory.
"""
macro addwhereparams!(ff,whereparams::FExpr...) :(addwhereparams!($(esc(ff)),$whereparams...)) end

"""
    TypeFactory(name::Symbol,mutable::Bool,params::Vector{Parameter},supertype::Inference,fields::Vector{Field},constructors::Vector{FunctionFactory})
    TypeFactory(    ;name::Symbol,
                    mutable::Bool=false,
                    params::Vector{Parameter}=Parameter[],
                    supertype::Inference=Inference(:Any),
                    fields::Vector{Field}=Field[],
                    constructors::Vector{FunctionFactory}=FunctionFactory[],
                    )
    TypeFactory(expr::Expr)

The struct to describe a `struct`.
"""
mutable struct TypeFactory <: AbstractFactory
    name::Symbol
    mutable::Bool
    params::Vector{Parameter}
    supertype::Inference
    fields::Vector{Field}
    constructors::Vector{FunctionFactory}
end
function TypeFactory(   ;name::Symbol,
                        mutable::Bool=false,
                        params::Vector{Parameter}=Parameter[],
                        supertype::Inference=Inference(:Any),
                        fields::Vector{Field}=Field[],
                        constructors::Vector{FunctionFactory}=FunctionFactory[],
                        )
    TypeFactory(name,mutable,params,supertype,fields,constructors)
end
function TypeFactory(expr::Expr)
    dict=splitstructdef(expr)
    TypeFactory(    dict[:name],
                    dict[:mutable],
                    Parameter.(dict[:params]),
                    Inference(dict[:supertype]),
                    [Field(field[1],Inference(field[2])) for field in dict[:fields]],
                    FunctionFactory.(dict[:constructors])
    )
end

"""
    @typefactory expr::Expr

Construct a `TypeFactory` directly from a type definition.
"""
macro typefactory(expr::Expr) expr=[expr];:(TypeFactory($expr...)) end

"""
    (tf::TypeFactory)(em::RawExpr) -> Expr
    (tf::TypeFactory)(em::MixEscaped) -> Expr

Convert a `TypeFactory` to the `Expr` representation of the `struct` it describes.
"""
function (tf::TypeFactory)(em::Union{RawExpr,<:MixEscaped})
    combinestructdef(Dict(
        :name           =>      escape(tf.name,isa(em,RawExpr) ? em : em.escaped),
        :mutable        =>      tf.mutable,
        :params         =>      [param(em) for param in tf.params],
        :supertype      =>      tf.supertype(em),
        :fields         =>      [(field.name,field.type(em)) for field in tf.fields],
        :constructors   =>      [constructor(em) for constructor in tf.constructors]
        ))
end

"""
    addfields!(tf::TypeFactory,fields::Field...) -> TypeFactory
    addfields!(tf::TypeFactory,fields::FExpr...) -> TypeFactory

Add a couple of fields to a type factory.
"""
addfields!(tf::TypeFactory)=tf
addfields!(tf::TypeFactory,fields::Field...)=(push!(tf.fields,fields...);tf)
addfields!(tf::TypeFactory,fields::FExpr...)=addfields!(tf,Field.(fields)...)

"""
    @addfields! tf fields::FExpr...

Add a couple of fields to a type factory.
"""
macro addfields!(tf,fields::FExpr...) :(addfields!($(esc(tf)),$fields...)) end

"""
    addconstructors!(tf::TypeFactory,constructors::FunctionFactory...) -> TypeFactory
    addconstructors!(tf::TypeFactory,constructors::Expr...) -> TypeFactory

Add a couple of constructors to a type factory.
"""
addconstructors!(tf::TypeFactory)=tf
addconstructors!(tf::TypeFactory,constructors::FunctionFactory...)=(push!(tf.constructors,constructors...);tf)
addconstructors!(tf::TypeFactory,constructors::Expr...)=addconstructors!(tf,FunctionFactory.(constructors)...)

"""
    @addconstructors! tf constructors::Expr...

Add a couple of constructors to a type factory.
"""
macro addconstructors!(tf,constructors::Expr...) :(addconstructors!($(esc(tf)),$constructors...)) end

"""
    addparams!(f::FunctionFactory,params::Inference...) ->FunctionFactory
    addparams!(f::FunctionFactory,params::FExpr...) -> FunctionFactory
    addparams!(f::TypeFactory,params::Parameter...) -> TypeFactory
    addparams!(f::TypeFactory,params::FExpr...) -> TypeFactory

Add a couple of parameters to a function factory or a type factory.
"""
addparams!(f::FunctionFactory)=f
addparams!(f::FunctionFactory,params::Inference...)=(push!(f.params,params...);f)
addparams!(f::FunctionFactory,params::FExpr...)=addparams!(f,Inference.(params)...)
addparams!(f::TypeFactory)=f
addparams!(f::TypeFactory,params::Parameter...)=(push!(f.params,params...);f)
addparams!(f::TypeFactory,params::FExpr...)=addparams!(f,Parameter.(params)...)

"""
    @addparams! f params::FExpr...

Add a couple of method parameters to a function factory or a type factory.
"""
macro addparams!(f,params::FExpr...) :(addparams!($(esc(f)),$params...)) end

end #module


text:
using Plots
include("model.jl")
include("second_order_cone.jl")
include("step.jl")
include("simulate.jl")
include("visualize.jl")

# horizon
T = 201

# time step
h = 0.01
t = range(0, stop = h * (T - 1), length = T)

mrp = MRP(UnitQuaternion(RotY(0.0 * π / 6.0) * RotX(0.0 * π / 10.0)))

# initial conditions
v1 = [2.5; 5.0; 0.0; 0.0; 0.0; 0.0]
q1 = [0.0; 0.0; 1.0; mrp.x; mrp.y; mrp.z]

v2 = v1 - gravity(model) * h
q2 = q1 + 0.5 * (v1 + v2) * h

# simulate
q_sol, y_sol, b_sol, Δq1, Δq2, Δu1 = simulate(q1, q2, T, h)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

visualize!(vis, model,
    q_sol,
    Δt = h)

settransform!(vis["/Cameras/default"],
    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(pi))))

settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, 0.0, 3.0),LinearMap(RotY(-pi/2.5))))

## trajectory optimization
T = 101
h = 0.01
ul, uu = control_bounds(model, T)

mrp_init = MRP(UnitQuaternion(RotY(0.0) * RotX(0.0)))
mrp_target = MRP(UnitQuaternion(RotZ(π / 2.0)))
q1 = [model.r, model.r, model.r + 1.0, mrp_init.x, mrp_init.y, mrp_init.z]
# qT = [0.0, 0.0, model.r * sqrt(3.0), mrp_corner.x, mrp_corner.y, mrp_corner.z]
qT = copy(q1)
qT[1] = 2.5
# qT[2] = 1.0
qT[4:6] = [mrp_target.x; mrp_target.y; mrp_target.z]
visualize!(vis, model,
    [qT],
    Δt = h)
x1 = [q1; q1]
xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(1.0 * ones(model.n)) : 10.0 * Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [[qT; qT] for t = 1:T],
        [zeros(model.m) for t = 1:T])

# Constraints
include("simulator_dynamics.jl")
con_sim = simulator_constraints(model, T)

# Problem

prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_sim,
			   dynamics = false)

# Trajectory initialization
x0 = configuration_to_state([q1, linear_interpolation(q1, qT, T)...]) #linear_interpolation(x1, x1, T) # linear interpolation on state
# u0 = [1.0 * [0.0; 0.0; 0.0; t == 13 ? 0.0 : 0.0; t == 1 ? 0.0 : 0.0; 0.0] for t = 1:T-1] # random controls
u0 = [[1.0 * [5.0; 0.0; 0.0]; 0.01 * rand(3)] for t = 1:T-1]
# step(q1, q1, u0[1], h;
#     tol = 1.0e-8, max_iter = 100, z_init = 10.0)
q_sol, y_sol, b_sol, Δq1, Δq2, Δu1 = simulate(q1, q1, T-1, h, u1 = u0, r_tol = 1.0e-5, μ_tol = 1.0e-5, z_init = 1.0)

visualize!(vis, model,
    q_sol,
    Δt = h)

# Pack trajectories into vector
x0 = configuration_to_state(q_sol) #linear_interpolation(x1, x1, T) # linear interpolation on state
z0 = pack(x0, u0, prob)
prob.num_var
prob.num_con
#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()
@time z̄, info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-2, c_tol = 1.0e-2)

cc = zeros(con_sim.n)
constraints!(cc, z̄, con_sim, model, prob.prob.idx, h, T)
norm(cc, Inf)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)
t = range(0, stop = h * (T - 1), length = T)

plot(t, hcat(q̄[2:end]...)', xlabel = "time (s)", ylabel = "", label = "")
plot(t, hcat(ū...,ū[end])',
	xlabel = "time (s)" , ylabel = "", linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

visualize!(vis, model,
    [[q̄[1] for i = 1:10] ..., q̄...],
    Δt = h)

@show q̄[end]

settransform!(vis["/Cameras/default"],
    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(-pi))))

settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, 0.0, 1.0),LinearMap(RotY(-pi/2.5))))

setobject!(vis["particle_goal"],
	Rect(Vec(0, 0, 0),Vec(0.2, 0.2, 0.2)),
	MeshPhongMaterial(color = RGBA(0.0, 1.0, 1.0, 0.5)))

settransform!(vis["particle_goal"], Translation(qT...))

open(vis)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the API to convert from Euler angles to rotation
#   representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_rot

"""
    angle_to_rot([T,] θ₁::Number[, θ₂::Number[, θ₃::Number]], rot_seq::Symbol)
    angle_to_rot([T,] Θ::EulerAngles)

Create a rotation description of type `T` that perform a set of rotations (`θ₁`,
`θ₂`, `θ₃`) about the coordinate axes specified in `rot_seq`.

The input values of the origin Euler angles can also be passed inside the
structure `Θ` (see [`EulerAngles`](@ref)).

The rotation sequence is defined by a `Symbol` specifing the rotation axes. The
possible values depends on the number of rotations as follows:

- **1 rotation** (`θ₁`): `:X`, `:Y`, or `:Z`.
- **2 rotations** (`θ₁`, `θ₂`): `:XY`, `:XZ`, `:YX`, `:YZ`, `:ZX`, or `:ZY`.
- **3 rotations** (`θ₁`, `θ₂`, `θ₃`): `:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`,
    `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`, `:ZYX`, or `:ZYZ`

# Example

```julia-repl
julia> dcm = angle_to_rot(pi / 5, :Z)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.809017  0.587785  0.0
 -0.587785  0.809017  0.0
  0.0       0.0       1.0

julia> quat = angle_to_rot(Quaternion, pi / 5, :Z)
Quaternion{Float64}:
  + 0.951057 + 0.0⋅i + 0.0⋅j + 0.309017⋅k

julia> dcm = angle_to_rot(pi / 5, pi / 7, :YZ)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.728899  0.433884  -0.529576
 -0.351019  0.900969   0.25503
  0.587785  0.0        0.809017

julia> quat = angle_to_rot(Quaternion, pi / 5, pi / 7, :YZ)
Quaternion{Float64}:
  + 0.927212 + 0.0687628⋅i + 0.301269⋅j + 0.21163⋅k

julia> dcm = angle_to_rot(pi / 2, pi / 3, pi / 4, :ZYX)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553

julia> q = angle_to_rot(Quaternion, pi / 2, pi / 3, pi / 4, :ZYX)
Quaternion{Float64}:
  + 0.701057 - 0.092296⋅i + 0.560986⋅j + 0.430459⋅k
```
"""
@inline function angle_to_rot(θ::Number, rot_seq::Symbol)
    return angle_to_dcm(θ, rot_seq)
end

@inline function angle_to_rot(::Type{DCM}, θ::Number, rot_seq::Symbol)
    return angle_to_dcm(θ, rot_seq)
end

@inline function angle_to_rot(::Type{Quaternion}, θ::Number, rot_seq::Symbol)
    return angle_to_quat(θ, rot_seq)
end

@inline function angle_to_rot(
    θ₁::Number,
    θ₂::Number,
    rot_seq::Symbol
)
    return angle_to_dcm(θ₁, θ₂, rot_seq)
end

@inline function angle_to_rot(
    ::Type{DCM},
    θ₁::Number,
    θ₂::Number,
    rot_seq::Symbol
)
    return angle_to_dcm(θ₁, θ₂, rot_seq)
end

@inline function angle_to_rot(
    ::Type{Quaternion},
    θ₁::Number,
    θ₂::Number,
    rot_seq::Symbol
)
    return angle_to_quat(θ₁, θ₂, rot_seq)
end

@inline function angle_to_rot(
    θ₁::Number,
    θ₂::Number,
    θ₃::Number,
    rot_seq::Symbol
)
    return angle_to_dcm(θ₁, θ₂, θ₃, rot_seq)
end

@inline function angle_to_rot(
    ::Type{DCM},
    θ₁::Number,
    θ₂::Number,
    θ₃::Number,
    rot_seq::Symbol
)
    return angle_to_dcm(θ₁, θ₂, θ₃, rot_seq)
end

@inline function angle_to_rot(
    ::Type{Quaternion},
    θ₁::Number,
    θ₂::Number,
    θ₃::Number,
    rot_seq::Symbol
)
    return angle_to_quat(θ₁, θ₂, θ₃, rot_seq)
end

@inline function angle_to_rot(Θ::EulerAngles)
    return angle_to_rot(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)
end

@inline function angle_to_rot(
    T::Union{Type{DCM}, Type{Quaternion}},
    Θ::EulerAngles
)
    return angle_to_rot(T, Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)
end


text:
# MOI wrapper

"""
    MOIEvaluator <: MOI.AbstractNLPEvaluator

Bridge from a `ExaPF.AbstractNLPEvaluator` to a `MOI.AbstractNLPEvaluator`.

## Attributes

* `nlp::AbstractNLPEvaluator`: the underlying `ExaPF` problem.
* `hash_x::UInt`: hash of the last evaluated variable `x`
* `has_hess::Bool` (default: `false`): if `true`, pass a Hessian structure to MOI.

"""
mutable struct MOIEvaluator <: MOI.AbstractNLPEvaluator
    nlp::AbstractNLPEvaluator
    hash_x::UInt
    has_hess::Bool
end
MOIEvaluator(nlp) = MOIEvaluator(nlp, UInt64(0), false)

function _update!(ev::MOIEvaluator, x)
    hx = hash(x)
    if hx != ev.hash_x
        update!(ev.nlp, x)
        ev.hash_x = hx
    end
end

MOI.features_available(ev::MOIEvaluator) = ev.has_hess ? [:Grad, :Hess] : [:Grad]
MOI.initialize(ev::MOIEvaluator, features) = nothing

function MOI.jacobian_structure(ev::MOIEvaluator)
    n = n_variables(ev.nlp)
    m = n_constraints(ev.nlp)
    jnnz = n * m
    rows = zeros(Int, jnnz)
    cols = zeros(Int, jnnz)
    jacobian_structure!(ev.nlp, rows, cols)
    return Tuple{Int, Int}[(r, c) for (r, c) in zip(rows, cols)]
end

function MOI.hessian_lagrangian_structure(ev::MOIEvaluator)
    n = n_variables(ev.nlp)
    return Tuple{Int, Int}[(i, i) for i in 1:n]
end

function MOI.eval_objective(ev::MOIEvaluator, x)
    _update!(ev, x)
    obj = objective(ev.nlp, x)
    return obj
end

function MOI.eval_objective_gradient(ev::MOIEvaluator, g, x)
    _update!(ev, x)
    gradient!(ev.nlp, g, x)
end

function MOI.eval_constraint(ev::MOIEvaluator, cons, x)
    _update!(ev, x)
    constraint!(ev.nlp, cons, x)
end

function MOI.eval_constraint_jacobian(ev::MOIEvaluator, ∂cons, x)
    _update!(ev, x)
    n = n_variables(ev.nlp)
    m = n_constraints(ev.nlp)
    jac = zeros(m, n)
    jacobian!(ev.nlp, jac, x)
    k = 1
    for i = 1:m
        for j = 1:n
            ∂cons[k] = jac[i, j]
            k += 1
        end
    end
end

function MOI.eval_hessian_lagrangian(ev::MOIEvaluator, H, x, σ, μ)
    fill!(H, 1.0)
end

function MOI.NLPBlockData(nlp::AbstractNLPEvaluator)
    lb, ub = bounds(nlp, Constraints())
    ev = MOIEvaluator(nlp)
    return MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), ev, true)
end


text:
using Test
using Printf
using Logging
using Oceananigans
using Oceananigans.Fields
using Oceananigans.OutputWriters
using LESbrary

Logging.global_logger(OceananigansLogger())

architectures = [CPU()]

@testset "LESbrary" begin
    include("test_fit_cubic.jl")
    include("test_interpolated_profiles.jl")
    include("test_diagnose_buoyancy_flux.jl")
    include("test_turbulence_statistics.jl")
    include("test_examples.jl")
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: dynamodb_streams
using AWS.Compat
using AWS.UUIDs

"""
    describe_stream(stream_arn)
    describe_stream(stream_arn, params::Dict{String,<:Any})

Returns information about a stream, including the current status of the stream, its Amazon
Resource Name (ARN), the composition of its shards, and its corresponding DynamoDB table.
You can call DescribeStream at a maximum rate of 10 times per second.  Each shard in the
stream has a SequenceNumberRange associated with it. If the SequenceNumberRange has a
StartingSequenceNumber but no EndingSequenceNumber, then the shard is still open (able to
receive more stream records). If both StartingSequenceNumber and EndingSequenceNumber are
present, then that shard is closed and can no longer receive more data.

# Arguments
- `stream_arn`: The Amazon Resource Name (ARN) for the stream.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ExclusiveStartShardId"`: The shard ID of the first item that this operation will
  evaluate. Use the value that was returned for LastEvaluatedShardId in the previous
  operation.
- `"Limit"`: The maximum number of shard objects to return. The upper limit is 100.
"""
function describe_stream(StreamArn; aws_config::AbstractAWSConfig=global_aws_config())
    return dynamodb_streams(
        "DescribeStream",
        Dict{String,Any}("StreamArn" => StreamArn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_stream(
    StreamArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dynamodb_streams(
        "DescribeStream",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("StreamArn" => StreamArn), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_records(shard_iterator)
    get_records(shard_iterator, params::Dict{String,<:Any})

Retrieves the stream records from a given shard. Specify a shard iterator using the
ShardIterator parameter. The shard iterator specifies the position in the shard from which
you want to start reading stream records sequentially. If there are no stream records
available in the portion of the shard that the iterator points to, GetRecords returns an
empty list. Note that it might take multiple calls to get to a portion of the shard that
contains stream records.   GetRecords can retrieve a maximum of 1 MB of data or 1000 stream
records, whichever comes first.

# Arguments
- `shard_iterator`: A shard iterator that was retrieved from a previous GetShardIterator
  operation. This iterator can be used to access the stream records in this shard.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Limit"`: The maximum number of records to return from the shard. The upper limit is
  1000.
"""
function get_records(ShardIterator; aws_config::AbstractAWSConfig=global_aws_config())
    return dynamodb_streams(
        "GetRecords",
        Dict{String,Any}("ShardIterator" => ShardIterator);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_records(
    ShardIterator,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dynamodb_streams(
        "GetRecords",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ShardIterator" => ShardIterator), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_shard_iterator(shard_id, shard_iterator_type, stream_arn)
    get_shard_iterator(shard_id, shard_iterator_type, stream_arn, params::Dict{String,<:Any})

Returns a shard iterator. A shard iterator provides information about how to retrieve the
stream records from within a shard. Use the shard iterator in a subsequent GetRecords
request to read the stream records from the shard.  A shard iterator expires 15 minutes
after it is returned to the requester.

# Arguments
- `shard_id`: The identifier of the shard. The iterator will be returned for this shard ID.
- `shard_iterator_type`: Determines how the shard iterator is used to start reading stream
  records from the shard:    AT_SEQUENCE_NUMBER - Start reading exactly from the position
  denoted by a specific sequence number.    AFTER_SEQUENCE_NUMBER - Start reading right after
  the position denoted by a specific sequence number.    TRIM_HORIZON - Start reading at the
  last (untrimmed) stream record, which is the oldest record in the shard. In DynamoDB
  Streams, there is a 24 hour limit on data retention. Stream records whose age exceeds this
  limit are subject to removal (trimming) from the stream.    LATEST - Start reading just
  after the most recent stream record in the shard, so that you always read the most recent
  data in the shard.
- `stream_arn`: The Amazon Resource Name (ARN) for the stream.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"SequenceNumber"`: The sequence number of a stream record in the shard from which to
  start reading.
"""
function get_shard_iterator(
    ShardId, ShardIteratorType, StreamArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dynamodb_streams(
        "GetShardIterator",
        Dict{String,Any}(
            "ShardId" => ShardId,
            "ShardIteratorType" => ShardIteratorType,
            "StreamArn" => StreamArn,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_shard_iterator(
    ShardId,
    ShardIteratorType,
    StreamArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dynamodb_streams(
        "GetShardIterator",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ShardId" => ShardId,
                    "ShardIteratorType" => ShardIteratorType,
                    "StreamArn" => StreamArn,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_streams()
    list_streams(params::Dict{String,<:Any})

Returns an array of stream ARNs associated with the current account and endpoint. If the
TableName parameter is present, then ListStreams will return only the streams ARNs for that
table.  You can call ListStreams at a maximum rate of 5 times per second.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ExclusiveStartStreamArn"`: The ARN (Amazon Resource Name) of the first item that this
  operation will evaluate. Use the value that was returned for LastEvaluatedStreamArn in the
  previous operation.
- `"Limit"`: The maximum number of streams to return. The upper limit is 100.
- `"TableName"`: If this parameter is provided, then only the streams associated with this
  table name are returned.
"""
function list_streams(; aws_config::AbstractAWSConfig=global_aws_config())
    return dynamodb_streams(
        "ListStreams"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_streams(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dynamodb_streams(
        "ListStreams", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end


text:
import CSV
import DataFrames

symbols = [:SOVEREIGNT, :GDP_MD_EST, :GDP_YEAR]
path = "data/natural_earth_vector"
attributesname = "temp2-attributes.csv"
nodesname = "temp1-nodes.csv"
attributes = DataFrames.DataFrame(CSV.File(joinpath(path, attributesname)))
nodes = DataFrames.DataFrame(CSV.File(joinpath(path, nodesname)))
dict = Dict()
for symbol in symbols
    column = DataFrames.groupby(attributes, symbol)
    dict[symbol] = column
end
gdp_year = DataFrames.groupby(attributes, :GDP_YEAR)
subdataframe = groupdataframe[(SOVEREIGNT="Iran",)]
gdp = subdataframe.GDP_MD_EST

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# Doing epsilon greedy
# https://youtu.be/PnHCvfgC_ZA?t=3013

using Game2048
using Flux

const EPSILON = 0.5
const LAMBDA = 0.9 # the discount fcator

function mirrors_rotations(board)
    tboard = transpose(board) |> collect

    return board,  rotl90(board),  rot180(board),  rotr90(board),
           tboard, rotl90(tboard), rot180(tboard), rotr90(tboard)
end

actor_model = Chain(
    board -> Float32.(board ./ maximum(board)),
    normalised_board -> reshape(normalised_board, 4, 4, 1, 1),
    Conv((2,2), 1=>64, relu),
    Conv((2,2), 64=>128, relu),
    x->reshape(x, 2*2*128),
    Dense(128*2*2, 4),
    softmax
)

critic_model1 = Chain(
    # boardmove -> (boardmove[1] ./ maximum(boardmove[1]), boardmove[2]),
    boardmove -> (reshape(boardmove[1], 4, 4, 1, 1), boardmove[2]),
    boardmove -> (Conv((2,2), 1=>64, relu)(boardmove[1]), boardmove[2]),
    boardmove -> (Conv((2,2), 64=>128, relu)(boardmove[1]), boardmove[2]),
    boardmove -> vcat(reshape(boardmove[1], 2*2*128), boardmove[2]),
    Dense(2*2*128+4, 256, relu),
    # Dense(256, 1, exp)
)

critic_model2 = Chain(
    Dense(2*2*128+4, 256, relu),
    Dense(256, 1, exp)
)

function action2vec(dir)
    a = zeros(Float32, 4)
    a[Int(dir)+1] = 1
    a
end

board = initboard()

p = Flux.params(critic_model1)
normalised_board(board) = Float32.(board ./ maximum(board))
critic_model1([board |> normalised_board, action2vec(left)])

loss_qa(boardmove, expected_value) = begin
    #sum(sum.(boardmove))
    oldSA = critic_model1(boardmove)
    sum(oldSA)
    #sum(oldSA)
    #sum((oldSA .- expected_value).^2)
    # sum((oldSA .- (reward .+ LAMBDA.*newSA)).^2)

end

opt = ADAM()
normalised_board(board) = Float32.(board ./ maximum(board))
X = [board |> normalised_board, action2vec(left)]
Y = Float32[10.0]
loss_qa(X, Y)

Flux.train!(loss_qa, p, [(X, Y)], opt)

# S,A
# R = Reward
# S'
# A'
# = SARSA

res = move.(Ref(board), DIRS)
possible = map(r->r[2][2], res)
prob_dir = actor_model(board)[possible]
_, greedy_pos = findmax(prob_dir)
reward = map(r->r[2][1], res)[possible][greedy_pos]
dir = rand() < EPSILON ? rand(DIRS[possible]) : DIRS[possible][prob_dir]

features = critic_model1([board, action2vec(dir)])
SA = critic_model2(vcat(features, action2vec(dir)))

#  up to here I have the SAR
new_state = res[possible][greedy_pos][1]
gen_new_tile!(new_state)
res = move.(Ref(new_state), DIRS)
possible = map(r->r[2][2], res)
prob_dir = actor_model(new_state)[possible]
_, greedy_pos = findmax(prob_dir)
reward = map(r->r[2][1], res)[possible][greedy_pos]
dir = rand() < EPSILON ? rand(DIRS[possible]) : DIRS[possible][greedy_pos]
features = critic_model1(new_state)
newSA = critic_model2(vcat(features, action2vec(dir)))

action = dir
loss_qa(boardmove, expected_value) = begin
    oldSA = critic_model1(boardmove)
    sum((oldSA .- estimate_reward).^2)
    # sum((oldSA .- (reward .+ LAMBDA.*newSA)).^2)
end

estimate_reward = reward + LAMBDA*newSA[1]

p = Flux.params(critic_model1)
opt = ADAM()
loss_qa(([board, action2vec(dir)], estimate_reward)...)
Flux.train!(loss_qa, p, [([board, action2vec(dir)], estimate_reward)], opt)

@time actor_model(board)
@time critic_model.(map(res->res[1], res))

using Flux:crossentropy
crossentropy([0.97, 0.01, 0.01, 0.01],[ 1, 0, 0, 0])


text:
module HumanReadableIDs

export encode_proquint, decode_proquint

using Distributions
using Random

_normalize(x) = x./sum(x)

#-------------------------------------------------------------------------------
# Tools to generate transition probabilities
function letter_transition_frequencies(words)
    # Build (right) stochastic matrix
    # https://en.wikipedia.org/wiki/Stochastic_matrix<Paste>
    char_to_ind(c) = c - 'a' + 1

    nletters = 26
    nstates = nletters
    #end_ind = nletters+1
    start_letter_freqs = zeros(Int, nletters)
    transition_freqs = zeros(Int, nstates,nstates)
    for word in words
        isempty(word) && continue
        start_letter_freqs[char_to_ind(word[1])] += 1
        for i = 1:length(word)-1
            i1 = char_to_ind(word[i])
            i2 = char_to_ind(word[nextind(word,i)])
            transition_freqs[i1,i2] += 1
        end
        # End state
        # transition_freqs[char_to_ind(last(word)),end_ind] += 1
    end

    start_letter_freqs, transition_freqs
end

function letter_transition_probabilites(words)
    start_letter_freqs, transition_freqs = letter_transition_frequencies(words)
    (_normalize(start_letter_freqs),
     [_normalize(transition_freqs[i,:]) for i=1:size(transition_freqs, 1)])
end

function read_name_corpus()
    names=String[]
    # Name corups from
    # https://www.cs.cmu.edu/Groups/AI/areas/nlp/corpora/names/male.txt
    # https://www.cs.cmu.edu/Groups/AI/areas/nlp/corpora/names/female.txt
    #
    # Also interesting, much larger corpus:
    # https://www.cs.cmu.edu/Groups/AI/areas/nlp/corpora/names/other/names.txt
    for line in vcat(readlines(joinpath(@__DIR__, "../notebooks/male.txt")),
                     readlines(joinpath(@__DIR__, "../notebooks/female.txt")))
        if isempty(line) || first(line) == "#"
            continue
        end
        append!(names, lowercase.(split(line, r"[^a-z]"i)))
    end
    names
end

#-------------------------------------------------------------------------------
# Markov chain sampling for Human-readable-ids (HRIDs)

struct MarkovHridSampler
    start_letter_sampler::Distributions.AliasTable
    pairwise_samplers::Vector{Distributions.AliasTable}
end

function MarkovHridSampler(start_letter_probs::AbstractVector,
                           pairwise_probs::AbstractVector)
    MarkovHridSampler(sampler(Categorical(start_letter_probs)),
                      sampler.(Categorical.(pairwise_probs)))
end

"""
Create a word-like ID of total length `len`, by running a Markov chain
with the initial letter probability `start_letter_probs` and transition
probabilities `pairwise_probs`.
"""
function Random.randstring(rng::AbstractRNG, m::MarkovHridSampler, len::Integer)
    nameinds = Int[]
    cind = rand(rng, m.start_letter_sampler)
    push!(nameinds, cind)
    for i=2:len
        cind = rand(rng, m.pairwise_samplers[cind])
        push!(nameinds, cind)
    end
    name = join('a' .+ nameinds .- 1)
    name = uppercase(name[1])*name[2:end]
    name
end

# As these are meant to be random IDs, they're securely generated with
# RandomDevice() by default.
Random.randstring(m::MarkovHridSampler, len::Integer) = randstring(RandomDevice(), m, len)

include("name_model.jl")

include("proquints.jl")

end


text:
using ArgParse
using JSON
using DataStructures

function main()
    s = ArgParseSettings()

    @add_arg_table s begin
        "dir"
            help="path to directory that contains the log files"
        "pattern"
            help="regex"

        "--n_chunks"
            default=50
            arg_type=Int64
    end

    args = parse_args(s, as_symbols=true)

    runs = Dict{Int64, Dict{Symbol, Any}}()

    n_chunks = args[:n_chunks]

    @info "Dir: $(args[:dir])"
    @info "Pattern: $(args[:pattern])"
    @info "Num of chunks: $(n_chunks)"

    pattern = Regex(args[:pattern])

    dataset_name = ""
    validation_name = ""
    test_name = ""
    space_script = ""
    method_name = ""

    for file_name in readdir(args[:dir])
        file_path = joinpath(args[:dir], file_name)

        !isfile(file_path) && continue
        !occursin(pattern, file_name) && continue

        info =  Dict{Symbol, Any}()

        info[:results_exec] = Vector()
        info[:validation_results] = Vector()

        @info file_name

        job_id = nothing

        open(file_path, "r") do file
            for l in eachline(file)
                l = replace(l, "[ Info: " => "")
                l = replace(l, r"^.*best +loss: +[0-9.]+\]" => "")

                o = nothing

                try
                    o = eval(Meta.parse(l))

                    if haskey(o, :x)
                        push!(info[:results_exec], o)
                        new_job_id = parse(Int64, match(r"validation_chunk_([0-9]+).txt", o[:args][:evaluation])[1])

                        if !isnothing(job_id) && job_id != new_job_id
                            @error("Different chunk ids!!")
                        end
                        job_id = new_job_id
                    end

                    if haskey(o, :vals)
                        push!(info[:validation_results], o)
                    end

                    if haskey(o, :evaluation)
                        info[:best_par] = o

                        if length(dataset_name) > 0 &&  dataset_name != o[:bug_dataset]
                            @error("Dataset is different")
                        end

                        validation_nm_tmp = replace(o[:evaluation], r"[0-9]+.txt" => ".txt")

                        if length(validation_name) > 0 &&  validation_name != validation_nm_tmp
                            @error("Validation is different")
                        end

                        test_nm_tmp = replace(o[:test], r"[0-9]+.txt" => ".txt")
                        if length(test_name) > 0 &&  test_name != test_nm_tmp
                            @error("Test is different")
                        end

                        if length(space_script) > 0 &&  space_script != o[:search_space_script]
                            @error("space_script is different")
                        end

                        if length(method_name) > 0 &&  method_name != o[:method_name]
                            @error("method is different")
                        end

                        dataset_name = o[:bug_dataset]
                        validation_name = validation_nm_tmp
                        test_name = test_nm_tmp
                        space_script = o[:search_space_script]
                        method_name = o[:method_name]
                    end

                    continue
                catch
                end

                try
                    o = JSON.parse(l)

                    if isa(o, Vector)
                        info[:test_results] = o
                    end
                catch
                end
            end
        end

        isnothing(job_id) && continue

        runs[job_id] = info
    end

    missing_ids = []
    best_parameters = Vector{Any}([nothing for i=1:n_chunks])

    test_results = OrderedDict("map"=>OrderedDict{Int64,Float64}(), "auc"=> OrderedDict{Int64,Float64}())

    for rr_k in 1:20
        test_results["rr@$(rr_k)"] = Dict{Int64,Float64}()
    end

    for i=0:n_chunks-1
        if !haskey(runs, i)
            push!(missing_ids, i)
            continue
        end

        validation_results = runs[i][:validation_results]

        if length(validation_results) == 0
            push!(missing_ids, i)
            continue
        end

        sort!(validation_results, by= j->j[:loss])
        best_parameters[i+1] = validation_results[1][:args]

        if !haskey(runs[i], :test_results)
            push!(missing_ids, i)
            continue
        end

        for r in runs[i][:test_results]
            if r["label"] == "MAP_RecallRate"
                test_results["map"][i] = r["map"]

                for (k,v) in r["rr"]
                    test_results["rr@$(k)"][i] = v
                end
            elseif r["label"] == "BinaryPredictionROC"
                test_results["auc"][i] = r["auc"]
            end
        end
    end

    @info [ i-1=>best_parameters[i] for i=1:length(best_parameters)]

    sort!(missing_ids)

    if length(missing_ids) > 0
        @info "Missing chunks (test): $(missing_ids)"

        println("Start commands:")
        for chunk_id in missing_ids

            if !isdefined(best_parameters, chunk_id+1)
                println("Chunk $(chunk_id) is undefined")

                continue
            end

            chunk_param = best_parameters[chunk_id+1]

            command = "julia -Jimage.so experiments/run_method.py $(chunk_param[:bug_dataset]) $(chunk_param[:test]) $(chunk_param[:method_name])"

            for (k, v) in chunk_param
                if k ∈ [:method_name, :bug_dataset, :evaluation, :test, :search_space_script, :max_evals]
                    continue
                end

                isnothing(v) && continue

                if isa(v, Bool)
                    if v
                        command = command * " --$(k)"
                    end
                else
                    command = command * " --$(k) $(v)"
                end

            end
            println(command)
        end
    end

    @info "$(method_name)"
    @info "$(dataset_name)\t$(validation_name)\t$(test_name)\t$(space_script)"

    @info JSON.json(test_results)
end

main()

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST