text:
# m6.7s.jl
# Clip-06-06-09s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

N = 100
df = DataFrame(
    :h0 => rand(Normal(10,2 ), N),
    :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
);
df[!, :fungus] = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
df[!, :h1] = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]

stan6_7 = "
data {
    int <lower=1> N;
    vector[N] h0;
    vector[N] h1;
    vector[N] treatment;
    vector[N] fungus;
}
parameters{
    real a;
    real bt;
    real bf;
    real<lower=0> sigma;
}
model {
    vector[N] mu;
    vector[N] p;
    a ~ lognormal(0, 0.2);
    bt ~ normal(0, 0.5);
    bf ~ normal(0, 0.5);
    sigma ~ exponential(1);
    for ( i in 1:N ) {
    p[i] = a + bt*treatment[i] + bf*fungus[i];
    mu[i] = h0[i] * p[i];
    }
    h1 ~ normal(mu, sigma);
}
"

m6_7_data = Dict(
    :N => nrow(df),
    :h0 => df[:, :h0],
    :h1 => df[:, :h1],
    :fungus => df[:, :fungus],
    :treatment => df[:, :treatment]
)

m6_7s = SampleModel("m6.7", stan6_7)
rc6_7s = stan_sample(m6_7s; data=m6_7_data)

if success(rc6_7s)
    part6_7s = read_samples(m6_7s, :particles);
    part6_7s |> display
end

# End of m6.7s.jl


text:
# Common regression problems

abstract Problem{T}

### Univariate Prediction Problem

immutable UnivariateRegression{L<:UnivariateLoss,
                               T<:AbstractFloat,
                               XT<:StridedMatrix,
                               YT<:StridedVector} <: Problem{T}
    loss::L
    d::Int
    n::Int
    bias::T
    X::XT
    Y::YT
end

function UnivariateRegression{T<:AbstractFloat}(loss::UnivariateLoss,
                                                X::StridedMatrix{T},
                                                Y::StridedVector;
                                                bias::Real=0.0)
    d, n = size(X)
    length(Y) == n || throw(DimensionMismatch())
    UnivariateRegression{typeof(loss), T, typeof(X), typeof(Y)}(
        loss, d, n, convert(T, bias), X, Y)
end

nsamples(pb::UnivariateRegression) = pb.n
inputs(pb::UnivariateRegression) = pb.X
outputs(pb::UnivariateRegression) = pb.Y

loss(pb::UnivariateRegression) = pb.loss

has_bias(pb::UnivariateRegression) = (pb.bias != 0)
pred_without_bias(pb::UnivariateRegression) = LinearPred(pb.d)
pred_with_bias(pb::UnivariateRegression) = AffinePred(pb.d, pb.bias)

initsol{L,T}(pb::UnivariateRegression{L,T}) = zeros(T, has_bias(pb) ? pb.d + 1: pb.d)

# specific problems

linearreg{T<:AbstractFloat}(X::StridedMatrix{T}, y::StridedVector{T}; bias::Real=0.0) =
    UnivariateRegression(SqrLoss(), X, y; bias=bias)

logisticreg{T<:AbstractFloat}(X::StridedMatrix{T}, y::StridedVector; bias::Real=0.0) =
    UnivariateRegression(LogisticLoss(), X, convert(Vector{T}, y); bias=bias)

### Multivariate Prediction Problem

immutable MultivariateRegression{L<:MultivariateLoss,
                                 T<:AbstractFloat,
                                 XT<:StridedMatrix,
                                 YT<:StridedVecOrMat} <: Problem{T}
    loss::L
    d::Int
    k::Int
    n::Int
    bias::T
    X::XT
    Y::YT
end

function MultivariateRegression{T<:AbstractFloat}(loss::MultivariateLoss,
                                                  X::StridedMatrix{T},
                                                  Y::StridedVecOrMat,
                                                  k::Int;
                                                  bias::Real=0.0)
    d, n = size(X)
    size(Y, ndims(Y)) == n || throw(DimensionMismatch())
    MultivariateRegression{typeof(loss), T, typeof(X), typeof(Y)}(
        loss, d, k, n, convert(T, bias), X, Y)
end

nsamples(pb::MultivariateRegression) = pb.n
inputs(pb::MultivariateRegression) = pb.X
outputs(pb::MultivariateRegression) = pb.Y

loss(pb::MultivariateRegression) = pb.loss

has_bias(pb::MultivariateRegression) = (pb.bias != 0)
pred_without_bias(pb::MultivariateRegression) = MvLinearPred(pb.d, pb.k)
pred_with_bias(pb::MultivariateRegression) = MvAffinePred(pb.d, pb.k, pb.bias)

initsol{L,T}(pb::MultivariateRegression{L,T}) = zeros(T, pb.k, has_bias(pb) ? pb.d + 1: pb.d)

# specific problems

linearreg{T}(X::StridedMatrix{T}, Y::StridedMatrix{T}; bias::Real=0.0) =
    MultivariateRegression(SumSqrLoss(), X, Y, size(Y,1); bias=bias)

mlogisticreg{T}(X::StridedMatrix{T}, Y::StridedVector{Int}, k::Int; bias::Real=0.0) =
    MultivariateRegression(MultiLogisticLoss(), X, Y, k; bias=bias)

ext: jl


text:
# functions for projecting back and forth between x-y and n-t
# (normal-tangential) coordinate systems
# in 3D, x-y-z to n-t-b (normal-tangential-binormal)
"""
  This function populates the matrix P that project from x-y-z to n-t-b.
  Methods are available for 2D and 3D, determined from the AbstractParamType
  object.  This is a somewhat specialized routine in that it only works for
  vectors like the state vector of the Euler or Navier-Stokes equations,
  where the first and last equations are coordinate system invarient and
  only the middle 2 (in 2D) or 3 (in 3D) equations need to be rotated.
  Specialized multiplication routines are provied as projectToXY and
  projectToNT.

  Inputs:

    params:  An AbstractParamType{Tdim}
    nrm: the normal direction in x-y coordinate system.  Must be a unit vector

  Inputs/Outputs

    P:  matrix to be populated with projection matrix

  Aliasing restrictions: none
"""
function getProjectionMatrix(params::AbstractParamType{2}, nrm::AbstractVector,
                             P::AbstractMatrix)
# test: P.' = inv(P)

  t1, t2 = getOrthogonalVector(params, nrm)

  # zero out P
  for i=1:length(P)
    P[i] = 0
  end
  P[1,1] = 1
  P[4,4] = 1
  P[2,2] = nrm[1]
  P[2,3] = nrm[2]
  P[3,2] = t1
  P[3,3] = t2

  return nothing
end

function getProjectionMatrix(params::AbstractParamType{3}, nrm::AbstractVector,
                             P::AbstractMatrix)
# test: P.' = inv(P)
  n1 = nrm[1]; n2 = nrm[2]; n3 = nrm[3]
  t1, t2, t3 = getOrthogonalVector(params, nrm)
  b1, b2, b3 = getBinormalVector(n1, n2, n3, t1, t2, t3)

  for i=1:length(P)
    P[i] = 0
  end

  P[1,1] = 1
  P[5,5] = 1
  P[2,2] = n1
  P[2,3] = n2
  P[2,4] = n3
  P[3,2] = t1
  P[3,3] = t2
  P[3,4] = t3
  P[4,2] = b1
  P[4,3] = b2
  P[4,4] = b3

  return nothing
end

"""
  This function generates a unit vector orthogonal to the input vector nrm and
  returns its components.  Methods are availble for 2D and 3D

  Inputs:

    params: an AbstractParamType, used to dispatch to the 2D or 3D method
    nrm: the input vector

  Outputs:

    t1, t2, (and t3 in 3D): the components of the unit vector orthogonal to
                            nrm

  Aliasing restrictions: none
"""
function getOrthogonalVector(params::AbstractParamType{2}, nrm::AbstractVector)
# this code uses a clever trick to deal with the cases where one of the
# components of nrm is zero in a branch free manner

  # add this factor to the denominators to avoid NaNs
  # the round(nrm[2]) below ensure only the result with
  # denominator > 0.5 gets chosen, so adding 1e-50 doesn't
  # affect the results
  add_fac = 1e-50

  # use the fact that dot(nrm, tangent) = 0 to produce a single equation
  # in 2 unknowns
  # pick 2 solutions, one with the first component equal to 2, one with the
  # second component equal to 1.  Below the code figures out which one to use
  v1 = 1.0
  v2 = -nrm[1]/(nrm[2] + add_fac)
  w1 = -nrm[2]/(nrm[1] + add_fac)
  w2 = 1.0

  # if nrm[2] > 0.5, use v, otherwise use w
  fac = round(absvalue(nrm[2]))
  t1 = fac*v1 + (1-fac)*w1
  t2 = fac*v2 + (1-fac)*w2

  len = sqrt(t1*t1 + t2*t2)

  t1 = t1/len
  t2 = t2/len

  return t1, t2
end

function getOrthogonalVector(params::AbstractParamType{3}, nrm::AbstractVector)
# this code uses a clever trick to deal with the cases where one of the
# components of nrm is zero in a branch free manner

  # unpack normal vector
  n1 = nrm[1]
  n2 = nrm[2]
  n3 = nrm[3]

  # add this factor to the denominators to avoid NaNs
  # the round(nrm[3]) below ensure only the result with
  # denominator > 0.5 gets chosen, so adding 1e-50 doesn't
  # affect the results
  add_fac = 1e-50

  # compute the three possible vectors
  v1 = 1.0
  v2 = 1.0
  v3 = -(n1 + n2)/(n3 + add_fac)

  w1 = 1.0
  w2 = -(n1 + n3)/(n2 + add_fac)
  w3 = 1.0

  x1 = -(n2 + n3)/(n1 + add_fac)
  x2 = 1.0
  x3 = 1.0

  # pick either w or x
  fac = round(absvalue(nrm[3]))

  z1 = fac*x1 + (1-fac)*w1
  z2 = fac*x2 + (1-fac)*w2
  z3 = fac*x3 + (1-fac)*w3

  # pick either v or z
  fac = round(absvalue(nrm[1]))
  t1 = fac*v1 + (1-fac)*z1
  t2 = fac*v2 + (1-fac)*z2
  t3 = fac*v3 + (1-fac)*z3

  len = sqrt(t1*t1 + t2*t2 + t3*t3)

  t1 = t1/len
  t2 = t2/len
  t3 = t3/len

  return t1, t2, t3
end

"""
  This function computes a vector normal to the 2 supplied vectors and
  returns the components.

  Inputs:

    n1, n2, n3: the components of the first vector
    t1, t2, t3: the components of the second vector

  Outputs:

    b1, b2, b3: the components of the binormal vector

  Aliasing restrictions: none
"""
@inline function getBinormalVector(n1::Number, n2::Number, n3::Number,
                           t1::Number, t2::Number, t3::Number)

  b1 = n2*t3 - n3*t2
  b2 = -(n1*t3 - n3*t1)
  b3 = n1*t2 - n2*t1

  return b1, b2, b3
end

"""
  Calculates the length of a vector, using a manually unrolled loop.
  Methods are available for 2D and 3D.

  Inputs:

    params: an AbstractParamType{Tdim}, used to dispatch to the right method
    nrm: the vector to calculate the length of.  Must have length 2 in 2D
         and 3 in 3D

  Outputs:

    length of nrm
"""
@inline function calcLength(params::AbstractParamType{2}, nrm::AbstractVector)
  return sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2])
end

@inline function calcLength(params::AbstractParamType{3}, nrm::AbstractVector)
  return sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3])
end

@inline function calcLength_rev(params::AbstractParamType{2}, nrm::AbstractVector,
                            nrm_bar::AbstractVector, len_bar::Number)

  len = sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2])
  fac = len_bar/len
  nrm_bar[1] += nrm[1]*fac
  nrm_bar[2] += nrm[2]*fac

  return nothing
end

@inline function calcLength_rev(params::AbstractParamType{3}, nrm::AbstractVector,
                            nrm_bar::AbstractVector, len_bar::Number)

  len = sqrt(nrm[1]*nrm[1] + nrm[2]*nrm[2] + nrm[3]*nrm[3])
  fac = len_bar/len
  nrm_bar[1] += nrm[1]*fac
  nrm_bar[2] += nrm[2]*fac
  nrm_bar[3] += nrm[3]*fac

  return nothing
end

function normalize_vec(params::AbstractParamType{T}, nrm::AbstractVector) where {T}
  len = calcLength(params, nrm)
  for d=1:T
    nrm[d] /= len
  end
end

"""
  This function projects a vector x from x-y coordinates to normal-tangential
  (n-t) coordinates, where A is a projection matrix from x-y to n-t, obtained
  from getProjectionMatrix.  This function is a specialized matrix-vector
  multiplication routine and only works for matrices with the particular
  sparsity pattern created by getProjectionMatrix.
  Methods are available for 2D and 3D

  Inputs:

    params: an AbstractParamType{{Tdim} used to dispatch to the 2D or 3D
            method
    P: the projection matrix
    x: vector to be projected

  Inputs/Outputs:

    b: the result of the projection

  Aliasing restrictions: x and b cannot alias
"""

function projectToNT(params::AbstractParamType{2}, P::AbstractMatrix,
                     x::AbstractVector, b::AbstractVector)
# regular multiplication
  m, n = size(P)

  # first loop: overwrite b
  b[1] = x[1]*P[1, 1]
  b[2] = 0
  b[3] = 0
  b[4] = 0

  b[2] += P[2,2]*x[2]
  b[3] += P[3,2]*x[2]

  b[2] += P[2,3]*x[3]
  b[3] += P[3,3]*x[3]

  b[4] += P[4,4]*x[4]

  return nothing
end

function projectToNT(params::AbstractParamType{3}, P::AbstractMatrix,
                     x::AbstractVector, b::AbstractVector)
# regular multiplication

  m, n = size(P)

  # first loop: overwrite b
  b[1] = x[1]*P[1, 1]
  b[2] = 0
  b[3] = 0
  b[4] = 0
  b[5] = 0

  b[2] += P[2,2]*x[2]
  b[3] += P[3,2]*x[2]
  b[4] += P[4,2]*x[2]

  b[2] += P[2,3]*x[3]
  b[3] += P[3,3]*x[3]
  b[4] += P[4,3]*x[3]

  b[2] += P[2,4]*x[4]
  b[3] += P[3,4]*x[4]
  b[4] += P[4,4]*x[4]

  b[5] += P[5,5]*x[5]

  return nothing
end

"""
  This function is similar to projectToNT, except it project from n-t
  to x-y.  Note that P is still the projection matrix from getProjectionMatrix
  that projects from x-y to n-t.
"""
function projectToXY(params::AbstractParamType{2}, P::AbstractMatrix,
                     x::AbstractVector, b::AbstractVector)
# transpose multiplication

  m, n = size(P)

  b[1] = P[1,1]*x[1]

  b[2] = 0
  b[2] += P[2,2]*x[2]
  b[2] += P[3,2]*x[3]

  b[3] = 0
  b[3] += P[2,3]*x[2]
  b[3] += P[3,3]*x[3]

  b[4] = 0
  b[4] += P[4,4]*x[4]

  return nothing
end

function projectToXY(params::AbstractParamType{3}, P::AbstractMatrix,
                     x::AbstractVector, b::AbstractVector)
# transposed multiplication

  m, n = size(P)

  b[1] = P[1,1]*x[1]

  b[2] = 0
  b[2] += P[2,2]*x[2]
  b[2] += P[3,2]*x[3]
  b[2] += P[4,2]*x[4]

  b[3] = 0
  b[3] += P[2,3]*x[2]
  b[3] += P[3,3]*x[3]
  b[3] += P[4,3]*x[4]

  b[4] = 0
  b[4] += P[2,4]*x[2]
  b[4] += P[3,4]*x[3]
  b[4] += P[4,4]*x[4]

  b[5] = 0
  b[5] += P[5,5]*x[5]

  return nothing
end


text:
# cudnn wrapper

"""

`conv_f(x,w;kwargs...)` executes convolutions or cross-correlations
using filters specified with `w` over tensor `x`.
Here is a description of all available keyword arguments:
* padding: the number of extra zeros implicitly concatenated at the start and at the end of each dimension. Default=floor((filterSize-1)/2) which preserves the input size when filterSize is odd and stride=1.
* stride: the number of elements to slide to reach the next filtering window. Default=1.
* upscale: upscale factor for each dimension. Default=1.
* mode: 0 for convolution and 1 for cross-correlation.  Default=0.
* alpha: can be used to scale the result. Default=1.
* algo: specifies which convolution algorithm shoud be used to compute the results. Default=0. See the CUDNN User Guide for details.
* workSpace: data pointer to GPU memory to a workspace needed to able to execute the specified algorithm. Default=C_NULL.
* workSpaceSizeInBytes: the size in bytes of the provided workSpace. Default=0.
* handle: handle to a previously created cuDNN context. Default=Dex allocated context.

"""
function conv_f{T}(x::DexArray{T},w::DexArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T),
                  algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    y = similar(x, cdims(x,w;o...))
    @cuda(cudnn, cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,FD(w),w,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(y),y)
    return y
end

"""

`conv_bx(x,w,dy;kwargs...)` This function computes the convolution gradient with
respect to the output tensor, returning results in gradDesc.

"""
function conv_bx{T}(x::DexArray{T},w::DexArray{T},dy::DexArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dx = similar(x)
    @cuda(cudnn,cudnnConvolutionBackwardData,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),FD(w),w,TD(dy),dy,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(dx),dx)
    return dx
end

"""

`conv_bw(x,w,dy;kwargs...)` This function computes the convolution gradient with
 respect to filter coefficients, returning results in gradDesc.

"""
function conv_bw{T}(x::DexArray{T},w::DexArray{T},dy::DexArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dw = similar(w)
    @cuda(cudnn,cudnnConvolutionBackwardFilter,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,TD(dy),dy,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),FD(dw),dw)
    return dw
end

"""

`pool_f(x;kwargs...)` computes pooling of input values (i.e., the
maximum or average of several adjacent values) to produce an output
with smaller height and/or width.

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

"""

function pool_f{T}(x::DexArray{T}; handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    y = similar(x, pdims(x; o...))
    @cuda(cudnn, cudnnPoolingForward,
          (Cptr, Cptr,      Ptr{T},    Cptr,Ptr{T},Ptr{T},   Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(x),x,    Ref(beta),TD(y),y)
    return y
end

function pool_b{T}(x::DexArray{T},y::DexArray{T},dy::DexArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    dx = similar(x)
    @cuda(cudnn,cudnnPoolingBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(y),y,TD(dy),dy,TD(x),x,Ref(beta),TD(dx),dx)
    return dx
end

"""

`dropout_f(x;kwargs...)` computes the forward dropout operation over input values
to produce an output.

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

"""

function dropout_f{T}(x::DexArray, dropout;
                      handle=cudnnhandle)
    y = similar(x)
    @cuda(cudnn, cudnnDropoutForward,
          (Cptr,Cptr,        Cptr,Ptr{T}
          handle,DD(dropout),TD(x)),x,

end

# cudnn datatype
# an enumerated type indicating the data type to which a tensor descriptor
# or filter descriptor refers.
DT(::DexArray{Float32})=UInt32(0)
DT(::DexArray{Float64})=UInt32(1)
DT(::DexArray{Float16})=UInt32(2)

# accquire dimension
function cdims{T,N}(x::DexArray{T,N},w::DexArray{T,N}; padding=0, stride=1, o...)
    ntuple(N) do i
        if i < N-1
            pi = if isa(padding,Number); padding; else padding[i]; end
            si = if isa(stride,Number); stride; else stride[i]; end
            1 + div(size(x,i) - size(w,i) + 2*pi, si)
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

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

# cudnn descriptors

type TD; ptr
    function TD(a::DexArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateTensorDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]    # why dim reversed?
        st = [Cint(stride(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetTensorNdDescriptor,
              (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint}),
              d[1], DT(a), n, sz, st)
        td = new(d[1])
        finalizer(td, x->@cuda(cudnn,cudnnDestroyTensorDescriptor,(Cptr,),x.ptr))
        return td
    end
end

type FD; ptr
    function FD(a::DexArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateFilterDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetFilterNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
              d[1], DT(a), 0,     n,   sz)
        fd = new(d[1])
        finalizer(fd, x->@cuda(cudnn,cudnnDestroyFilterDescriptor,(Cptr,),x.ptr))
        return fd
    end
end

type CD; ptr
    function CD(x::DexArray,w::DexArray; padding=0, stride=1, upscale=1, mode=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
              (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
              d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        cd = new(d[1])
        finalizer(cd, x->@cuda(cudnn,cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr))
        return cd
    end
end

type PD; ptr
    function PD(x::DexArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreatePoolingDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetPoolingNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
              d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        pd = new(d[1])
        finalizer(pd, x->@cuda(cudnn,cudnnDestroyPoolingDescriptor,(Cptr,),x.ptr))
        return pd
    end
end

type DD; ptr
    function DD(dropout;handle=cudnnhandle)
        d = Cptr[0]
        @cuda(cudnn, cudnnCreateDropoutDescriptor, (Ptr{Cptr},),d)
        statessize_p = Csize_t[0]
        @cuda(cudnn, cudnnDropoutGetStatesSize, (Cptr,Ptr{Csize_t}), handle,statessize_p)
        statessize = statessize_p[1]
        states = DexArray(Int8, Int(statessize))
        @cude(cudnn, cudnnSetDropoutDescriptor,
              (Cptr,Cptr, Cfloat,Cptr,  Cint,Culonglong),
              d[1],handle,Cfloat(dropout),state,statessize, 0)
        dd = new(d[1])
        finalizer(dd, x->@cuda(cudnn, cudnnDestroyDropoutDescriptor, (Cptr,),x.ptr)
        return dd
    end
end

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

function cdsize(w, nd)
    if isa(w,Integer)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch())
    end
end

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


text:
#=============================================
Validation methods for dictionaries of options
==============================================#

available_keywords = Set{Symbol}([])

"""
    check_aliases!(dict, dict_copy, aliases)

This function has several purposes:

- translate aliases to the option that is used internally
- check that not several aliases were used at the same time

INPUT:

- `dict` -- a dictionary of options
- `dict_copy` -- a copy of the dictionary of options for internal names
- `aliases` -- option aliases; the first name is the one we use internally
"""
function check_aliases!(dict::Dict{Symbol,Any}, dict_copy::Dict{Symbol,Any}, aliases::Vector{Symbol})
    # find aliases and check consistency in case several aliases have been used
    print_warning::Bool = false
    for alias in aliases
        # add alias to global keyword set
        push!(available_keywords, alias)

        if haskey(dict, alias)
            # alias was used
            if haskey(dict_copy, aliases[1])
                # several aliases were used
                if (dict_copy[aliases[1]] != dict[alias])
                    error("Several option aliases were used with different values for aliases $aliases.")
                else
                    print_warning = true
                end
            else
                dict_copy[aliases[1]] = dict[alias]
            end
        end
    end
    if print_warning
        warn("Several option aliases were used for aliases $aliases.")
    end
end

"""
    check_aliases_and_add_default_value!(dict, dict_copy, aliases, default_value, [modify_dict])

This function has several purposes:

- translate aliases to the option that is used internally
- check that not several aliases were used at the same time
- assign default values for undefined options

### Input

- `dict`          -- a dictionary of options
- `dict_copy`     -- a copy of the dictionary of options for internal names
- `aliases`       -- option aliases; the first name is the one we use internally
- `default_value` -- the default value for the option
- `modify_dict`   -- (optional, default: `false`) indicates if `dict` should be
                     modified
"""
function check_aliases_and_add_default_value!(dict::Dict{Symbol,Any}, dict_copy::Dict{Symbol,Any}, aliases::Vector{Symbol}, default_value::Any, modify_dict::Bool=false)
    check_aliases!(dict, dict_copy, aliases)

    if !haskey(dict_copy, aliases[1])
        # no alias and no value
        # TODO<notification> print message to user for each default option, depending on verbosity level
        if (modify_dict)
            dict[aliases[1]] = default_value
        end
        dict_copy[aliases[1]] = default_value
    end
end

"""
    check_valid_option_keywords(dict)

Validate that all input keywords are recognized.

### Input

- `dict` -- input dictionary
"""
function check_valid_option_keywords(dict)
    for kv_pair in dict
        key::Symbol = kv_pair.first
        if !in(key, available_keywords)
            error(get_unrecognized_key_message(key))
        end
    end
end

"""
    get_unrecognized_key_message(key)

Create an error message for an unrecognized option key.

### Input

- `key` -- unrecognized option key

### Output

The error message.
"""
function get_unrecognized_key_message(key)
    return "Unrecognized option '$key' found. See " *
        "`keys(Reachability.available_keywords.dict)` for all valid keywords."
end


text:
#sigmoid typedef - type definition for sigmoid-valued numbers.
#environment bits parameter.

#for now, put this here
__sigmoid_settings = Dict{Symbol, Any}()
getsetting(k) = haskey(__sigmoid_settings, k) ? __sigmoid_settings[k] : nothing

if (Int == Int32) || getsetting(:basebits) == 32
  const __BITS = 32
elseif getsetting(:basebits) == 16
  const __BITS = 16
elseif getsetting(:basebits) == 8
  const __BITS = 8
else #default to a 64-bit environment.
  const __BITS = 64
end

primitive type Sigmoid{N,ES,mode} <: AbstractFloat __BITS end

#these are deliberately made incompatible with the standard rounding modes types
#found in the julia std library.

const roundingmodes = [:guess,
  :ubit,
  :exact,
  :upper,
  :lower]

#set some type aliases.
Posit{N, ES} = Sigmoid{N, ES, :guess}
Vnum{N, ES} = Sigmoid{N, ES, :ubit}

#there's a couple of dummy types that we'll use for syntatical sugar purposes.
Exact{N,ES} = Sigmoid{N, ES, :EXACT}
ULP{N,ES} = Sigmoid{N, ES, :ULP}
#trampoline their constructor against the Vnum constructor.
(::Type{Exact{N,ES}}){N,ES}(n::Unsigned)::Vnum{N,ES} = iseven(n) ? Vnum{N,ES}(n) : throw(ArgumentError("Exact numbers must have an even int representation!"))
(::Type{ULP{N,ES}}){N,ES}(n::Unsigned)::Vnum{N,ES}   = isodd(n) ? Vnum{N,ES}(n) : throw(ArgumentError("ULP numbers must have an odd int representation!"))

struct Valid{N, ES} <: AbstractFloat
  lower::Vnum{N, ES}
  upper::Vnum{N, ES}
end

#create a right arrow function representation for the construction of valids.
→{N,ES}(lower::Vnum{N,ES}, upper::Vnum{N,ES}) = Valid{N,ES}(lower, upper)

#create a tile constructor that builds a tile from a vnum.
tile{N,ES}(x::Sigmoid{N,ES,:ubit}) = Valid{N,ES}(reinterpret(Vnum{N,ES},x), reinterpret(Vnum{N,ES},x))
tile{N,ES}(x::Sigmoid{N,ES,:exact}) = Valid{N,ES}(reinterpret(Vnum{N,ES},x), reinterpret(Vnum{N,ES},x))

export Sigmoid, Posit, Vnum, Valid, Exact, ULP, →

#sigmoid numbers don't natively have NaN, so NaNs should all be noisy.
struct NaNError <: Exception
  operand::Function
  parameters::Array{Any,1}
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using QuantLib.Time, QuantLib.Math

const BASIS_POINT = 0.0001

## TYPES ##
## Coupon pricers
type BlackIborCouponPricer{O <: OptionletVolatilityStructure} <: IborCouponPricer
  discount::Float64
  spreadLegValue::Float64
  accrual_period::Float64
  initialized::Bool
  capletVolatility::O

  # function call(::Type{BlackIborCouponPricer})
  #   new{OptionletVolatilityStructure}(0.0, 0.0, 0.0, false)
  # end

  # function call{O}(::Type{BlackIborCouponPricer}, capletVolatility::O)
  #   new{O}(0.0, 0.0, 0.0, false, capletVolatility)
  # end
end

BlackIborCouponPricer() = BlackIborCouponPricer{NullOptionletVolatilityStructure}(0.0, 0.0, 0.0, false, NullOptionletVolatilityStructure())
BlackIborCouponPricer(capVol::OptionletVolatilityStructure) = BlackIborCouponPricer{typeof(capVol)}(0.0, 0.0, 0.0, false, capVol)

type IborCoupon{DC <: DayCount, X <: InterestRateIndex, ICP <: IborCouponPricer} <: Coupon
  couponMixin::CouponMixin{DC}
  paymentDate::Date
  nominal::Float64
  fixingDate::Date
  fixingValueDate::Date
  fixingEndDate::Date
  fixingDays::Int
  iborIndex::X
  gearing::Float64
  spread::Float64
  isInArrears::Bool
  spanningTime::Float64
  pricer::ICP
end

function IborCoupon{X <: InterestRateIndex, DC <: DayCount, ICP <: IborCouponPricer}(paymentDate::Date, nominal::Float64, startDate::Date, endDate::Date, fixingDays::Int, iborIndex::X,
                    gearing::Float64, spread::Float64, refPeriodStart::Date, refPeriodEnd::Date, dc::DC, isInArrears::Bool, pricer::ICP)
  # TODO check if right fixing days
  _fixing_date = isInArrears ? fixing_date(iborIndex, endDate) : fixing_date(iborIndex, startDate)
  fixing_cal = iborIndex.fixingCalendar
  idx_fixing_days = iborIndex.fixingDays
  fixing_val_date = advance(Base.Dates.Day(idx_fixing_days), fixing_cal, _fixing_date, iborIndex.convention)

  if isInArrears
    fixing_end_date = maturity_date(iborIndex, fixing_val_date)
  else
    next_fixing = advance(-Base.Dates.Day(fixingDays), fixing_cal, endDate, iborIndex.convention)
    fixing_end_date = advance(Base.Dates.Day(idx_fixing_days), fixing_cal, next_fixing, iborIndex.convention)
  end

  spanning_time = year_fraction(iborIndex.dc, fixing_val_date, fixing_end_date)

  ## TODO ensure positive (> 0) spanning_time

  return IborCoupon{DC, X, ICP}(CouponMixin{DC}(startDate, endDate, refPeriodStart, refPeriodEnd, dc, -1.0), paymentDate, nominal, _fixing_date, fixing_val_date, fixing_end_date, fixingDays, iborIndex, gearing, spread,
                    isInArrears, spanning_time, pricer)
end

# Coupon methods
amount(coup::IborCoupon) = calc_rate(coup) * accrual_period(coup) * coup.nominal

get_pay_dates(coups::Vector{Union{IborCoupon, SimpleCashFlow}}) = Date[date(coup) for coup in filter(check_coupon, coups)]
get_reset_dates(coups::Vector{Union{IborCoupon, SimpleCashFlow}}) = Date[accrual_start_date(coup) for coup in filter(check_coupon, coups)]
get_gearings(coups::Vector{Union{IborCoupon, SimpleCashFlow}}) = Float64[coup.gearing for coup in filter(check_coupon, coups)]

function calc_rate(coup::IborCoupon)
  QuantLib.initialize!(coup.pricer, coup)
  # if coup.isInArrears
  #   println("sr ", sr)
  #   error("BREAK")
  # end
  return swaplet_rate(coup.pricer, coup)
end

function index_fixing(coupon::IborCoupon)
  today = settings.evaluation_date

  if coupon.fixingDate >= today
    return forecast_fixing(coupon.iborIndex, coupon.iborIndex.ts, coupon.fixingValueDate, coupon.fixingEndDate, coupon.spanningTime)
  end

  pastFix = get(coupon.iborIndex.pastFixings, coupon.fixingDate, -1.0)

  if pastFix == -1.0
    return forecast_fixing(coupon.iborIndex, coupon.iborIndex.ts, coupon.fixingValueDate, coupon.fixingEndDate, coupon.spanningTime)
  else
    return pastFix
  end

  # error("Fixing date on or before eval date")
end

function accrued_amount(coup::IborCoupon, settlement_date::Date)
  if settlement_date <= accrual_start_date(coup) || settlement_date > coup.paymentDate
    return 0.0
  end

  return coup.nominal * calc_rate(coup) * year_fraction(get_dc(coup), accrual_start_date(coup), min(settlement_date, accrual_end_date(coup)), ref_period_start(coup), ref_period_end(coup))
end

# Floating legs
type IborLeg <: Leg
  coupons::Vector{Union{IborCoupon, SimpleCashFlow}}

  function IborLeg{X <: InterestRateIndex, DC <: DayCount, C <: BusinessDayConvention, ICP <: IborCouponPricer}(schedule::Schedule, nominal::Float64, iborIndex::X, paymentDC::DC, paymentAdj::C,
                   fixingDays::Vector{Int} = fill(iborIndex.fixingDays, length(schedule.dates) - 1),
                   gearings::Vector{Float64} = ones(length(schedule.dates) - 1), spreads::Vector{Float64} = zeros(length(schedule.dates) - 1),
                   caps::Vector{Float64} = Vector{Float64}(), floors::Vector{Float64} = Vector{Float64}(), isInArrears::Bool = false,
                   isZero::Bool = false, pricer::ICP = BlackIborCouponPricer();
                   add_redemption::Bool = true, cap_vol::OptionletVolatilityStructure = NullOptionletVolatilityStructure())
    n = add_redemption ? length(schedule.dates) : length(schedule.dates) - 1
    if isa(pricer.capletVolatility, NullOptionletVolatilityStructure) && ~isa(cap_vol, NullOptionletVolatilityStructure)
      pricer = BlackIborCouponPricer(cap_vol)
    end
    coups = Vector{Union{IborCoupon, SimpleCashFlow}}(n)
    last_payment_date = adjust(schedule.cal, paymentAdj, schedule.dates[end])

    _start = ref_start = schedule.dates[1]
    _end = ref_end = schedule.dates[2]
    payment_date = adjust(schedule.cal, paymentAdj, _end)
    #TODO: setup payment adjustments and the like
    coups[1] = IborCoupon(payment_date, nominal, _start, _end, fixingDays[1], iborIndex, gearings[1], spreads[1], ref_start, ref_end,
                          paymentDC, isInArrears, pricer)

    # ref_date = _start = end_date = _end = Date()
    # build coupons
    count = 2
    ref_start = _start = _end
    ref_end = _end = count == length(schedule.dates) ? schedule.dates[end] : schedule.dates[count + 1]
    payment_date = adjust(schedule.cal, paymentAdj, _end)

    while _start < schedule.dates[end]
      @inbounds coups[count] = IborCoupon(payment_date, nominal, _start, _end, fixingDays[count], iborIndex, gearings[count], spreads[count], ref_start, ref_end,
                            paymentDC, isInArrears, pricer)

      count += 1
      ref_start = _start = _end
      ref_end = _end = count == length(schedule.dates) ? schedule.dates[end] : schedule.dates[count + 1]
      payment_date = adjust(schedule.cal, paymentAdj, _end)
    end

    # for i = 1:n
    #   ref_start = _start = schedule.dates[i]
    #   ref_end = _end = schedule.dates[i + 1]
    #   payment_date = isZero ? last_payment_date : adjust(schedule.cal, paymentAdj, _end)
    #
    #   # Just Floating Rate Coupons right now
    #   coups[i] = IborCoupon(payment_date, nominal, _start, _end, fixingDays[i], iborIndex, gearings[i], spreads[i], ref_start, ref_end,
    #                         paymentDC, isInArrears, pricer)
    # end

    if add_redemption
      @inbounds coups[end] = SimpleCashFlow(nominal, _end)
    end

    new(coups)
  end
end

## Pricer Methods ##
function initialize!(pricer::BlackIborCouponPricer, coup::IborCoupon)
  idx = coup.iborIndex
  yts = idx.ts

  payment_date = date(coup)
  if payment_date > yts.referenceDate
    pricer.discount = discount(yts, payment_date)
  else
    pricer.discount = 1.0
  end

  pricer.accrual_period = accrual_period(coup)

  pricer.spreadLegValue = coup.spread * pricer.accrual_period * pricer.discount

  return pricer
end

function swaplet_price(pricer::BlackIborCouponPricer, coup::IborCoupon)
  _swaplet_price = adjusted_fixing(pricer, coup) * pricer.accrual_period * pricer.discount
  # if coup.isInArrears
  #   println("gearing ", coup.gearing)
  #   error("BREAK")
  # end
  return coup.gearing * _swaplet_price + pricer.spreadLegValue
end

swaplet_rate(pricer::BlackIborCouponPricer, coup::IborCoupon) = swaplet_price(pricer, coup) / (pricer.accrual_period * pricer.discount)

function adjusted_fixing(pricer::BlackIborCouponPricer, coup::IborCoupon, fixing::Float64 = -1.0)
  if fixing == -1.0
    fixing = index_fixing(coup)
  end

  if !coup.isInArrears
    return fixing
  end

  d1 = coup.fixingDate
  ref_date = pricer.capletVolatility.referenceDate
  if d1 <= ref_date
    return fixing
  end

  # See Hull, 4th ed., page 550
  idx = coup.iborIndex
  d2 = value_date(idx, d1)
  d3 = maturity_date(idx, d2)
  tau = year_fraction(idx.dc, d2, d3)
  variance = black_variance(pricer.capletVolatility, d1, fixing)
  # println("fixing ", fixing)
  # println("d1 ", d1)
  # println("d2 ", d2)
  # println("d3 ", d3)
  # println("tau ", tau)
  # println("variance ", variance)

  adj = fixing * fixing * variance * tau / (1.0 + fixing * tau)
  # println("adj ", adj)
  # error("BREAK")
  return fixing + adj
end

function update_pricer!{O <: OptionletVolatilityStructure}(leg::IborLeg, opt::O)
  for coup in leg.coupons
    if isa(coup, IborCoupon)
      coup.pricer.capletVolatility = opt
    end
  end

  return leg
end


text:
ENV["GKSwstype"] = "100"

using Pkg; Pkg.instantiate()

using JLD2
using FourierFlows
using Plots
using Printf
using FFTW: irfft

filename = "./data/euler_longtime.jld2"
filename_diags = "./data/euler_longtime_diags.jld2"

normalize_diags_with_initial = false

withoutgif(path) = (length(path)>3 && path[end-3:end] == ".gif") ? path[1:end-4] : path

"""
    uniquepath(path)

Returns `path` with a number appended if `isfile(path)`, incremented until `path` does not exist.
"""
function uniquegifpath(path)
  n = 1
  if isfile(path)
    path = withoutgif(path) * "_$n.gif"
  end
  while isfile(path)
    n += 1
    path = withoutgif(path)[1:end-length("_$(n-1)")] * "_$n.gif"
  end

  return path
end

moviegif_filename = "./movies/euler_longtime.gif"
moviegif_filename = uniquegifpath(moviegif_filename)
moviemp4_filename = moviegif_filename[1:end-4]*".mp4"

# ## Visualizing the simulation

function plot_output(x, y, ζ, ψ, t, k₀, ν, t_final)
    tν = ν * k₀^2 * t

    p_ζ = heatmap(x, y, ζ',
             aspectratio = 1,
                       c = :balance,
                    #clim = (-2, 2),
                   xlims = (x[1], x[end]),
                   ylims = (y[1], y[end]),
                  xlabel = "x",
                  ylabel = "y",
                   title = "vorticity, tν = "*@sprintf("%.2f", tν),
              framestyle = :box)

    p_ψ = contourf(x, y, ψ',
               aspectratio = 1,
                    levels = 11,
                     xlims = (x[1], x[end]),
                     ylims = (y[1], y[end]),
                    xlabel = "x",
                    ylabel = "y",
                     title = "streamfunction",
                framestyle = :box)

    p_diags1 = plot(4, # this means "a plot with two series"
                   label = ["|u|₂ / |u|₂(t=0)"  "|ζ|₂ / |ζ|₂(t=0)" "|ζ|₄ / |ζ|₄(t=0)" "|∇ζ|₂ / |∇ζ|₂(t=0)"],
                  legend = :topright,
               linewidth = 2,
                   alpha = 0.7,
                  xlabel = "ν k₀² t",
                   xlims = (0, 1.01 * ν * k₀^2 * t_final),
                   # ylims = (0, 3),
                   # yscale = :log10
                   )

    p_diags2 = plot(2, # this means "a plot with two series"
                   label = ["∂²ψ/∂x∂y(0, 0)" "(∂²/∂x²-∂²/∂y²)ψ(0, 0)"],
                  legend = :topright,
               linewidth = 2,
                   alpha = 0.7,
                  xlabel = "ν k₀² t",
                   xlims = (0, 1.01 * ν * k₀^2 * t_final),
                   # yscale = :log10,
                   # ylims = (0, 3),
                   )

    p_diags3 = plot(3, # this means "a plot with two series"
                   label = ["max(|∂²ψ/∂x²|)" "max(|∂²ψ/∂y²|)" "max(|∂²ψ/∂x∂y|)"],
                  legend = :topright,
               linewidth = 2,
                   alpha = 0.7,
                  xlabel = "ν k₀² t",
                   xlims = (0, 1.01 * ν * k₀^2 * t_final),
                   # yscale = :log10,
                   # ylims = (0, 3),
                   )

    l = @layout [ Plots.grid(1, 2)
                  c{0.22h}
                  d{0.22h}
                  e{0.22h} ]

    p = plot(p_ζ, p_ψ, p_diags1, p_diags2, p_diags3, layout = l, size = (900, 1000))

    return p
end

diags = jldopen(filename_diags)

# now let's make a movie of the flow fields

file = jldopen(filename)

nx, ny = file["grid/nx"], file["grid/ny"]
Lx, Ly = file["grid/Lx"], file["grid/Ly"]

global ν = file["params/ν"]
global k₀ = file["params/k₀"]
global nsubs = file["params/nsubs"]

grid = TwoDGrid(nx, Lx, ny, Ly)

x, y = grid.x, grid.y

iterations = parse.(Int, keys(file["snapshots/t"]))
final_iteration = iterations[end]

ζh_initial = file["snapshots/zetah/0"]
t_final = file["snapshots/t/$final_iteration"]

ζh = zeros(Complex{Float64}, (grid.nkr, grid.nl))
ψh = zeros(Complex{Float64}, (grid.nkr, grid.nl))

@. ζh = ζh_initial
ζ = irfft(deepcopy(ζh), grid.nx)
@. ψh = @. -grid.invKrsq * ζh
ψ = irfft(deepcopy(ψh), grid.nx)

energyh = @. 1 / 2 * grid.invKrsq * abs2(ζh_initial)
energy_initial = 1 / (grid.Lx * grid.Ly) * FourierFlows.parsevalsum(energyh, grid)
U = sqrt(2*energy_initial)

global Re = U * grid.Lx / ν

global p = plot_output(x, y, ζ, ψ, 0.0, k₀, ν, t_final)

global total_iterations = length(iterations)

startwalltime = time()

global ∂²ψ∂x² = zeros(grid.nx, grid.ny)
global ∂²ψ∂y² = zeros(grid.nx, grid.ny)
global ∂²ψ∂x∂y = zeros(grid.nx, grid.ny)

global ∂²ψ∂x²h = zeros(Complex{Float64}, grid.nkr, grid.nl)
global ∂²ψ∂y²h = zeros(Complex{Float64}, grid.nkr, grid.nl)
global ∂²ψ∂x∂yh = zeros(Complex{Float64}, grid.nkr, grid.nl)

anim = @animate for (i, iteration) in enumerate(iterations[1:end-1])

  if i%25== 0
    estimated_remaining_walltime = (time()-startwalltime)/60 / i * (total_iterations - i)
    log = @sprintf("estimated remaining walltime: %.2f min", estimated_remaining_walltime)
    println(log)
  end

  local t = file["snapshots/t/$iteration"]
  local tν = ν * k₀^2 * t

  local ζh = file["snapshots/zetah/$iteration"]

  @. ψh = -grid.invKrsq * ζh

  hypeh = @. (grid.kr^2 - grid.l^2) * ψh
  hype = irfft(deepcopy(hypeh), grid.nx)
  hype_00 = hype[Int(grid.nx/2), Int(grid.ny/2)]

  @. ∂²ψ∂x²h = -grid.kr^2 * ψh
  @. ∂²ψ∂y²h = -grid.l^2 * ψh
  @. ∂²ψ∂x∂yh = -grid.kr * grid.l * ψh

  local ∂²ψ∂x² = irfft(∂²ψ∂x²h, grid.nx)
  local ∂²ψ∂y² = irfft(∂²ψ∂y²h, grid.nx)
  local ∂²ψ∂x∂y = irfft(∂²ψ∂x∂yh, grid.nx)

  ∂²ψ∂x∂y_00 = ∂²ψ∂x∂y[Int(grid.nx/2), Int(grid.ny/2)]

  local ζ = irfft(deepcopy(ζh), grid.nx)
  local ψ = irfft(deepcopy(ψh), grid.nx)

  p[1][1][:z] = ζ
  p[1][:title] = "vorticity, tν = " * @sprintf("%.4f", tν)
  p[2][1][:z] = ψ
  p[2][:title] = "streamfunction"
  p[3][:title] = "Re = " * @sprintf("%.2f", Re)

  t_diags = diags["diags/energy/t"][(i-1)*nsubs+1:i*nsubs]
  t_diags = diags["diags/energy/t"][(i-1)*nsubs+1:i*nsubs]
  tν_diags = ν * k₀^2 * t_diags

  ΔE  = (diags["diags/energy/data"][(i-1)*nsubs+1:i*nsubs]diags["diags/energy/data"][1]).^(1/2)
  ΔΖ₂ = (diags["diags/enstrophyL2/data"][(i-1)*nsubs+1:i*nsubs]).^(1/2)
  ΔΖ₄ = (diags["diags/enstrophyL4/data"][(i-1)*nsubs+1:i*nsubs]).^(1/4)
  ΔP  = (diags["diags/palinstrophy/data"][(i-1)*nsubs+1:i*nsubs]).^(1/2)

  if normalize_diags_with_initial == true
    ΔE  /= (diags["diags/palinstrophy/data"][1])^(1/2)
    ΔΖ₂ /= (diags["diags/enstrophyL2/data"][1])^(1/2)
    ΔΖ₄ /= (diags["diags/enstrophyL4/data"][1])^(1/4)
    ΔP  /= (diags["diags/palinstrophy/data"][1])^(1/2)
  end

  push!(p[3][1], tν_diags, log10.(ΔE))
  push!(p[3][2], tν_diags, log10.(ΔΖ₂))
  push!(p[3][3], tν_diags, log10.(ΔΖ₄))
  push!(p[3][4], tν_diags, log10.(ΔP))
  push!(p[4][1], tν, ∂²ψ∂x∂y_00)
  push!(p[4][2], tν, hype_00)
  push!(p[5][1], tν, maximum(abs.(∂²ψ∂x²)))
  push!(p[5][2], tν, maximum(abs.(∂²ψ∂y²)))
  push!(p[5][3], tν, maximum(abs.(∂²ψ∂x∂y)))
end

gif(anim, moviegif_filename, fps=14)
mp4(anim, moviemp4_filename, fps=14)

max_stars_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_issues_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_forks_repo_name: navidcy/TwoDTurbulence-EulerLimit


text:
@testset "show" begin
    function show_str(d)
        io = IOBuffer()
        io_context = IOContext(io, :limit => true, :displaysize => (10, 20))
        show(io_context, MIME"text/plain"(), d)
        return String(take!(io))
    end

    if VERSION < v"1.6-"
        @test show_str(Indices{Int64}()) == "0-element Indices{Int64}"
        @test show_str(Indices{Int64}(1:3)) == "3-element Indices{Int64}\n 1\n 2\n 3"
        @test show_str(Indices{Int64}(1:100)) == "100-element Indices{Int64}\n 1\n 2\n 3\n ⋮\n 99\n 100"
        @test show_str(filterview(iseven, Indices{Int64}(1:100))) == "Greater than 6-element FilteredIndices{Int64,Indices{Int64},typeof(iseven)}\n 2\n 4\n 6\n ⋮\n 98\n 100"
        if Int === Int64
            @test show_str(Indices{Vector{Int64}}([collect(1:100)])) == "1-element Indices{Array{Int64,1}}\n [1, 2, 3, 4, 5, 6,…"
            @test show_str(Indices{Vector{Int64}}([collect(1:(100+i)) for i in 1:100])) == "100-element Indices{Array{Int64,1}}\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…\n ⋮\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…"
        else
            @test show_str(Indices{Vector{Int64}}([collect(1:100)])) == "1-element Indices{Array{Int64,1}}\n Int64[1, 2, 3, 4, …"
            @test show_str(Indices{Vector{Int64}}([collect(1:(100+i)) for i in 1:100])) == "100-element Indices{Array{Int64,1}}\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …\n ⋮\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …"
        end

        @test show_str(Dictionary{Int64,Int64}()) == "0-element Dictionary{Int64,Int64}"
        @test show_str(Dictionary{Int64,Int64}(1:3,101:103)) == "3-element Dictionary{Int64,Int64}\n 1 │ 101\n 2 │ 102\n 3 │ 103"
        @test show_str(Dictionary{Int64,Int64}(1:100,101:200)) == "100-element Dictionary{Int64,Int64}\n   1 │ 101\n   2 │ 102\n   3 │ 103\n   ⋮ │ ⋮\n  99 │ 199\n 100 │ 200"
        @test show_str(filterview(iseven, Dictionary{Int64,Int64}(1:100,101:200))) == "Greater than 6-element FilteredDictionary{Int64,Int64,Dictionary{Int64,Int64},typeof(iseven)}\n   2 │ 102\n   4 │ 104\n   6 │ 106\n   ⋮ │ ⋮\n  98 │ 198\n 100 │ 200"
        @test show_str(Dictionary{Int64,Int64}(collect(1:100), collect(101:200))) == "100-element Dictionary{Int64,Int64}\n   1 │ 101\n   2 │ 102\n   3 │ 103\n   ⋮ │ ⋮\n  99 │ 199\n 100 │ 200"
        if Int === Int64
            @test show_str(Dictionary{Vector{Int64},Vector{Int64}}([collect(1:100+i) for i in 1:100],[collect(101:200+i) for i in 1:100])) == "100-element Dictionary{Array{Int64,1},Array{Int64,1}}\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…\n        ⋮ │ ⋮\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…"
        else
            @test show_str(Dictionary{Vector{Int64},Vector{Int64}}([collect(1:100+i) for i in 1:100],[collect(101:200+i) for i in 1:100])) == "100-element Dictionary{Array{Int64,1},Array{Int64,1}}\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…\n        ⋮ │ ⋮\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…"
        end
    else
        @test show_str(Indices{Int64}()) == "0-element Indices{Int64}"
        @test show_str(Indices{Int64}(1:3)) == "3-element Indices{Int64}\n 1\n 2\n 3"
        @test show_str(Indices{Int64}(1:100)) == "100-element Indices{Int64}\n 1\n 2\n 3\n ⋮\n 99\n 100"
        @test show_str(filterview(iseven, Indices{Int64}(1:100))) == "Greater than 6-element FilteredIndices{Int64, Indices{Int64}, typeof(iseven)}\n 2\n 4\n 6\n ⋮\n 98\n 100"
        if Int === Int64
            @test show_str(Indices{Vector{Int64}}([collect(1:100)])) == "1-element Indices{Vector{Int64}}\n [1, 2, 3, 4, 5, 6,…"
            @test show_str(Indices{Vector{Int64}}([collect(1:(100+i)) for i in 1:100])) == "100-element Indices{Vector{Int64}}\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…\n ⋮\n [1, 2, 3, 4, 5, 6,…\n [1, 2, 3, 4, 5, 6,…"
        else
            @test show_str(Indices{Vector{Int64}}([collect(1:100)])) == "1-element Indices{Vector{Int64}}\n Int64[1, 2, 3, 4, …"
            @test show_str(Indices{Vector{Int64}}([collect(1:(100+i)) for i in 1:100])) == "100-element Indices{Vector{Int64}}\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …\n ⋮\n Int64[1, 2, 3, 4, …\n Int64[1, 2, 3, 4, …"
        end

        @test show_str(Dictionary{Int64,Int64}()) == "0-element Dictionary{Int64, Int64}"
        @test show_str(Dictionary{Int64,Int64}(1:3,101:103)) == "3-element Dictionary{Int64, Int64}\n 1 │ 101\n 2 │ 102\n 3 │ 103"
        @test show_str(Dictionary{Int64,Int64}(1:100,101:200)) == "100-element Dictionary{Int64, Int64}\n   1 │ 101\n   2 │ 102\n   3 │ 103\n   ⋮ │ ⋮\n  99 │ 199\n 100 │ 200"
        @test show_str(filterview(iseven, Dictionary{Int64,Int64}(1:100,101:200))) == "Greater than 6-element FilteredDictionary{Int64, Int64, Dictionary{Int64, Int64}, typeof(iseven)}\n   2 │ 102\n   4 │ 104\n   6 │ 106\n   ⋮ │ ⋮\n  98 │ 198\n 100 │ 200"
        @test show_str(Dictionary{Int64,Int64}(collect(1:100), collect(101:200))) == "100-element Dictionary{Int64, Int64}\n   1 │ 101\n   2 │ 102\n   3 │ 103\n   ⋮ │ ⋮\n  99 │ 199\n 100 │ 200"
        if Int === Int64
            @test show_str(Dictionary{Vector{Int64},Vector{Int64}}([collect(1:100+i) for i in 1:100],[collect(101:200+i) for i in 1:100])) == "100-element Dictionary{Vector{Int64}, Vector{Int64}}\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…\n        ⋮ │ ⋮\n [1, 2, … │ [101, 1…\n [1, 2, … │ [101, 1…"
        else
            @test show_str(Dictionary{Vector{Int64},Vector{Int64}}([collect(1:100+i) for i in 1:100],[collect(101:200+i) for i in 1:100])) == "100-element Dictionary{Vector{Int64}, Vector{Int64}}\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…\n        ⋮ │ ⋮\n Int64[1… │ Int64[1…\n Int64[1… │ Int64[1…"
        end
    end

end


text:
using Roots

# Given
P = 170000
F = 813000
n = 15

# find i
i = (F / P)^(1/n) - 1

@assert isapprox(F, P * (1 + i)^n; atol=0.01)


text:
include("utils.jl")
include("types.jl")
include("context.jl")


text:
module Monitor

export monitor, attach_monitor
using Caching

const init_size = 1000000

type monitor{TL,TV}
    fetch_function::Function
    label_function::Function

    labels::Vector{TL}
    values::Vector{TV}
    unique_only::Bool
end

function monitor(TL::Type,TV::Type,fetch_function::Function,label_function::Function;unique_only=true)
    new = monitor(fetch_function,label_function,Vector{TL}(),Vector{TV}(),unique_only)
    sizehint!(new.labels,init_size)
    sizehint!(new.values,init_size)
end

function monitor(TV::Type,fetch_function::Function;unique_only=true)
    monitor(fetch_function,Monitor.counter(),Vector{Int}(),Vector{TV}(),unique_only)
end

function monitor(TV::Type;unique_only=true)
    monitor(()->nothing,Monitor.counter(),Vector{Int}(),Vector{TV}(),unique_only)
end

function attach_monitor{TL,TV}(m::monitor{TL,TV},original_function::Function)
    prev_calls = Vector{UInt}()
    sizehint!(prev_calls,init_size)
    function(args...)
        if m.unique_only && !(hashkeys((args...)) in prev_calls)
            push!(m.labels,m.label_function()::TL)
            push!(m.values,m.fetch_function()::TV)
            push!(prev_calls,hashkeys((args...)))
        end
        original_function(args...)
    end
end

function counter(init::Int=0)
    function()
        init::Int+=1
        return init
    end
end

function hashkeys(tup::Tuple)
    hash(map(tup) do xi
        hash(xi)
    end)
end

end

max_stars_repo_name: sebastianpech/Monitor

max_issues_repo_name: sebastianpech/Monitor

max_forks_repo_name: sebastianpech/Monitor


text:
@enum Suit ♣ ♦ ♥ ♠
@enum Palo Bastos Espadas Oros Copas

cadena( n::Int, s::Suit ) = string( n, " ", s )
cadena( n::Int, p::Palo ) = string( n, " de ", p )

println( cadena( 3, ♦) )
println( cadena( 7, Bastos ) )


text:
using CEnum

function LLVMInitializeAllTargetInfos()
    ccall((:LLVMInitializeAllTargetInfos, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeAllTargets()
    ccall((:LLVMInitializeAllTargets, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeAllTargetMCs()
    ccall((:LLVMInitializeAllTargetMCs, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeAllAsmPrinters()
    ccall((:LLVMInitializeAllAsmPrinters, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeAllAsmParsers()
    ccall((:LLVMInitializeAllAsmParsers, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeAllDisassemblers()
    ccall((:LLVMInitializeAllDisassemblers, libLLVMExtra), Cvoid, ())
end

function LLVMInitializeNativeTarget()
    ccall((:LLVMInitializeNativeTarget, libLLVMExtra), LLVMBool, ())
end

function LLVMInitializeNativeAsmParser()
    ccall((:LLVMInitializeNativeAsmParser, libLLVMExtra), LLVMBool, ())
end

function LLVMInitializeNativeAsmPrinter()
    ccall((:LLVMInitializeNativeAsmPrinter, libLLVMExtra), LLVMBool, ())
end

function LLVMInitializeNativeDisassembler()
    ccall((:LLVMInitializeNativeDisassembler, libLLVMExtra), LLVMBool, ())
end

@cenum(LLVMDebugEmissionKind,
    LLVMDebugEmissionKindNoDebug = 0,
    LLVMDebugEmissionKindFullDebug = 1,
    LLVMDebugEmissionKindLineTablesOnly = 2,
    LLVMDebugEmissionKindDebugDirectivesOnly = 3,
)

function LLVMAddBarrierNoopPass(PM)
    ccall((:LLVMAddBarrierNoopPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddDivRemPairsPass(PM)
    ccall((:LLVMAddDivRemPairsPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddLoopDistributePass(PM)
    ccall((:LLVMAddLoopDistributePass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddLoopFusePass(PM)
    ccall((:LLVMAddLoopFusePass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddLoopLoadEliminationPass(PM)
    ccall((:LLVMAddLoopLoadEliminationPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddLoadStoreVectorizerPass(PM)
    ccall((:LLVMAddLoadStoreVectorizerPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMAddVectorCombinePass(PM)
    ccall((:LLVMAddVectorCombinePass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

if version() < v"12"
function LLVMAddInstructionSimplifyPass(PM)
    ccall((:LLVMAddInstructionSimplifyPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end
end

mutable struct LLVMOpaquePass end

const LLVMPassRef = Ptr{LLVMOpaquePass}

function LLVMAddPass(PM, P)
    ccall((:LLVMAddPass, libLLVMExtra), Cvoid, (LLVMPassManagerRef, LLVMPassRef), PM, P)
end

# typedef LLVMBool ( * LLVMPassCallback ) ( void * Ref , void * Data )
const LLVMPassCallback = Ptr{Cvoid}

function LLVMCreateModulePass2(Name, Callback, Data)
    ccall((:LLVMCreateModulePass2, libLLVMExtra), LLVMPassRef, (Cstring, LLVMPassCallback, Ptr{Cvoid}), Name, Callback, Data)
end

function LLVMCreateFunctionPass2(Name, Callback, Data)
    ccall((:LLVMCreateFunctionPass2, libLLVMExtra), LLVMPassRef, (Cstring, LLVMPassCallback, Ptr{Cvoid}), Name, Callback, Data)
end

function LLVMGetDebugMDVersion()
    ccall((:LLVMGetDebugMDVersion, libLLVMExtra), Cuint, ())
end

function LLVMGetValueContext(V)
    ccall((:LLVMGetValueContext, libLLVMExtra), LLVMContextRef, (LLVMValueRef,), V)
end

function LLVMAddTargetLibraryInfoByTriple(T, PM)
    ccall((:LLVMAddTargetLibraryInfoByTriple, libLLVMExtra), Cvoid, (Cstring, LLVMPassManagerRef), T, PM)
end

function LLVMAddInternalizePassWithExportList(PM, ExportList, Length)
    ccall((:LLVMAddInternalizePassWithExportList, libLLVMExtra), Cvoid, (LLVMPassManagerRef, Ptr{Cstring}, Csize_t), PM, ExportList, Length)
end

function LLVMExtraAppendToUsed(Mod, Values, Count)
    ccall((:LLVMExtraAppendToUsed, libLLVMExtra), Cvoid, (LLVMModuleRef, Ptr{LLVMValueRef}, Csize_t), Mod, Values, Count)
end

function LLVMExtraAppendToCompilerUsed(Mod, Values, Count)
    ccall((:LLVMExtraAppendToCompilerUsed, libLLVMExtra), Cvoid, (LLVMModuleRef, Ptr{LLVMValueRef}, Csize_t), Mod, Values, Count)
end

function LLVMExtraAddGenericAnalysisPasses(PM)
    ccall((:LLVMExtraAddGenericAnalysisPasses, libLLVMExtra), Cvoid, (LLVMPassManagerRef,), PM)
end

function LLVMExtraSetInitializer(GlobalVar, ConstantVal)
    ccall((:LLVMExtraSetInitializer, libLLVMExtra), Cvoid, (LLVMValueRef, LLVMValueRef), GlobalVar, ConstantVal)
end

function LLVMExtraSetPersonalityFn(Fn, PersonalityFn)
    ccall((:LLVMExtraSetPersonalityFn, libLLVMExtra), Cvoid, (LLVMValueRef, LLVMValueRef), Fn, PersonalityFn)
end

function LLVMExtraDIScopeGetName(Scope, Len)
    ccall((:LLVMExtraDIScopeGetName, libLLVMExtra), Cstring, (LLVMMetadataRef, Ptr{Cuint}), Scope, Len)
end

function LLVMAddCFGSimplificationPass2(PM, BonusInstThreshold, ForwardSwitchCondToPhi,
                                       ConvertSwitchToLookupTable, NeedCanonicalLoop,
                                       HoistCommonInsts, SinkCommonInsts,
                                       SimplifyCondBranch, FoldTwoEntryPHINode)
    ccall((:LLVMAddCFGSimplificationPass2, libLLVMExtra), Cvoid,
          (LLVMPassManagerRef, Cint, LLVMBool, LLVMBool, LLVMBool, LLVMBool, LLVMBool, LLVMBool, LLVMBool),
          PM, BonusInstThreshold, ForwardSwitchCondToPhi, ConvertSwitchToLookupTable, NeedCanonicalLoop,
          HoistCommonInsts, SinkCommonInsts, SimplifyCondBranch, FoldTwoEntryPHINode)
end

# bug fixes

# TODO: upstream

function LLVMExtraDumpMetadata(MD)
    ccall((:LLVMExtraDumpMetadata, libLLVMExtra), Cvoid, (LLVMMetadataRef,), MD)
end

function LLVMExtraPrintMetadataToString(MD)
    ccall((:LLVMExtraPrintMetadataToString, libLLVMExtra), Cstring, (LLVMMetadataRef,), MD)
end

function LLVMExtraGetMDNodeNumOperands2(MD)
    ccall((:LLVMExtraGetMDNodeNumOperands2, libLLVMExtra), Cuint, (LLVMMetadataRef,), MD)
end

function LLVMExtraGetMDString2(MD, Len)
    ccall((:LLVMExtraGetMDString2, libLLVMExtra), Cstring, (LLVMMetadataRef, Ptr{Cuint}), MD, Len)
end

function LLVMExtraGetMDNodeOperands2(MD, Dest)
    ccall((:LLVMExtraGetMDNodeOperands2, libLLVMExtra), Cvoid, (LLVMMetadataRef, Ptr{LLVMMetadataRef}), MD, Dest)
end

function LLVMExtraGetNamedMetadataNumOperands2(NMD)
    ccall((:LLVMExtraGetNamedMetadataNumOperands2, libLLVMExtra), Cuint, (LLVMNamedMDNodeRef,), NMD)
end

function LLVMExtraGetNamedMetadataOperands2(NMD, Dest)
    ccall((:LLVMExtraGetNamedMetadataOperands2, libLLVMExtra), Cvoid, (LLVMNamedMDNodeRef, Ptr{LLVMMetadataRef}), NMD, Dest)
end

function LLVMExtraAddNamedMetadataOperand2(NMD, Val)
    ccall((:LLVMExtraAddNamedMetadataOperand2, libLLVMExtra), Cvoid, (LLVMNamedMDNodeRef, LLVMMetadataRef), NMD, Val)
end

if version().major == 12

function LLVMCreateTypeAttribute(C, KindID, type_ref)
    ccall((:LLVMCreateTypeAttribute, libLLVMExtra), LLVMAttributeRef, (LLVMContextRef, Cuint, LLVMTypeRef), C, KindID, type_ref)
end

function LLVMGetTypeAttributeValue(A)
    ccall((:LLVMGetTypeAttributeValue, libLLVMExtra), LLVMTypeRef, (LLVMAttributeRef,), A)
end

function LLVMIsTypeAttribute(A)
    ccall((:LLVMIsTypeAttribute, libLLVMExtra), LLVMBool, (LLVMAttributeRef,), A)
end

struct LLVMOrcCSymbolFlagsMapPair
    Name::LLVMOrcSymbolStringPoolEntryRef
    Flags::LLVMJITSymbolFlags
end

const LLVMOrcCSymbolFlagsMapPairs = Ptr{LLVMOrcCSymbolFlagsMapPair}

struct LLVMOrcCSymbolAliasMapEntry
    Name::LLVMOrcSymbolStringPoolEntryRef
    Flags::LLVMJITSymbolFlags
end

struct LLVMOrcCSymbolAliasMapPair
    Name::LLVMOrcSymbolStringPoolEntryRef
    Entry::LLVMOrcCSymbolAliasMapEntry
end

const LLVMOrcCSymbolAliasMapPairs = Ptr{LLVMOrcCSymbolAliasMapPair}

mutable struct LLVMOrcOpaqueMaterializationResponsibility end

const LLVMOrcMaterializationResponsibilityRef = Ptr{LLVMOrcOpaqueMaterializationResponsibility}

# typedef void ( * LLVMOrcMaterializationUnitMaterializeFunction ) ( void * Ctx , LLVMOrcMaterializationResponsibilityRef MR )
const LLVMOrcMaterializationUnitMaterializeFunction = Ptr{Cvoid}

# typedef void ( * LLVMOrcMaterializationUnitDiscardFunction ) ( void * Ctx , LLVMOrcJITDylibRef JD , LLVMOrcSymbolStringPoolEntryRef Symbol )
const LLVMOrcMaterializationUnitDiscardFunction = Ptr{Cvoid}

# typedef void ( * LLVMOrcMaterializationUnitDestroyFunction ) ( void * Ctx )
const LLVMOrcMaterializationUnitDestroyFunction = Ptr{Cvoid}

# typedef LLVMErrorRef ( * LLVMOrcGenericIRModuleOperationFunction ) ( void * Ctx , LLVMModuleRef M )
const LLVMOrcGenericIRModuleOperationFunction = Ptr{Cvoid}

mutable struct LLVMOrcOpaqueIRTransformLayer end

const LLVMOrcIRTransformLayerRef = Ptr{LLVMOrcOpaqueIRTransformLayer}

# typedef LLVMErrorRef ( * LLVMOrcIRTransformLayerTransformFunction ) ( void * Ctx , LLVMOrcThreadSafeModuleRef * ModInOut , LLVMOrcMaterializationResponsibilityRef MR )
const LLVMOrcIRTransformLayerTransformFunction = Ptr{Cvoid}

mutable struct LLVMOrcOpaqueIndirectStubsManager end

const LLVMOrcIndirectStubsManagerRef = Ptr{LLVMOrcOpaqueIndirectStubsManager}

mutable struct LLVMOrcOpaqueLazyCallThroughManager end

const LLVMOrcLazyCallThroughManagerRef = Ptr{LLVMOrcOpaqueLazyCallThroughManager}

function LLVMOrcCreateCustomMaterializationUnit(Name, Ctx, Syms, NumSyms, InitSym, Materialize, Discard, Destroy)
    ccall((:LLVMOrcCreateCustomMaterializationUnit, libLLVMExtra), LLVMOrcMaterializationUnitRef, (Cstring, Ptr{Cvoid}, LLVMOrcCSymbolFlagsMapPairs, Csize_t, LLVMOrcSymbolStringPoolEntryRef, LLVMOrcMaterializationUnitMaterializeFunction, LLVMOrcMaterializationUnitDiscardFunction, LLVMOrcMaterializationUnitDestroyFunction), Name, Ctx, Syms, NumSyms, InitSym, Materialize, Discard, Destroy)
end

function LLVMOrcLazyReexports(LCTM, ISM, SourceRef, CallableAliases, NumPairs)
    ccall((:LLVMOrcLazyReexports, libLLVMExtra), LLVMOrcMaterializationUnitRef, (LLVMOrcLazyCallThroughManagerRef, LLVMOrcIndirectStubsManagerRef, LLVMOrcJITDylibRef, LLVMOrcCSymbolAliasMapPairs, Csize_t), LCTM, ISM, SourceRef, CallableAliases, NumPairs)
end

function LLVMOrcMaterializationResponsibilityGetTargetDylib(MR)
    ccall((:LLVMOrcMaterializationResponsibilityGetTargetDylib, libLLVMExtra), LLVMOrcJITDylibRef, (LLVMOrcMaterializationResponsibilityRef,), MR)
end

function LLVMOrcMaterializationResponsibilityGetExecutionSession(MR)
    ccall((:LLVMOrcMaterializationResponsibilityGetExecutionSession, libLLVMExtra), LLVMOrcExecutionSessionRef, (LLVMOrcMaterializationResponsibilityRef,), MR)
end

function LLVMOrcMaterializationResponsibilityGetInitializerSymbol(MR)
    ccall((:LLVMOrcMaterializationResponsibilityGetInitializerSymbol, libLLVMExtra), LLVMOrcSymbolStringPoolEntryRef, (LLVMOrcMaterializationResponsibilityRef,), MR)
end

function LLVMOrcMaterializationResponsibilityGetRequestedSymbols(MR, NumSymbols)
    ccall((:LLVMOrcMaterializationResponsibilityGetRequestedSymbols, libLLVMExtra), Ptr{LLVMOrcSymbolStringPoolEntryRef}, (LLVMOrcMaterializationResponsibilityRef, Ptr{Csize_t}), MR, NumSymbols)
end

function LLVMOrcDisposeSymbols(Symbols)
    ccall((:LLVMOrcDisposeSymbols, libLLVMExtra), Cvoid, (Ptr{LLVMOrcSymbolStringPoolEntryRef},), Symbols)
end

function LLVMOrcMaterializationResponsibilityNotifyResolved(MR, Symbols, NumPairs)
    ccall((:LLVMOrcMaterializationResponsibilityNotifyResolved, libLLVMExtra), LLVMErrorRef, (LLVMOrcMaterializationResponsibilityRef, LLVMOrcCSymbolMapPairs, Csize_t), MR, Symbols, NumPairs)
end

function LLVMOrcMaterializationResponsibilityNotifyEmitted(MR)
    ccall((:LLVMOrcMaterializationResponsibilityNotifyEmitted, libLLVMExtra), LLVMErrorRef, (LLVMOrcMaterializationResponsibilityRef,), MR)
end

function LLVMOrcMaterializationResponsibilityFailMaterialization(MR)
    ccall((:LLVMOrcMaterializationResponsibilityFailMaterialization, libLLVMExtra), Cvoid, (LLVMOrcMaterializationResponsibilityRef,), MR)
end

function LLVMOrcThreadSafeModuleWithModuleDo(TSM, F, Ctx)
    ccall((:LLVMOrcThreadSafeModuleWithModuleDo, libLLVMExtra), LLVMErrorRef, (LLVMOrcThreadSafeModuleRef, LLVMOrcGenericIRModuleOperationFunction, Ptr{Cvoid}), TSM, F, Ctx)
end

function LLVMOrcIRTransformLayerEmit(IRTransformLayer, MR, TSM)
    ccall((:LLVMOrcIRTransformLayerEmit, libLLVMExtra), Cvoid, (LLVMOrcIRTransformLayerRef, LLVMOrcMaterializationResponsibilityRef, LLVMOrcThreadSafeModuleRef), IRTransformLayer, MR, TSM)
end

function LLVMOrcCreateLocalIndirectStubsManager(TargetTriple)
    ccall((:LLVMOrcCreateLocalIndirectStubsManager, libLLVMExtra), LLVMOrcIndirectStubsManagerRef, (Cstring,), TargetTriple)
end

function LLVMOrcDisposeIndirectStubsManager(ISM)
    ccall((:LLVMOrcDisposeIndirectStubsManager, libLLVMExtra), Cvoid, (LLVMOrcIndirectStubsManagerRef,), ISM)
end

function LLVMOrcCreateLocalLazyCallThroughManager(TargetTriple, ES, ErrorHandlerAddr, Result)
    ccall((:LLVMOrcCreateLocalLazyCallThroughManager, libLLVMExtra), LLVMErrorRef, (Cstring, LLVMOrcExecutionSessionRef, LLVMOrcJITTargetAddress, Ptr{LLVMOrcLazyCallThroughManagerRef}), TargetTriple, ES, ErrorHandlerAddr, Result)
end

function LLVMOrcDisposeLazyCallThroughManager(LCM)
    ccall((:LLVMOrcDisposeLazyCallThroughManager, libLLVMExtra), Cvoid, (LLVMOrcLazyCallThroughManagerRef,), LCM)
end

function LLVMOrcLLJITGetIRTransformLayer(J)
    ccall((:LLVMOrcLLJITGetIRTransformLayer, libLLVMExtra), LLVMOrcIRTransformLayerRef, (LLVMOrcLLJITRef,), J)
end

function LLVMOrcLLJITApplyDataLayout(J, Mod)
    ccall((:LLVMOrcLLJITApplyDataLayout, libLLVMExtra), LLVMErrorRef, (LLVMOrcLLJITRef, LLVMModuleRef), J, Mod)
end

end # version

@cenum LLVMCloneFunctionChangeType::UInt32 begin
    LLVMCloneFunctionChangeTypeLocalChangesOnly = 0
    LLVMCloneFunctionChangeTypeGlobalChanges = 1
    LLVMCloneFunctionChangeTypeDifferentModule = 2
    LLVMCloneFunctionChangeTypeClonedModule = 3
end

function LLVMCloneFunctionInto(NewFunc, OldFunc, ValueMap, ValueMapElements, Changes, NameSuffix, TypeMapper, TypeMapperData, Materializer, MaterializerData)
    ccall((:LLVMCloneFunctionInto, libLLVMExtra), Cvoid, (LLVMValueRef, LLVMValueRef, Ptr{LLVMValueRef}, Cuint, LLVMCloneFunctionChangeType, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), NewFunc, OldFunc, ValueMap, ValueMapElements, Changes, NameSuffix, TypeMapper, TypeMapperData, Materializer, MaterializerData)
end

function LLVMFunctionDeleteBody(Func)
    ccall((:LLVMFunctionDeleteBody, libLLVMExtra), Cvoid, (LLVMValueRef,), Func)
end

function LLVMDestroyConstant(Const)
    ccall((:LLVMDestroyConstant, libLLVMExtra), Cvoid, (LLVMValueRef,), Const)
end


text:
"""
HistoryReplayer uses the generative PMDP interface to reaplay requests in history h.

Usage:

"""
struct HistoryReplayer <: PMDP{State, Action}
    pp::PMDPProblem
    empty_product::Product
    empty_product_id::Int64
    h::AbstractVector{<:NamedTuple}
    request_timesteps::Array{Int64}

    function HistoryReplayer(m::PMDP, h::AbstractVector{<:NamedTuple})
        return new(m.pp, empty_product(m), empty_product_id(m), h, collect([s.t for s in h[:s]]))
    end
end

POMDPs.initialstate(hr::HistoryReplayer) = Deterministic(hr.h[1].s)

"""
Returns next requested product. If in given timestep one of the prodcuts has selling period end, update the product request probs.
"""
function sample_request(hr::HistoryReplayer, t::Timestep, rng::AbstractRNG)::Int64
    history_index = findfirst(x->x==t, hr.request_timesteps)
    if history_index === nothing
        iₚ = empty_product_id(hr)
    else
        iₚ = hr.h[history_index].s.iₚ
    end
    return iₚ
end

"""
Sample user budget Budget for product requested in state s.
"""
function sample_customer_budget(hr::HistoryReplayer, s::State, rng::AbstractRNG)::Action
    history_index = findfirst(x->x==s.t, hr.request_timesteps)
    if history_index === nothing
        budget = EMPTY_PRODUCT_USER_BUDGET
    else
        budget = hr.h[history_index].info.b
    end
    return budget
end


text:
@testset "Common Algorithms" begin
    @testset "Constrained Output Elimination" begin
        constrained_output_elimination = Deductive.constrained_output_elimination

        # most of these test cases are stolen from the Pluto notebook linked in #21
        f1 = constrained_output_elimination(Dict(
            :a => Set([1]),
            :b => Set([1, 2])
        ))
        @test f1[:a] == Set([1])
        @test f1[:b] == Set([2])

        f2 = constrained_output_elimination(Dict(
            :a => Set([1]),
            :b => Set([2, 1]),
            :c => Set([2, 3])
        ))
        @test f2[:a] == Set([1])
        @test f2[:b] == Set([2])
        @test f2[:c] == Set([3])

        f3 = constrained_output_elimination(Dict(
            :a => Set([1]),
            :b => Set([1, 2, 3]),
            :c => Set([2, 3]),
            :d => Set([2, 3, 4]),
            :e => Set([4, 5])
        ))
        @test f3[:a] == Set([1])
        @test f3[:b] == Set([2, 3])
        @test f3[:c] == Set([2, 3])
        @test f3[:d] == Set([4])
        @test f3[:e] == Set([5])
    end
end


text:
########################################
## File Name: nominal_policies.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/14
## Description: Nominal policies for SACBP
########################################

abstract type NominalPolicy end

struct GradientMultiTargetLocalizationPolicy <: NominalPolicy end

function control_nominal(policy::GradientMultiTargetLocalizationPolicy,
                         model::UKFPosRange,
                         s::AugState{PhysPos{T},VecBelMvNormal{T}},
                         u_param_min::Vector{<:Real},
                         u_param_max::Vector{<:Real},
                         dt::Real,
                         Q::Matrix{T}) where T <: Real
    costModel = CostPosRangeLocalization();
    grad = sum(transpose.(trans_jacobi_auto_p(model,s,dt,zeros(length(s.b)),Q)).*termCost_grad_b(costModel,trans(model,s,dt,zeros(length(s.b)),Q)));
    u_val = -1.0*(grad/(norm(grad)+1e-10));
    for ii = 1:length(u_val)
        if u_val[ii] < u_param_min[ii]
            u_val[ii] = u_param_min[ii];
        elseif u_val[ii] > u_param_max[ii]
            u_val[ii] = u_param_max[ii];
        end
    end
    return PosControl(s.t,u_val)
end

struct ManipulatePositionControlPolicy <: NominalPolicy end

function control_nominal(policy::ManipulatePositionControlPolicy,
                         b::BelMvNormal{<:Real},
                         u_param_min::Vector{<:Real},
                         u_param_max::Vector{<:Real},
                         posGain::Real,
                         rotGain::Real)
    μ,Σ = b.params[1:b.dim],reshape(b.params[b.dim+1:end],b.dim,b.dim);
    fVec = -posGain.*μ[1:2];
    torque = -rotGain*(μ[3] - pi); # Target θ is pi.
    u_val = [fVec;torque];
    for ii = 1:length(u_val)
        if u_val[ii] < u_param_min[ii]
            u_val[ii] = u_param_min[ii];
        elseif u_val[ii] > u_param_max[ii]
            u_val[ii] = u_param_max[ii];
        end
    end
    return MControl2D(b.t,u_val)
end


text:
@testset "load_iris" begin
    xtmp, ytmp, vars = load_iris()

    @test typeof(xtmp) <: Matrix{Float64}
    @test typeof(ytmp) <: Vector{String}
    @test typeof(vars) <: Vector{String}
    @test size(xtmp) == (4, 150)
    @test length(ytmp) == 150
    @test length(vars) == size(xtmp, 1)
    @test mean(xtmp, dims=2) ≈ [5.843333333333333333, 3.05733333333333333, 3.758000, 1.199333333333333]
    @test mean(xtmp[:,1:50], dims=2) ≈ [5.006, 3.428, 1.462, 0.246]
    @test mean(xtmp[:,51:100], dims=2) ≈ [5.936, 2.77, 4.26, 1.326]
end


text:
using Base.Test
@testset "Broadcast and vectorization" begin
@testset "Base.broadcast" begin
    A::Vector{Float64} = 1:5
    # A = [1.0, 2.0, 3.0, 4.0, 5.0]
    B::Array{Float64, 2} = reshape(1:10,(5,2))
    # B = [1.0 6.0; 2.0 7.0; 3.0 8.0; 4.0 9.0; 5.0 10.0]
    C = broadcast(+, A, B)
    # C = [2.0 7.0; 4.0 9.0; 6.0 11.0; 8.0 13.0; 10.0 15.0]
    @test C[1,1] == 2
    @test C[5,2] == 15
end
# ============================================= #
@testset "Base.broadcast!" begin
    X::Vector{Float64} = 1:5
    # A = [1.0, 2.0, 3.0, 4.0, 5.0]
    Y::Array{Float64, 2} = reshape(1:10,(5,2))
    # B = [1.0 6.0; 2.0 7.0; 3.0 8.0; 4.0 9.0; 5.0 10.0]
    broadcast!(+, Y, X, Y)
    # B = [2.0 7.0; 4.0 9.0; 6.0 11.0; 8.0 13.0; 10.0 15.0]
    @test Y[1,1] == 2
    @test Y[5,2] == 15

    X = 1:5
    Y = reshape(1:10, (5,2))
    a::Float64 = 0.5
    broadcast!((x,y)->begin
        a*x + y
    end, Y, X, Y)

    @test Y[1,1] == 1.5
    @test Y[5,2] == 12.5
end
# ============================================= #
@testset "Dot operator" begin
    # Vectorized dot operator #
    a::Float64 = 0.5
    X::Vector{Float64} = 1:5
    Y::Array{Float64, 2} = reshape(1:10, (5,2))
    Y .= a.*X .+ Y

    @test Y[1,1] == 1.5
    @test Y[5,2] == 12.5
# --------------------------------------------- #
    # Another kind of dot expression #
    X = 1:5
    Y = reshape(1:10, (5,2))
    @. Y = a*X + Y

    @test Y[1,1] == 1.5
    @test Y[5,2] == 12.5
end

end

# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #


text:
module LineFinder

using RvSpectMLBase
using DataFrames, Query
using Statistics
using LsqFit
#using Optim
#using Plots
#using ThreadedIterables
using ..TemporalGPInterpolation

default_min_deriv2 = 3 # 1.5e5
default_smooth_factor= 4.0
default_min_pixels_in_line = 5
default_use_logλ = true
default_use_logflux = false

" Struct containing parameters for running LineFinder "
struct LineFinderPlan
  min_deriv2::Float64
  smooth_factor::Float64
  min_pixels_in_line::Int
  use_logλ::Bool
  use_logflux::Bool
end

""" Create a [LineFinderPlan](@ref)
Optional Inputs:
- min_deriv2 (default_min_deriv2)
- smooth_factor (default_smooth_factor)
- min_pixels_in_line (default_min_pixels_in_line)
- use_logλ (default_use_logλ),
- use_logflux (default_use_logflux)
"""
function LineFinderPlan(; min_deriv2::Real = default_min_deriv2, smooth_factor::Real = default_smooth_factor,
                          min_pixels_in_line::Int = default_min_pixels_in_line,
                          use_logλ::Bool = default_use_logλ, use_logflux::Bool = default_use_logflux )
  LineFinderPlan(min_deriv2, smooth_factor, min_pixels_in_line, use_logλ, use_logflux)
end

"""  `find_lines_candidates_in_chunk( chunk, plan )`
Convenience function to find links in one chunk of spectrum.
# Inputs:
- chunk: ChunkOfSpectrum to be analyzed
- deriv2: array of second derivative with respect to log λ
# Optional Arguments:
- plan: LineFinderPlan
# Returns DataFrame containing keys:
- pixels: array of ranges of pixels wihtin chunk
- pixel_max_d2fdlnλ2: pixel with maximum 2nd derivative
"""
function find_line_candidates_in_chunk(chunk::AbstractChunkOfSpectrum, deriv2::AbstractVector{T}; plan::LineFinderPlan = LineFinderPlan() ) where { T<:Real }
  idx_d2_gt_0 = findall(d2->isless(zero(d2),d2), deriv2)
  line_candidates = UnitRange[]
  first_idx_in_line = last_idx_in_line = first(idx_d2_gt_0)
    for i in idx_d2_gt_0[2:end]
      if i-last_idx_in_line > 1   # No longer in last set of contiguous indices, process last line candidate
        if last_idx_in_line - first_idx_in_line + 1 < plan.min_pixels_in_line
          # noop continue
        elseif maximum(deriv2[first_idx_in_line:last_idx_in_line]) < plan.min_deriv2 #* norm_in_chunk
          # noop continue
        else  # Add last range of pixels as line candidate
          push!(line_candidates,first_idx_in_line:last_idx_in_line)
        end
        first_idx_in_line = i
        last_idx_in_line = i
      else # ! i-last_idx_in_line>1   # Still in same line
        last_idx_in_line = i
      end
    end # for i in idx_d2_gt_0
    # See if we were in the middle of a line candidate at the end of the chunk
    if last_idx_in_line - first_idx_in_line + 1 < plan.min_pixels_in_line
      # noop continue
    elseif maximum(deriv2[first_idx_in_line:last_idx_in_line]) < plan.min_deriv2 #* norm_in_chunk
      # noop continue
    else
      push!(line_candidates,first_idx_in_line:last_idx_in_line)
    end
  pixel_to_report = map(r-> r[argmax(deriv2[r])], line_candidates)

  return DataFrame(:pixels=>line_candidates, :pixel_max_d2fdlnλ2=>pixel_to_report )
  #return line_candidates
end

" Utility function for `line_model`.  Don't assume will remain unchanged in future versions."
function pack(; a::T1, b::T2, λc::T3, depth::T4, σ²::T5) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real  }
  return [a, λc-6500, log(depth), log(σ²), b]
end

" Utility function for `line_model`.  Don't assume will remain unchanged in future versions."
function unpack(θ::AbstractVector{T} ) where {T<:Real}
  return (a=θ[1], λc=θ[2]+6500, depth=exp(θ[3]), σ²=exp(θ[4]), b=θ[5]  )
end

""" `line_model( λ, θ )`
Calculate model for one Gaussian absorption line evaluated at wavelengths λ
Input is an array that will be passed to [unpack](@ref).
Currently assumed to be chunk normalization, central wavelength, depth, σ² and slope of normalization within chunk.
"""
function line_model(λ::Union{T1,AbstractVector{T1} }, θ::AbstractVector{T2}) where { T1<:Real, T2<:Real}
  a, λc, d, σ², b = unpack(θ)
  return @. a*(1.0+b*(λ-λc))*(1.0 - d * exp(-0.5*(λ-λc)^2 /σ²) )
end

function fit_line(λ::AbstractArray{T1,1}, flux::AbstractArray{T2,1}, var::AbstractArray{T3,1}, idx::UnitRange; show_trace::Bool = false) where { T1<:Real, T2<:Real, T3<:Real }
  @assert length(λ) == length(flux) == length(var)
    λc0 = λ[idx[argmin(flux[idx])]]
  continuum0 = maximum(flux[idx])
  depth0 = 1-minimum(flux[idx])/continuum0
  width0 = 0.05^2   # might need to make more general to deal with different v sin i's
  slope0 = 0.0
  θ0 = pack(a=continuum0,b=slope0,λc=λc0,depth=depth0,σ²=width0)
  lower_bound = pack(a=0.0, b=-1.0, λc=minimum(λ[idx]), depth=0.0, σ²=0.01^2)
  upper_bound = pack(a=2.0, b= 1.0, λc=maximum(λ[idx]), depth=1.0, σ²=0.15^2)
  #upper_bound = pack(a=2.0, b= 1.0, λc=maximum(λ[idx]), depth=1.0, σ²=0.1)
  try
    res = curve_fit(line_model, λ[idx], flux[idx], 1.0 ./(var[idx]), θ0, lower=lower_bound, upper=upper_bound, show_trace=show_trace, autodiff=:forwarddiff)
    covar = estimate_covar(res)
    param = unpack(res.param)
    chi2perdof = sum(abs2.(res.resid).*res.wt) / length(idx)
    #=
    println("idx = ", idx)
    println("lambda = ",λ[idx])
    println("flux = ",flux[idx])
    println("var = ", var[idx])
    println("result = ", res)
    println("param = ", param)
    println("chi2perdof = ", chi2perdof)
    println("res.converged = ",  res.converged)
    @error("good luck!")
    =#
    return (param=param, covar=covar, χ²_per_dof=chi2perdof, is_converged = res.converged)
  catch
    param = (a=continuum0, λc=λc0, depth=depth0, σ²=width0, b=slope0  )
    covar = Array{Float64,2}(undef,0,0)
    return (param=param, covar=covar , χ²_per_dof=Inf, is_converged = false)
  end
end

""" fit_line( λ, flux, var )
Fits a basic Gaussian absorption line times a line (variable slope) to data.
Returns a tuple of (param, χ²_per_dof, is_converged)
Warning:  Has some hardcoded parameters that likely need to be generalized for different instruments.
"""
function fit_line(λ::AbstractArray{T1,1}, flux::AbstractArray{T2,1}, var::AbstractArray{T3,1} ) where { T1<:Real, T2<:Real, T3<:Real }
  @assert length(λ) == length(flux) == length(var)
  @warn "Need to check why this results in fits not converging!"
  λc0 = λ[argmin(flux)]
  continuum0 = maximum(flux)
  depth0 = 1-minimum(flux)/continuum0
  width0 = 0.05^2   # might need to make more general to deal with different v sin i's
  slope0 = 0.0
  θ0 = pack(a=continuum0,b=slope0,λc=λc0,depth=depth0,σ²=width0)
  lower_bound = pack(a=0.0, b=-1.0, λc=minimum(λ), depth=0.0, σ²=0.01^2)
  upper_bound = pack(a=2.0, b= 1.0, λc=maximum(λ), depth=1.0, σ²=0.15^2)
  try
    #res = curve_fit(line_model, λ, flux, 1.0 ./(var), θ0, lower=lower_bound, upper=upper_bound, show_trace=false, autodiff=:forwarddiff)
    λtmp = copy(λ)
    ftmp = copy(flux)
    invvar = 1.0 ./ var
    res = curve_fit(line_model, λtmp, ftmp,invvar, θ0, lower=lower_bound, upper=upper_bound, show_trace=false, autodiff=:forwarddiff)
    #res = curve_fit(line_model, λ, flux, 1.0 ./(var), θ0, lower=lower_bound, upper=upper_bound, show_trace=false, autodiff=:forwarddiff)
    param = unpack(res.param)
    covar = estimate_covar(fit)
    chi2perdof = sum(abs2.(res.resid).*res.wt) / length(λ)
    return (param=param, covar=covar, χ²_per_dof=chi2perdof, is_converged = res.converged)
  catch
    param = (a=continuum0, λc=λc0, depth=depth0, σ²=width0, b=slope0  )
    covar = Array{Float64,2}(undef,0,0)
    return (param=param, covar=covar, χ²_per_dof=Inf, is_converged = false)
  end
end

"""  find_lines_in_chunk( chunk, plan )
Convenience function to find lines in one chunk of spectrum.
# Inputs:
- chunk
# Optional Arguments:
- plan: LineFinderPlan
- chunk_id: If provide non-zero values, then add column `chunk_id` to resulting DataFrame
- obs_id: If provide non-zero values, then add column `obs_id` to resulting DataFrame
- keep_bad_fits:  If true, then does not filter out non-converged and other bad line fits from resulting dataframe
- verbose: If true, prints extra debugging info.
# Return:
- line_fit_list:  DataFrame with results of fit to each line
"""
function find_lines_in_chunk(chunk::AbstractChunkOfSpectrum; plan::LineFinderPlan = LineFinderPlan(),
                              chunk_id::Integer = 0, obs_id::Integer = 0, keep_bad_fits::Bool = false, verbose::Bool = false) where {T <:Real}

  function fit_line(idx::UnitRange)
    #RvSpectML.LineFinder.fit_line(view(chunk.λ,idx), view(flux,idx), view(chunk.var, idx) )   # some bug in this version of function causes non-convergence
    LineFinder.fit_line( chunk.λ, flux, chunk.var, idx)
  end
  norm_in_chunk = mean(chunk.flux)
  @assert 0.9<=norm_in_chunk<=1.1  # Just making sure we're passing roughly normalized chunks
  # Fit flux with extra smoothing for measuring derivatives well
  #(flux_smooth, deriv_smooth, deriv2_smooth) = predict_mean_and_derivs(chunk.λ, chunk.flux./norm_in_chunk, chunk.λ, sigmasq_obs=chunk.var./norm_in_chunk^2,
  (flux_smooth, deriv_smooth, deriv2_smooth) = predict_mean_and_derivs(chunk.λ, chunk.flux, chunk.λ, sigmasq_obs=chunk.var,
                                  smooth_factor=plan.smooth_factor, use_logx=plan.use_logλ, use_logy=plan.use_logflux )

  # Find pixel ranges that might be lines
  line_candidates = find_line_candidates_in_chunk(chunk, deriv2_smooth, plan=plan)
  if verbose    println("# Found ", size(line_candidates,1), " line candidates.")   end
  if chunk_id != 0    line_candidates[!,:chunk_id] .= chunk_id  end
  if obs_id != 0      line_candidates[!,:obs_id] .= obs_id      end
  #line_candidates[!,:pixels] = line_candidates[!,:pixels]
  line_candidates[!,:fit_min_λ] = chunk.λ[first.(line_candidates[!,:pixels])]
  line_candidates[!,:fit_max_λ] = chunk.λ[ last.(line_candidates[!,:pixels])]

  # Refit flux without extra smoothing, so depths and widths won't be biased by extra smoothing
  #(flux, deriv, deriv2) = predict_mean_and_derivs(chunk.λ, chunk.flux./norm_in_chunk, chunk.λ, sigmasq_obs=chunk.var./norm_in_chunk^2,
  (flux, deriv, deriv2) = predict_mean_and_derivs(chunk.λ, chunk.flux, chunk.λ, sigmasq_obs=chunk.var,
                                  smooth_factor=1, use_logx=plan.use_logλ, use_logy=plan.use_logflux )
  # Compute model fit parameters and add results to dataframe
  result = map(fit_line, line_candidates.pixels)
  #line_candidates[!,:fit_param] = [result[i].param for i in 1:length(result) ]
  line_candidates[!,:fit_λc] = [result[i].param.λc for i in 1:length(result) ]
  line_candidates[!,:fit_depth] = [result[i].param.depth for i in 1:length(result) ]
  line_candidates[!,:fit_σ²] = [result[i].param.σ² for i in 1:length(result) ]
  line_candidates[!,:fit_a] = [result[i].param.a for i in 1:length(result) ]
  line_candidates[!,:fit_b] = [result[i].param.b for i in 1:length(result) ]
  line_candidates[!,:fit_covar] = [result[i].covar for i in 1:length(result) ]
  line_candidates[!,:fit_χ²_per_dof] = [result[i].χ²_per_dof for i in 1:length(result) ]
  line_candidates[!,:fit_converged] = [result[i].is_converged for i in 1:length(result) ]
  line_candidates[!,:λ_near_center] = map(i->chunk.λ[line_candidates[i,:pixel_max_d2fdlnλ2]], 1:size(line_candidates,1) )
  line_candidates[!,:f_near_center] = map(i->flux_smooth[line_candidates[i,:pixel_max_d2fdlnλ2]], 1:size(line_candidates,1) )
  line_candidates[!,:dfdlnλ_near_center] = map(i->deriv_smooth[line_candidates[i,:pixel_max_d2fdlnλ2]], 1:size(line_candidates,1) )
  line_candidates[!,:df2dlnλ2_near_center] = map(i->deriv2_smooth[line_candidates[i,:pixel_max_d2fdlnλ2]], 1:size(line_candidates,1) )
  line_candidates[!,:λ_min_flux] = line_candidates[!,:λ_near_center].-line_candidates[!,:dfdlnλ_near_center]./(line_candidates[!,:df2dlnλ2_near_center])
  line_candidates[!,:dfdlnλ_min_flux] = line_candidates[!,:dfdlnλ_near_center].+log.(line_candidates[!,:λ_min_flux]./line_candidates[!,:λ_near_center]).*line_candidates[!,:df2dlnλ2_near_center]

  if ! keep_bad_fits
    # Avoid "lines" at edge of chunk, either due to contamination or GP fit using zero mean
    line_candidates = line_candidates |> @filter(minimum(_.pixels)>1) |> @filter(maximum(_.pixels)<length(flux)) |> #DataFrame
    #line_candidates = line_candidates |>
        @filter(_.fit_converged) |>
        #@filter(_.fit_σ²<0.018) |>
        @filter(-0.5 <_.fit_b <0.5) |> @filter(_.fit_a < 1.8) |>
        DataFrame
  end

  return line_candidates
end

"""  find_lines_in_chunklist ( chunklist, line_list )
Convenience function to find lines in each chunk of a spectrum.
# Inputs:
- chunklist
# Optional Arguments:
- plan:  LineFinderPlan
- keep_bad_fits:  If true, then does not filter out non-converged and other bad line fits from resulting dataframe
- verbose: If true, prints extra debugging info.
# Return:
- line_list: DataFrame results of fit to each line
"""
function find_lines_in_chunklist(chunk_list::AbstractChunkList ; obs_id::Integer = 0, plan::LineFinderPlan = LineFinderPlan(), keep_bad_fits::Bool = false, verbose::Bool = false )
  mapreduce(i->find_lines_in_chunk(chunk_list[i], plan=plan, chunk_id=i, obs_id=obs_id, keep_bad_fits=keep_bad_fits, verbose=verbose ), append!, 1:length(chunk_list) )
end

"""  `find_lines_in_chunklist_timeseries( chunklist_timeseries, line_list )`
Convenience function to find lines in each chunk of each spectrum in a timeseries.
# Inputs:
- chunklist_timeseries
# Optional Arguments:
- plan:  LineFinderPlan
- keep_bad_fits:  If true, then does not filter out non-converged and other bad line fits from resulting dataframe
- verbose: If true, prints extra debugging info.
# Return:
- line_list: Array of results from [find_lines_in_chunklist](@ref)
"""
function find_lines_in_chunklist_timeseries(clt::AbstractChunkListTimeseries ; plan::LineFinderPlan = LineFinderPlan(), keep_bad_fits::Bool = false, verbose::Bool = false  )
    #@threaded
    map(cl->find_lines_in_chunklist(clt.chunk_list[cl], obsid=cl, plan=plan, keep_bad_fits=keep_bad_fits, verbose=verbose ),1:num_obs(clt))
end

""" `find_pixels_for_all_lines_in_chunklist( chunk_list, lines)`
Return array giving range of pixels for each line within provided chunklist
Inputs:
- chunk_list
- lines: DataFrame containing keys: `:fit_min_λ`, `:fit_max_λ`, `:chunk_id`
Output:
- Vector of pixel ranges
"""
function find_pixels_for_all_lines_in_chunklist( chunk_list::AbstractChunkList, lines::DataFrame)
  @assert hasproperty(lines,:fit_min_λ)
  @assert hasproperty(lines,:fit_max_λ)
  @assert hasproperty(lines,:chunk_id)
  map(l->find_pixels_for_line_in_chunklist(chunk_list[1], lines[l,:fit_min_λ],  lines[l,:fit_max_λ], lines[l,:chunk_id] ), 1:size(lines,1) )
end

global fit_line_in_chunklist_timeseries_count_msgs = 0

""" `fit_line_in_chunklist_timeseries( chunk_list_timeseries, λmin, λmax, chunk_index)`
Return DataFrame with results of fits to each line in a given chunk of chunk_list_timeseries (including each observation time)
Inputs:
- chunk_list_timeseries: Data to fit
- λmin
- λmax
- chunk_index:  Restricts fitting to specified chunk
Outputs a DataFrame with keys:
- fit_a: Normalization of continuum around line
- fit_λc: Central wavelength
- fit_depth: Line depth
- fit_σ²:
- fit_b: Slope of continuum around line
- fit_covar: Covariance matrix for fit parameters
- χ²_per_dof: quality of fit to chunk
- fit_converged: Bool indicating if `LsqFit` converged
- obs_idx: index of observation in chunk_list_timeseries
- chunk_id: index of chunk in chunk_list_timeseries
- pixels: range of pixels that was fit
"""
function fit_line_in_chunklist_timeseries(clt::AbstractChunkListTimeseries, λmin::Real, λmax::Real, chid::Integer; show_trace::Bool = false)
  df = DataFrame()
  nobs = length(clt.chunk_list)
  ParamT = NamedTuple{(:a, :λc, :depth, :σ², :b),NTuple{5,Float64}}  # maybe generalize by making part of plan?
  param = Vector{ParamT}(undef,nobs)
  fit_a = Vector{Float64}(undef,nobs)
  fit_λc = Vector{Float64}(undef,nobs)
  fit_depth = Vector{Float64}(undef,nobs)
  fit_σ² = Vector{Float64}(undef,nobs)
  fit_b = Vector{Float64}(undef,nobs)
  fit_covar = Vector{Array{Float64,2}}(undef,nobs)
  fit_converged = falses(nobs)
  pixels = Vector{UnitRange{Int64}}(undef,nobs)
  gof = zeros(nobs)
  global count_msgs
  for t in 1:nobs
    pixels[t] = RvSpectMLBase.find_pixels_for_line_in_chunklist(clt.chunk_list[t], λmin, λmax, chid).pixels
    mean_flux = mean(clt.chunk_list[t][chid].flux[pixels[t]])
    flux = clt.chunk_list[t][chid].flux ./ mean_flux
    var = clt.chunk_list[t][chid].var ./ mean_flux^2
    (param_tmp, fit_covar[t], gof[t], fit_converged[t] ) = fit_line(clt.chunk_list[t][chid].λ, flux, var, pixels[t], show_trace=show_trace )
    #=
    if fit_line_in_chunklist_timeseries_count_msgs<5
      (param_tmp, fit_covar[t], gof[t], fit_converged[t] ) = fit_line(clt.chunk_list[t][chid].λ, flux, var, pixels[t], show_trace=true )
        println("lambda = ",extrema(clt.chunk_list[t][chid].λ[pixels[t]]))
      println("flux = ",extrema(clt.chunk_list[t][chid].flux[pixels[t]]))
      println("var = ", extrema(clt.chunk_list[t][chid].var[pixels[t]]))
      println("flux = ",extrema(flux))
      println("var = ", extrema(var))
      println("pixels = ",extrema(pixels[t]))
      println("param = ", param_tmp)
      println("gof = ", gof[t])
      println("fit_converged = ", fit_converged[t])
      fit_line_in_chunklist_timeseries_count_msgs += 1
    else
      (param_tmp, fit_covar[t], gof[t], fit_converged[t] ) = fit_line(clt.chunk_list[t][chid].λ, flux, var, pixels[t] )
    end
    =#
    fit_a[t] = param_tmp.a
    fit_λc[t] = param_tmp.λc
    fit_depth[t] = param_tmp.depth
    fit_σ²[t] = param_tmp.σ²
    fit_b[t] = param_tmp.b
  end

  return DataFrame(:fit_a=>fit_a, :fit_λc=>fit_λc, :fit_depth=>fit_depth, :fit_σ²=>fit_σ², :fit_b=>fit_b,
                    :fit_covar=>fit_covar, :χ²_per_dof=>gof, :fit_converged=>fit_converged, :obs_idx =>collect(1:nobs), :chunk_id=>chid, :pixels=>pixels  )
end

""" `fit_line_in_chunklist_timeseries( chunk_list_timeseries, λmin, λmax )`
Return DataFrame with results of fits to each line in a chunk_list_timeseries (including each observation time)
Inputs:
- chunk_list_timeseries: Data to fit
- λmin
- λmax
- chunk_index:  Restricts fitting to specified chunk
Outputs a DataFrame with keys:
- fit_a: Normalization of continuum around line
- fit_λc: Central wavelength
- fit_depth: Line depth
- fit_σ²:
- fit_b: Slope of continuum around line
- fit_covar: Covariance matrix for fit parameters
- χ²_per_dof: quality of fit to chunk
- fit_converged: Bool indicating if `LsqFit` converged
- obs_idx: index of observation in chunk_list_timeseries
- chunk_id: index of chunk in chunk_list_timeseries
- pixels: range of pixels that was fit
"""
function fit_all_lines_in_chunklist_timeseries(clt::AbstractChunkListTimeseries, lines::DataFrame ; plan::LineFinderPlan = LineFinderPlan(), show_trace::Bool = false )
  @assert size(lines,1) >= 2
  line_idx::Int64 = 1
  λmin = lines[line_idx,:fit_min_λ]
  λmax = lines[line_idx,:fit_max_λ]
  chid = lines[line_idx,:chunk_id]
  df = fit_line_in_chunklist_timeseries(clt, λmin, λmax, chid, show_trace=show_trace)
  df[!,:line_id] .= line_idx
  for line_idx in 2:size(lines,1)
    λmin = lines[line_idx,:fit_min_λ]
    λmax = lines[line_idx,:fit_max_λ]
    chid = lines[line_idx,:chunk_id]
    df_tmp = fit_line_in_chunklist_timeseries(clt, λmin, λmax, chid)
    df_tmp[!,:line_id] .= line_idx
    append!(df,df_tmp)
  end
  return df
end

""" `find_worst_telluric_in_each_line_fit!( df, clt, data)`
Patterned after find_worst_telluric_in_each_chunk in src/instruments/tellurics.jl
TODO: Move into EXPRES namespace?
"""
function find_worst_telluric_in_each_line_fit!( df::DataFrame, clt::AbstractChunkListTimeseries, data::AbstractArray{AS,1} )  where {AS<:AbstractSpectra}
  @assert hasproperty(df,:line_id)
  @warn "find_worst_telluric_in_each_line_fit! will likely be moved soon"
  min_telluric_model_this_obs = ones(size(df,1))
  min_telluric_model_all_obs  = ones(size(df,1))
  for (i, row) in enumerate(eachrow(df))
    t_idx = row.obs_idx
    ch_idx = row.chunk_id
    pixels = row.pixels
    view_indices = clt.chunk_list[t_idx].data[ch_idx].λ.indices
    cols = view_indices[1][pixels]
    order = view_indices[2]
    min_telluric_model_this_obs[i] = minimum(view(data[t_idx].metadata[:tellurics], cols, order))
  end
  df.min_telluric_model_this_obs = min_telluric_model_this_obs
  return df
  df_tmp = df |> @groupby(_.line_id) |> @map( { min_telluric_model_all_obs=minimum(_.min_telluric_model_this_obs), line_id=first(_.line_id) } ) |> DataFrame
  df.min_telluric_model_all_obs = df_tmp[telluric_info[!,:line_id],:min_telluric_model_all_obs]
  return df
end

end # module LineFinder


text:
"Default PDDL requirements."
const DEFAULT_REQUIREMENTS = Dict{Symbol,Bool}(
    Symbol("strips") => true,
    Symbol("typing") => false,
    Symbol("negative-preconditions") => false,
    Symbol("disjunctive-preconditions") => false,
    Symbol("equality") => false,
    Symbol("existential-preconditions") => false,
    Symbol("universal-preconditions") => false,
    Symbol("quantified-preconditions") => false,
    Symbol("conditional-effects") => false,
    Symbol("adl") => false,
    Symbol("derived-predicates") => false
)

const IMPLIED_REQUIREMENTS = Dict{Symbol,Vector{Symbol}}(
    Symbol("adl") => [
        Symbol("strips"),
        Symbol("typing"),
        Symbol("negative-preconditions"),
        Symbol("disjunctive-preconditions"),
        Symbol("equality"),
        Symbol("quantified-preconditions"),
        Symbol("conditional-effects")
    ],
    Symbol("quantified-preconditions") => [
        Symbol("existential-preconditions"),
        Symbol("universal-preconditions")
    ]
)