text:
# binary.jl: Elementwise broadcasting binary functions for arrays and scalars.
# uses binary_ops from broadcast.jl.

import Base.Broadcast: broadcasted

# binary_op defines the broadcast_func of a Julia function for KnetArrays.
# The corresponding kernel is defined in libknet8.
function binary_op(f, j=f, o...)
    J=Symbol(j)
    M = which(@__MODULE__, J)
    for S in (32,64)
        T = Symbol("Float$S")

        F01 = "$(f)_$(S)_01" # Scalar,Array->Array
        F11 = "$(f)_$(S)_11" # Array,Array->Array (same size) (not broadcast)
        F12 = "$(f)_$(S)_12" # Array,Array->Array (one have to be vector)
        F13_x_y = "$(f)_$(S)_13_x_y" # e.g. (A(x,y,z,w,t...), B(1,1,1,w,1...))
        F13_y_x = "$(f)_$(S)_13_y_x" # different versions for efficiency
        # F14_x_y = "$(f)_$(S)_14_x_y" # e.g. (M(x,y,z,w,t...), N(w,1,1,1...)
        # F14_y_x = "$(f)_$(S)_14_y_x" # different versions for efficiency
        # F15 reserved for another kernel, eliminated later and combined with F16
        F16_3 = "$(f)_$(S)_16_3" # multi-dimensional bcast unrolled up to 5 dims
        F16_4 = "$(f)_$(S)_16_4" # multi-dimensional bcast unrolled up to 5 dims
        F16_5 = "$(f)_$(S)_16_5" # multi-dimensional bcast unrolled up to 5 dims
        F17 = "$(f)_$(S)_17" # multi-dimensional bcast with loops

        @eval begin
            # Scalar,Array->Array
            function broadcasted(::typeof($J),x::$T,y::KnetArray{$T})
                z = similar(y)
                @knet8($F01,(Cint,$T,Ptr{$T},Ptr{$T}),length(z),x,y,z)
                return z
            end
            # Array,Array->Array
            function broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T})
                if size(x)==size(y)
                    z = similar(x)
                    @knet8($F11,(Cint,Ptr{$T},Ptr{$T},Ptr{$T}),length(z),x,y,z)
                    return z
                else
                    bs = vbroadcast_shape(x,y)
                    z = similar(x,bs[1])
                    _broadcasted($J,x,y,z,bs)
                end
            end
            # Helpers
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T,1},bs)
                if length(x) == 1
                    broadcasted($J,x[1],y)
                elseif length(y) == 1
                    broadcasted($J,x,y[1])
                else # length(x) == length(y) was handled above
                    throw(DimensionMismatch("$(map(size,(x,y,z)))"))
                end
            end
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T,2},bs)
                # xlast or ylast will be broadcasting dimension
                (dz,sx,nx,sy,ny,xlast,ylast,xdims,ydims,multi) = bs
                if (nx == 1
                    || ny == 1
                    || ((xdims == 1 && (xlast==1 || 512 < sx )) ||
                        (ydims == 1 && (ylast==1 || sy < 512 )))
                    || (xdims==1 && nx<704)
                    || (ydims==1 && (ny<704)))
                    @knet8($F12,
                           (Cint,Ptr{$T},Cint,Cint,Ptr{$T},Cint,Cint,Ptr{$T}),
                           length(z),x,sx,nx,y,sy,ny,z)
                elseif xdims == 1
                    dim_stride = strides(y)[xlast]
                    next_stride = (xlast+1) > ndims(y) ?
                        0 : strides(y)[xlast+1]
                    dim_size = prod(size(y)[xlast+1:end])
                    @knet8($F13_y_x,
                           (Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint),
                           y,x,z,dim_stride,next_stride,dim_size,
                           length(y),length(x))
                elseif ydims == 1
                    dim_stride = strides(x)[ylast]
                    next_stride = (ylast+1) > ndims(x) ?
                        0 : strides(x)[ylast+1]
                    dim_size = prod(size(x)[ylast+1:end])
                    @knet8($F13_x_y,
                           (Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint),
                           x,y,z,dim_stride,next_stride,dim_size,
                           length(x), length(y))
                else
                    throw(DimensionMismatch("$(map(size,(x,y,z)))"))
                end
                return z
            end
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T,3},bs)
                sx,sy,sz = get_strides(x,y,z)
                @knet8($F16_3,
                       (Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                       x,y,z, sx[1],sx[2],sx[3], sy[1],sy[2],sy[3], sz[1],sz[2],sz[3], length(z))
                return z
            end
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T,4},bs)
                sx,sy,sz = get_strides(x,y,z)
                @knet8($F16_4,
                       (Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                       x,y,z, sx[1],sx[2],sx[3],sx[4], sy[1],sy[2],sy[3],sy[4], sz[1],sz[2],sz[3],sz[4], length(z))
                return z
            end
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T,5},bs)
                sx,sy,sz = get_strides(x,y,z)
                @knet8($F16_5,
                       (Ptr{$T},Ptr{$T},Ptr{$T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
                       x,y,z, sx[1],sx[2],sx[3],sx[4],sx[5], sy[1],sy[2],sy[3],sy[4],sy[5], sz[1],sz[2],sz[3],sz[4],sz[5], length(z))
                return z
            end
            function _broadcasted(::typeof($J),x::KnetArray{$T},y::KnetArray{$T},z::KnetArray{$T},bs)
                # ndims(z) <= 5 handled above, this is for > 5
                sx,sy,sz = map(s->convert(KnetArray, s), get_strides(x,y,z))
                @knet8($F17,
                       (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint},Ptr{Cint},Ptr{Cint},Cint,Cint),
                       x,y,z, sx, sy, sz, length(z), ndims(z))
                return z
            end

            # Bcasted methods
            ($M).$J(x::Bcasted{<:KnetArray{$T}}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            ($M).$J(x::KnetArray{$T}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            ($M).$J(x::Bcasted{<:KnetArray{$T}}, y::KnetArray{$T}) = broadcasted($J, x.value, y) |> Bcasted
            ($M).$J(x::Bcasted{$T}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            ($M).$J(x::$T, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            ($M).$J(x::Bcasted{$T}, y::KnetArray{$T}) = broadcasted($J, x.value, y) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:KnetArray{$T}}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            broadcasted(::typeof($J),x::KnetArray{$T}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:KnetArray{$T}}, y::KnetArray{$T}) = broadcasted($J, x.value, y) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{$T}, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            broadcasted(::typeof($J),x::$T, y::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{$T}, y::KnetArray{$T}) = broadcasted($J, x.value, y) |> Bcasted
        end # @eval
    end # for
    @eval begin # so we do not trigger some default Base implementation
        ($M).$J(x::Bcasted, y::Bcasted) = throw(MethodError($J,(x,y)))
        ($M).$J(x, y::Bcasted) = throw(MethodError($J,(x,y)))
        ($M).$J(x::Bcasted, y) = throw(MethodError($J,(x,y)))
        broadcasted(::typeof($J),x::Bcasted, y::Bcasted) = throw(MethodError(broadcasted, ($J, x, y)))
        broadcasted(::typeof($J),x, y::Bcasted) = throw(MethodError(broadcasted, ($J, x, y)))
        broadcasted(::typeof($J),x::Bcasted, y) = throw(MethodError(broadcasted, ($J, x, y)))
    end
end # function binary_op

# vbroadcast_shape computes index/offset arguments for a broadcasting kernel call.
function vbroadcast_shape(x,y)
    nz = max(ndims(x),ndims(y))
    dz = ones(Int,nz)
    xdims = ydims = xsame = ysame = xlast = ylast = 0; zlen = 1;
    sx = sy = -1;
    nx = ny = -1;

    for i = 1:nz
        # xdims: number of xdims != 1
        # xlast: last index != 1

        if size(x,i) > 1
            xdims += 1; xlast = i
            dz[i] = size(x,i)
        end
        if size(y,i) > 1
            ydims += 1; ylast = i
            if dz[i] == 1
                dz[i] = size(y,i)
            elseif dz[i] != size(y,i)
                throw(DimensionMismatch(
                    "arrays could not be broadcast to a common size"))
            end
        end

        xsame += (dz[i] == size(x,i))
        ysame += (dz[i] == size(y,i))
        zlen *= dz[i]
    end

    multi = false
    if (xsame != nz && xdims > 1) || (ysame != nz && ydims > 1)
        multi = true
    end

    if xdims == 0
        sx = zlen; nx = 1
    elseif xdims == 1
        sx = prod(dz[1:xlast-1]); nx=dz[xlast]
    elseif xsame == nz
        sx = 1
        nx = zlen
    elseif !multi
        error("Broadcasting error")
    end

    if ydims == 0
        sy = zlen
        ny = 1
    elseif ydims == 1
        sy = prod(dz[1:ylast-1])
        ny = dz[ylast]
    elseif ysame == nz
        sy = 1
        ny = zlen
    elseif !multi
        error("Broadcasting error")
    end

    return (tuple(dz...), sx, nx, sy, ny,xlast,ylast,xdims,ydims,multi)
end

function get_strides(x,y,z)
    # x,y,z may have different ndims, work with the max
    n = max(ndims(x), ndims(y), ndims(z))
    stride_x = Int32[ stride(x,i) for i=1:n ]
    stride_y = Int32[ stride(y,i) for i=1:n ]
    stride_z = Int32[ stride(z,i) for i=1:n ]
    dims_x = Int32[ size(x,i) for i=1:n ]
    dims_y = Int32[ size(y,i) for i=1:n ]
    for i in 1:n
        dims_x[i] == dims_y[i] && continue
        if dims_x[i]==1
            stride_x[i]=0
        else
            stride_y[i]=0
        end
    end
    return stride_x, stride_y, stride_z
end

# Additional imports: fns in binary_ops are defined using broadcasted.
import Base: +, -, *, /, \

# Here we'll just define some functions that specifically do not have broadcasting.
(+)(x::KnetArray{T},y::KnetArray{T}) where {T} = (size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.+)(x,y))
(-)(x::KnetArray{T},y::KnetArray{T}) where {T} = (size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.-)(x,y))
#(*){T}(x::KnetArray{T},y::KnetArray{T})=(.*)(x,y) # This is matmul
#(/){T}(x::KnetArray{T},y::KnetArray{T})=(./)(x,y) # This is another linalg op

# Broadcast max/min haven't been defined in Base:
# max(a::Array,b::Array)=broadcast(max,a,b)
# min(a::Array,b::Array)=broadcast(min,a,b)
# tkelman: These two methods aren't necessary, and overwrite Base. You can get this behavior via max.(a,b), with @compat needed on 0.4.

# import Base: broadcast

# Scalar kernels are defined for scalar,array order only.
# For array,scalar we can get most for free.
@eval begin
    broadcasted(::typeof(+),a::KnetArray{T},s::Number) where {T} = (.+)(T(s),a)
    broadcasted(::typeof(+),s::Number,a::KnetArray{T}) where {T} = (.+)(T(s),a)
    broadcasted(::typeof(-),a::KnetArray{T},s::Number) where {T} = (.+)(T(-s),a)
    broadcasted(::typeof(-),s::Number,a::KnetArray{T}) where {T} = (.-)(T(s),a)
    broadcasted(::typeof(*),a::KnetArray{T},s::Number) where {T} = (.*)(T(s),a)
    broadcasted(::typeof(*),s::Number,a::KnetArray{T}) where {T} = (.*)(T(s),a)
    broadcasted(::typeof(/),a::KnetArray{T},s::Number) where {T} = (.*)(T(1/s),a)
    broadcasted(::typeof(/),s::Number,a::KnetArray{T}) where {T} = (./)(T(s),a)
    broadcasted(::typeof(max),a::KnetArray{T},s::Number) where {T} = max.(T(s),a)
    broadcasted(::typeof(max),s::Number,a::KnetArray{T}) where {T} = max.(T(s),a)
    broadcasted(::typeof(min),a::KnetArray{T},s::Number) where {T} = min.(T(s),a)
    broadcasted(::typeof(min),s::Number,a::KnetArray{T}) where {T} = min.(T(s),a)
    broadcasted(::typeof(^),s::Number,a::KnetArray{T}) where {T} = (.^)(T(s),a)
    # Pow is the one exception, we need to define a separate kernel:
    rpow(s,a)=a^s # only broadcast#rpow is defined above, we need rpow defined
    broadcasted(::typeof(^),a::KnetArray{T},s::Number) where {T} = rpow.(T(s),a)
end

@eval begin
    broadcasted(::typeof(==),a::KnetArray{T},s::Number) where {T} = (T(s).==a)
    broadcasted(::typeof(==),s::Number,a::KnetArray{T}) where {T} = (T(s).==a)
    broadcasted(::typeof(!=),a::KnetArray{T},s::Number) where {T} = (T(s).!=a)
    broadcasted(::typeof(!=),s::Number,a::KnetArray{T}) where {T} = (T(s).!=a)
    broadcasted(::typeof(>),a::KnetArray{T},s::Number) where {T} = (T(s).<a)
    broadcasted(::typeof(>),s::Number,a::KnetArray{T}) where {T} = (T(s).>a)
    broadcasted(::typeof(>=),a::KnetArray{T},s::Number) where {T} = (T(s).<=a)
    broadcasted(::typeof(>=),s::Number,a::KnetArray{T}) where {T} = (T(s).>=a)
    broadcasted(::typeof(<),a::KnetArray{T},s::Number) where {T} = (T(s).>a)
    broadcasted(::typeof(<),s::Number,a::KnetArray{T}) where {T} = (T(s).<a)
    broadcasted(::typeof(<=),a::KnetArray{T},s::Number) where {T} = (T(s).>=a)
    broadcasted(::typeof(<=),s::Number,a::KnetArray{T}) where {T} = (T(s).<=a)
end

# Bcasted methods

for f in Symbol.((+,-,*,/,max,min,^,==,!=,>,>=,<,<=))
    M = which(@__MODULE__, f)
    @eval begin
        broadcasted(::typeof($f),s::Bcasted{<:Number},a::Bcasted{<:KnetArray}) = broadcasted($f, s.value, a.value) |> Bcasted
        broadcasted(::typeof($f),s::Bcasted{<:Number},a::KnetArray) = broadcasted($f, s.value, a) |> Bcasted
        broadcasted(::typeof($f),s::Number,a::Bcasted{<:KnetArray}) = broadcasted($f, s, a.value) |> Bcasted
        broadcasted(::typeof($f),a::Bcasted{<:KnetArray},s::Bcasted{<:Number}) = broadcasted($f, a.value, s.value) |> Bcasted
        broadcasted(::typeof($f),a::KnetArray,s::Bcasted{<:Number}) = broadcasted($f, a, s.value) |> Bcasted
        broadcasted(::typeof($f),a::Bcasted{<:KnetArray},s::Number) = broadcasted($f, a.value, s) |> Bcasted
        ($M).$f(s::Bcasted{<:Number},a::Bcasted{<:KnetArray}) = broadcasted($f, s.value, a.value) |> Bcasted
        ($M).$f(s::Bcasted{<:Number},a::KnetArray) = broadcasted($f, s.value, a) |> Bcasted
        ($M).$f(s::Number,a::Bcasted{<:KnetArray}) = broadcasted($f, s, a.value) |> Bcasted
        ($M).$f(a::Bcasted{<:KnetArray},s::Bcasted{<:Number}) = broadcasted($f, a.value, s.value) |> Bcasted
        ($M).$f(a::KnetArray,s::Bcasted{<:Number}) = broadcasted($f, a, s.value) |> Bcasted
        ($M).$f(a::Bcasted{<:KnetArray},s::Number) = broadcasted($f, a.value, s) |> Bcasted
    end
end

# familiar aliases for broadcasting operations of array & scalar (#7226):
# (+)(a::KnetArray{T},s::Number) where {T} = (.+)(T(s),a)  -- deprecated
# (+)(s::Number,a::KnetArray{T}) where {T} = (.+)(T(s),a)  -- deprecated
# (-)(a::KnetArray{T},s::Number) where {T} = (.+)(T(-s),a) -- deprecated
# (-)(s::Number,a::KnetArray{T}) where {T} = (.-)(T(s),a)  -- deprecated
(*)(a::KnetArray{T},s::Number) where {T} = (.*)(T(s),a)
(*)(s::Number,a::KnetArray{T}) where {T} = (.*)(T(s),a)
(/)(a::KnetArray{T},s::Number) where {T} = (.*)(T(1/s),a)
(\)(s::Number,a::KnetArray{T}) where {T} = (.*)(T(1/s),a)

#(/)(s::Number,a::KnetArray{T}) where {T} = (.*)(T(1/s),a) # TODO: non-elementwise definition in linalg
#(^)(a::KnetArray{T},s::Number) where {T} = (.^)(a,T(s)) # non-elementwise definition in linalg
#(^)(s::Number,a::KnetArray{T}) where {T} = (.^)(T(s),a) # non-elementwise definition in linalg

# Define all overloaded Julia functions for KnetArrays:

for f in binary_ops
    if !isa(f,Tuple); f=(f,); end
    binary_op(f...)
end

# Fix #412 where KnetArray(randn(Float64,4,4,4,4)).^2 gives a 1-D result
broadcasted(::typeof(Base.literal_pow), ::typeof(^), k::KnetArray{T}, n::Val{N}) where {T,N} = broadcasted(^, k, N)


text:
module DynamicMovementPrimitives
using ODE, Requires

export DMPopts, centraldiff,fit, solve, force, acceleration, solve_canonical, plotdmp, plotdmpphase

function centraldiff(v::AbstractMatrix)
    dv = diff(v)/2
    a1 = [dv[1,:];dv]
    a2 = [dv;dv[end,:]]
    a = a1+a2
end

function centraldiff(v::AbstractVector)
    dv = diff(v)/2
    a1 = [dv[1];dv]
    a2 = [dv;dv[end]]
    a = a1+a2
end

"""Takes an n vector of m vectors and creates a n×m matrix"""
vv2m(x::Vector) = [x[i][j] for i in eachindex(x), j in eachindex(x[1])]

"""
`DMPopts(Nbasis,αx,αz) = DMPopts(Nbasis, αx, αz, βz = αz/4)`\n
Holds parameters for fitting a DMP
# Fields
`Nbasis,αx,αz,βz,sched_sig`\n
`sched_sig` can be chosen as `:canonical` (default), `:time` or `position`

See example file or the paper by Ijspeert et al. 2013
"""
immutable DMPopts
    Nbasis::Int
    αx::Real
    αz::Real
    βz::Real
    sched_sig::Symbol
end

DMPopts(Nbasis,αx,αz) = DMPopts(Nbasis,αx,αz,αz/4,:position)
DMPopts(Nbasis,αx,αz,sched_sig::Symbol) = DMPopts(Nbasis,αx,αz,αz/4,sched_sig)

"""
The result of fitting a DMP
#Fields
`opts,g,y,ẏ,ÿ,w,τ,c,σ2`\n
See example file or the paper by Ijspeert et al. 2013
"""
immutable DMP
    opts::DMPopts
    g::Vector{Float64}
    y::Matrix{Float64}
    ẏ::Matrix{Float64}
    ÿ::Matrix{Float64}
    t::AbstractVector
    w::Matrix{Float64}
    τ::Float64
    c::VecOrMat{Float64}
    σ2::VecOrMat{Float64}
end

function get_centers_linear(Nbasis,x)
    ma = maximum(x,1)
    mi = minimum(x,1)
    n = size(x,2)
    d = ma-mi
    Ni = d./(Nbasis+1)

    σ2 = zeros(Nbasis,n)
    c = zeros(Nbasis,n)
    for i = 1:n
        σ2[:,i] = d[i]*(0.5/Nbasis)^2 * ones(Nbasis)
        c[:,i]  = linspace(mi[i]+0Ni[i],ma[i]-0Ni[i],Nbasis)
    end
    return c, σ2
end

# function get_centers_log(Nbasis)
#     Ni = 1/(Nbasis+1)
#     return (logspace(log10(Ni),log10(1-Ni),Nbasis))[end:-1:1]
# end

comp(x) = (x)

function kernel_matrix(x::AbstractVecOrMat,c,σ2)
    Ψ = Float64[exp(-1/(2σ2[j])*(comp(x)-(c[j]))^2) for x in x, j in eachindex(c)]
    Ψ ./= sum(Ψ,2)
    Ψ
end

function kernel_vector(x::Real,c,σ2)
    Ψ = Float64[exp(-1/(2σ2[j])*(comp(x)-(c[j]))^2) for j in eachindex(c)]
    Ψ ./= sum(Ψ)
    Ψ
end

solve_canonical(αx,τ,T::AbstractVector) = exp(-αx/τ.*T)
solve_canonical(αx,τ,T::Real) = solve_canonical(αx,τ,(0:T-1))
solve_canonical(dmp::DMP,t) = solve_canonical(dmp.opts.αx, dmp.τ, t)
_1(y::VecOrMat) = y[1,:][:]
_1(dmp::DMP) = _1(dmp.y)
_T(dmp::DMP) = size(dmp.y,1)

function get_sched_sig(s,αx,τ,t,y,g)
    if s == :canonical
        return solve_canonical(αx,τ,t)
    elseif s == :position
        return g'.-y
    elseif s == :time
        return t
    end
    warn("Scheduling signal $s unknown")
    return solve_canonical(αx,τ,t)
end
get_sched_sig(dmp::DMP) = get_sched_sig(dmp.opts.sched_sig,dmp.opts.αx,dmp.τ,dmp.t,dmp.g)

"""
`fit(y, ẏ, ÿ, opts, g=y[end])`

Fits a DMP to data\n
`y, ẏ, ÿ` are position, velocity and acceleration respectively, `T×n` matrices where `T` is the number of time steps and `n` is the number of degrees of freedom.

see also `solve`, `plotdmp`
"""
function fit(y,ẏ,ÿ,t,opts,g=y[end,:][:])
    T       = t[end]
    N       = size(y,1)
    n       = isa(y,Matrix) ? size(y,2) : 1
    τ       = T/3 # After three time constants we have reached 1-exp(-3) ≈ 0.95
    Nbasis  = opts.Nbasis
    βz      = opts.βz
    αz      = opts.αz
    αx      = opts.αx
    y0      = _1(y)
    x       = get_sched_sig(opts.sched_sig,αx,τ,t,y,g)
    ft      = τ^2*ÿ - αz*(βz*(g'.-y)-τ*ẏ)
    ξ       = x.*(g-y0)'
    c,σ2    = get_centers_linear(Nbasis,x)
    if opts.sched_sig != :position
        Ψ = kernel_matrix(x,c,σ2)
    end
    w = zeros(Nbasis,n)
    for i = 1:n #joints
        if opts.sched_sig == :position
            Ψ = kernel_matrix(x[:,i],c[:,i],σ2[:,i])
        end
        for j = 1:Nbasis
            sTΓ = ξ[:,i].*Ψ[:,j]
            w[j,i] = vecdot(sTΓ,ft[:,i])/vecdot(sTΓ,ξ[:,i])
        end
    end
    return DMP(opts, g, y, ẏ, ÿ,t, w, τ,c,σ2)
end

"""
`force(dmp,x)`

Calculate the forcing term for `dmp` when the phase variable is `x`\n
The return value will be an `n` Vector or `T×n` Matrix depending on `typeof(x)`
"""
function force(d::DMP, x)
    if d.opts.sched_sig == :position
        force_multiple(d,x)
    else
        force_single(d,x)
    end
end

function force(d::DMP, x,i)
    if d.opts.sched_sig == :position
        force_multiple(d,x,i)
    else
        force_single(d,x,i)
    end
end

function force_single(d::DMP,x::Number, i)
    # ODE Point case
    y0 = _1(d)
    Ψ  = kernel_vector(x, d.c, d.σ2)
    f = vecdot(Ψ,d.w[:,i]) * x*(d.g[i]-y0[i])
end

function force_single(d::DMP,x::Number)
    # Point case
    y0 = _1(d)
    Ψ  = kernel_matrix([x], d.c, d.σ2)
    f  = (Ψ*d.w)[:] .* (x*(d.g-y0))
end

function force_single(d::DMP,x::AbstractVector)
    # Trajectory case
    y0 = _1(d)
    Ψ    = kernel_matrix(x, d.c, d.σ2)
    f = Ψ*d.w .* (x.*(d.g-y0)')
end

function force_multiple(d::DMP,x::Number, i)
    # ODE Point case
    y0 = _1(d)
    Ψ  = kernel_vector(x, d.c[:,i], d.σ2[:,i])
    f  = vecdot(Ψ,d.w[:,i]) * x*(d.g[i]-y0[i])
end

function force_multiple(d::DMP,x::AbstractVector)
    # Point case
    n = size(d.c,2) # Number of DOF
    y0 = _1(d)
    f = zeros(n)
    for i = 1:n
        Ψ    = kernel_matrix([x[i]], d.c[:,i], d.σ2[:,i])
        f[i] = vecdot(Ψ,d.w[:,i]) * x[i].*(d.g[i]-y0[i])
    end
    f
end

function force_multiple(d::DMP,x::AbstractMatrix)
    # Trajectory case
    T,n = size(x)
    y0 = _1(d)
    f = zeros(T,n)
    for i = 1:n
        Ψ    = kernel_matrix(x[:,i], d.c[:,i], d.σ2[:,i])
        f[:,i] = Ψ*d.w[:,i] .* (x[:,i].*(d.g[i]-y0[i])')
    end
    f
end

function acceleration(d::DMP, y::Number,ẏ::Number,x::Number,g::Number,i=1)
    f = force(d,x,i)
    d.opts.αz*(d.opts.βz*(g-y)-ẏ)+f
end

function acceleration(d::DMP, y::AbstractVector,ẏ::AbstractVector,x ,g::AbstractVector)
    f = force(d,x)
    d.opts.αz*(d.opts.βz*(g-y)-ẏ)+f
end

function acceleration(d::DMP, y::AbstractMatrix,ẏ::AbstractMatrix,x::AbstractVecOrMat,g::AbstractVector)
    f = force(d,x)
    d.opts.αz*(d.opts.βz*(g'.-y)-ẏ)+f
end

"""
`t,y,z,x = solve(dmp::DMP, t = 0:_T(dmp)-1; y0 = _1(dmp), g = dmp.g, solver=ode45)`

`t` time vector

## Keyword arguments: \n
`y0` start position, defaults to the initial point in training data from `dmp`
`g` goal, defaults to goal from `dmp`\n
`solver` the ode solver to use, see https://github.com/JuliaLang/ODE.jl \n

see also `plotdmp`
"""
function solve(dmp::DMP, t = dmp.t; y0 = _1(dmp), g = dmp.g, solver=ode45)
    if dmp.opts.sched_sig == :position
        return solve_position(dmp, t, y0, g, solver)
    elseif dmp.opts.sched_sig == :time
        return solve_time(dmp, t, y0, g, solver)
    end
    solve_canonical(dmp, t, y0, g, solver)
end

function solve_canonical(dmp, t, y0, g, solver)
    T,n     = size(dmp.y)
    αx      = dmp.opts.αx
    τ       = dmp.τ
    z       = zeros(T,n)
    y       = zeros(T,n)
    x       = zeros(T)
    for i = 1:n
        function time_derivative(t,state)
            local z   = state[1]
            local y   = state[2]
            local x   = state[3]
            zp  = acceleration(dmp, y, z, x,g[i],i)
            yp  = z
            xp  = -αx * x
            [zp;yp;xp] / τ
        end
        state0  = [0; y0[i]; 1.]
        tout,state_history = solver(time_derivative, state0, t,points=:specified)
        res = vv2m(state_history)
        z[:,i] = res[:,1]/τ
        y[:,i] = res[:,2]
        x = res[:,3] # TODO: se till att denna är samma för alla DOF
    end
    t,y,z,x
end

function solve_position(dmp, t, y0, g, solver)
    T,n     = size(dmp.y)
    αx      = dmp.opts.αx
    τ       = dmp.τ
    z       = zeros(T,n)
    y       = zeros(T,n)
    for i = 1:n
        function time_derivative(t,state)
            local z   = state[1]
            local y   = state[2]
            zp  = acceleration(dmp, y, z, g[i]-y,g[i],i)
            yp  = z
            [zp;yp] / τ
        end
        state0  = [0; y0[i]]
        tout,state_history = solver(time_derivative, state0, t,points=:specified)
        res = vv2m(state_history)
        z[:,i] = res[:,1]/τ
        y[:,i] = res[:,2]
    end
    t,y,z,y
end

function solve_time(dmp, t, y0, g, solver)
    T,n     = size(dmp.y)
    αx      = dmp.opts.αx
    τ       = dmp.τ
    z       = zeros(T,n)
    y       = zeros(T,n)
    for i = 1:n
        function time_derivative(t,state)
            local z   = state[1]
            local y   = state[2]
            zp  = acceleration(dmp, y, z, t ,g[i],i)
            yp  = z
            [zp;yp] / τ
        end
        state0  = [0; y0[i]]
        tout,state_history = solver(time_derivative, state0, t,points=:specified)
        res = vv2m(state_history)
        z[:,i] = res[:,1]/τ
        y[:,i] = res[:,2]
    end
    t,y,z,t
end

plotdmp(dmp::DMP, args...) = println("To plot a DMP, install package Plots.jl or call tout,yout,ẏout,xout = solve(dmp) to produce your own plot.")

@require Plots begin
math(sl) = map(s->string("\$",s,"\$") ,sl)
function plotdmp(dmp::DMP; kwargs...)
    tout,yout,ẏout,xout = solve(dmp)
    n = size(dmp.y,2)
    fig = Plots.subplot(n = n, nc = 1)
    for i = 1:n
        Plots.plot!(fig[i,1],tout,[yout[:,i] ẏout[:,i]],lab = ["y_{out}" "ẏ_{out}"] |> math; kwargs...)
        Plots.plot!(fig[i,1],tout,[dmp.y[:,i] dmp.ẏ[:,i]],l=:dash,lab = ["y" "ẏ"] |> math; kwargs...)
    end
end

function plotdmpphase(dmp::DMP; kwargs...)
    tout,yout,ẏout,xout = solve(dmp)
    Plots.plot(yout[:,1],yout[:,2],lab = ["y_{out}" "ẏ_{out}"] |> math; kwargs...)
    Plots.plot!(dmp.y[:,1],dmp.y[:,2],l=:dash,lab = ["y" "ẏ"] |> math; kwargs...)
end
end

end # module

ext: jl


text:
import Statistics
betterThanAverage(classPoints, yourPoints) = yourPoints > Statistics.mean(classPoints)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
    CloudOptics

Provides cloud optical properties as a function of effective radius for the RRTMGP bands
  Based on Mie calculations for liquid
    and results from doi:10.1175/JAS-D-12-039.1 for ice with variable surface roughness
  Can use either look-up tables or Pade approximates according to which data has been loaded
  Mike Iacono (AER) is the original author
The class can be used as-is but is also intended as an example of how to extend the RTE framework
"""
module CloudOptics

using DocStringExtensions
using ..OpticalProps
using ..Utilities

export CloudOpticalProps, CloudOpticalPropsPGP
export cloud_optics!
export get_min_radius, get_max_radius

export PadeMethod, LookUpTable
export CloudOpticsLUT
export CloudOpticsPade

abstract type AbstractInterpolationMethod{FT <: AbstractFloat} end

"""
    PadeMethod{FT} <: AbstractInterpolationMethod{FT}

Pade approximation coefficients

# Fields
$(DocStringExtensions.FIELDS)
"""
struct PadeMethod{FT} <: AbstractInterpolationMethod{FT}
    "Extinction coefficients"
    ext::Array{FT}
    "Single scattering albedo"
    ssa::Array{FT}
    "Asymmetry parameter"
    asy::Array{FT}
    "Extinction coefficient particle size regime boundaries"
    sizreg_ext::Array{FT}
    "Single scattering albedo particle size regime boundaries"
    sizreg_ssa::Array{FT}
    "Asymmetry parameter particle size regime boundaries"
    sizreg_asy::Array{FT}
    "particle size lower bound for interpolation"
    rad_lwr::FT
    "particle size upper bound for interpolation"
    rad_upr::FT
    function PadeMethod(
        ext::Array{FT},
        ssa::Array{FT},
        asy::Array{FT},
        sizreg_ext::Array{FT},
        sizreg_ssa::Array{FT},
        sizreg_asy::Array{FT},
    ) where {FT}

        # Error checking
        nbnd = size(ext, 1)
        nsizereg = size(ext, 2)
        ncoeff_ssa_g = size(ssa, 3)
        nbound = length(sizreg_ext)
        nrghice = size(ext, 4)

        @assert nsizereg == 3
        if length(size(ext)) == 4
            @assert all(size(ssa) .== (nbnd, nsizereg, ncoeff_ssa_g, nrghice))
            @assert all(size(asy) .== (nbnd, nsizereg, ncoeff_ssa_g, nrghice))
        else
            @assert all(size(ssa) .== (nbnd, nsizereg, ncoeff_ssa_g))
            @assert all(size(asy) .== (nbnd, nsizereg, ncoeff_ssa_g))
        end
        @assert length(sizreg_ssa) == nbound
        @assert length(sizreg_asy) == nbound
        # Consistency among size regimes
        rad_lwr = sizreg_ext[1]
        rad_upr = sizreg_ext[end]
        @assert !any([sizreg_ssa[1], sizreg_asy[1]] .< rad_lwr)
        @assert !any([sizreg_ssa[nbound], sizreg_asy[nbound]] .> rad_upr)

        return new{FT}(ext, ssa, asy, sizreg_ext, sizreg_ssa, sizreg_asy, rad_lwr, rad_upr)
    end
end

"""
    LookUpTable{FT,I<:Int} <: AbstractInterpolationMethod{FT}

Lookup table interpolation constants and coefficients

# Fields
$(DocStringExtensions.FIELDS)
"""
struct LookUpTable{FT, I <: Int} <: AbstractInterpolationMethod{FT}
    "particle size lower bound for interpolation"
    rad_lwr::FT
    "particle size upper bound for interpolation"
    rad_upr::FT
    "constant for calculating interpolation indexes"
    rad_fac::FT
    "Number of evenly-spaced elements for linear interpolation"
    nsteps::I
    "Size of step for linear interpolation"
    step_size::FT
    "Extinction coefficients"
    ext::Array{FT}
    "Single scattering albedo"
    ssa::Array{FT}
    "Asymmetry parameter"
    asy::Array{FT}
    function LookUpTable(
        rad_lwr::FT,
        rad_upr::FT,
        rad_fac::FT,
        ext::Array{FT},
        ssa::Array{FT},
        asy::Array{FT},
    ) where {FT <: AbstractFloat}

        # Error checking
        # LUT coefficient dimensions
        nsize = size(ext, 1)
        nbnd = size(ext, 2)

        # Error checking
        #   Can we check for consistency between table bounds and _fac?
        if length(size(ext)) == 3
            nrghice = size(ext, 3)
            @assert all(size(ssa) .== (nsize, nbnd, nrghice))
            @assert all(size(asy) .== (nsize, nbnd, nrghice))
        else
            @assert all(size(ssa) .== (nsize, nbnd))
            @assert all(size(asy) .== (nsize, nbnd))
        end

        nsteps = size(ext, 1)
        step_size = (rad_upr - rad_lwr) / FT(nsteps - 1)

        return new{FT, Int}(rad_lwr, rad_upr, rad_fac, nsteps, step_size, ext, ssa, asy)
    end
end
get_min_radius(aim::AbstractInterpolationMethod) = aim.rad_lwr
get_max_radius(aim::AbstractInterpolationMethod) = aim.rad_upr

get_num_roughness_types(this::PadeMethod) = size(this.ext, 4)
get_num_roughness_types(this::LookUpTable) = size(this.ext, 3)

"""
    CloudOpticsLUT{FT, I} <: AbstractOpticalProps{FT, I}

# Fields
$(DocStringExtensions.FIELDS)
"""
struct CloudOpticsLUT{FT, I} <: AbstractOpticalProps{FT, I}
    base::OpticalPropsBase{FT, I}
    "Ice surface roughness category - needed for Yang (2013) ice optics parameterization"
    icergh::I
    "Interpolation method for liquid properties"
    liq::LookUpTable{FT, I}
    "Interpolation method for ice properties"
    ice::LookUpTable{FT, I}
end

"""
    CloudOpticsPade{FT, I} <: AbstractOpticalProps{FT, I}

# Fields
$(DocStringExtensions.FIELDS)
"""
struct CloudOpticsPade{FT, I} <: AbstractOpticalProps{FT, I}
    "Base optical properties"
    base::OpticalPropsBase{FT, I}
    "Ice surface roughness category - needed for Yang (2013) ice optics parameterization"
    icergh::I
    "Interpolation method for liquid properties"
    liq::PadeMethod{FT}
    "Interpolation method for ice properties"
    ice::PadeMethod{FT}
end

"""
    CloudOpticalProps{FT<:AbstractFloat}

Cloud properties (liquid or ice) used to combine with
averaged optical properties

# Fields
$(DocStringExtensions.FIELDS)
"""
struct CloudOpticalProps{FT <: AbstractFloat}
    "Water path"
    wp::Array{FT}
    "Effective radius"
    re::Array{FT}
end
CloudOpticalProps(FT, ncol, nlay) = CloudOpticalProps{FT}(zeros(FT, ncol, nlay), zeros(FT, ncol, nlay))

"""
    CloudOpticalPropsPGP{FT<:AbstractFloat}

Cloud properties (liquid or ice) used to combine with
averaged optical properties per grid point

# Fields
$(DocStringExtensions.FIELDS)
"""
mutable struct CloudOpticalPropsPGP{FT <: AbstractFloat}
    "Water path"
    wp::FT
    "Effective radius"
    re::FT
end
CloudOpticalPropsPGP(FT) = CloudOpticalPropsPGP{FT}(FT(0), FT(0))

function Base.convert(::Type{CloudOpticalProps}, data::Array{CloudOpticalPropsPGP{FT}}) where {FT}
    s = size(data)
    return CloudOpticalProps{FT}(
        Array{FT}([data[i, j].wp for i in 1:s[1], j in 1:s[2]]),
        Array{FT}([data[i, j].re for i in 1:s[1], j in 1:s[2]]),
    )
end

Base.convert(::Type{Array{CloudOpticalPropsPGP}}, data::CloudOpticalProps{FT}) where {FT} =
    [CloudOpticalPropsPGP{FT}(data.wp[i, j], data.re[i, j]) for i in 1:size(data.wp, 1), j in 1:size(data.wp, 2)]

#####
##### Derive cloud optical properties from provided cloud physical properties
#####

"""
    combine_optical_props!(op::AbstractOpticalPropsArry, liq::TwoStream{FT}, ice::TwoStream{FT}) where {FT<:AbstractFloat}

Combine liquid and ice contributions into total cloud optical properties
   See also the `increment!` routines in `OpticalProps_kernels`
"""
function combine_optical_props!(op::OneScalar{FT}, liq::TwoStream{FT}, ice::TwoStream{FT}) where {FT <: AbstractFloat}
    # Absorption optical depth  = (1-ssa) * τ = τ - τssa
    nbnd = size(liq.τ, 3)
    op.τ[:, :, 1:nbnd] .= (liq.τ .- liq.ssa) + (ice.τ .- ice.ssa)
end
function combine_optical_props!(op::TwoStream{FT}, liq::TwoStream{FT}, ice::TwoStream{FT}) where {FT <: AbstractFloat}
    nbnd = size(liq.τ, 3)
    τ = liq.τ + ice.τ
    ssa = liq.ssa + ice.ssa
    op.g[:, :, 1:nbnd] .= (liq.g + ice.g) ./ max.(eps(FT), ssa)
    op.ssa[:, :, 1:nbnd] .= ssa ./ max.(eps(FT), τ)
    op.τ[:, :, 1:nbnd] .= τ
end
function combine_optical_props!(
    op::OneScalarPGP{FT},
    liq::TwoStreamPGP{FT},
    ice::TwoStreamPGP{FT},
) where {FT <: AbstractFloat}
    # Absorption optical depth  = (1-ssa) * τ = τ - τssa
    nbnd = size(liq.τ, 1)
    op.τ[1:nbnd] .= (liq.τ .- liq.ssa) + (ice.τ .- ice.ssa)
end
function combine_optical_props!(
    op::TwoStreamPGP{FT},
    liq::TwoStreamPGP{FT},
    ice::TwoStreamPGP{FT},
) where {FT <: AbstractFloat}
    nbnd = size(liq.τ, 1)
    τ = liq.τ + ice.τ
    ssa = liq.ssa + ice.ssa
    op.g[1:nbnd] .= (liq.g + ice.g) ./ max.(eps(FT), ssa)
    op.ssa[1:nbnd] .= ssa ./ max.(eps(FT), τ)
    op.τ[1:nbnd] .= τ
end

function validate_cloud_optics!(
    this::AbstractOpticalProps{FT},
    liq::CloudOpticalProps{FT},
    ice::CloudOpticalProps{FT},
    optical_props::AbstractOpticalPropsArry{FT},
) where {FT}
    liqmsk = BitArray(liq.wp .> FT(0))
    icemsk = BitArray(ice.wp .> FT(0))
    # Error checking
    @assert bands_are_equal(this, optical_props)
    @assert get_nband(optical_props) == get_ngpt(optical_props)
    @assert !(this.icergh < 1 || this.icergh > get_num_roughness_types(this.ice))
    @assert !any_vals_outside(liq.re, liqmsk, this.liq.rad_lwr, this.liq.rad_upr)
    @assert !any_vals_outside(ice.re, icemsk, this.ice.rad_lwr, this.ice.rad_upr)
    any(liqmsk) && @assert !any_vals_less_than(liq.wp, liqmsk, FT(0))
    any(icemsk) && @assert !any_vals_less_than(ice.wp, icemsk, FT(0))
    return nothing
end
function validate_cloud_optics!(
    this::AbstractOpticalProps{FT},
    liq::CloudOpticalPropsPGP{FT},
    ice::CloudOpticalPropsPGP{FT},
    optical_props::AbstractOpticalPropsPGP{FT},
) where {FT}
    liqmsk = liq.wp > FT(0)
    icemsk = ice.wp > FT(0)
    # Error checking
    @assert bands_are_equal(this, optical_props)
    @assert get_nband(optical_props) == get_ngpt(optical_props)
    @assert !(this.icergh < 1 || this.icergh > get_num_roughness_types(this.ice))
    @assert !any_vals_outside(liq.re, liqmsk, this.liq.rad_lwr, this.liq.rad_upr)
    @assert !any_vals_outside(ice.re, icemsk, this.ice.rad_lwr, this.ice.rad_upr)
    liqmsk && @assert !any_vals_less_than(liq.wp, liqmsk, FT(0))
    icemsk && @assert !any_vals_less_than(ice.wp, icemsk, FT(0))
    return nothing
end

"""
    cloud_optics!(this::AbstractOpticalProps{FT},
                  clouds_liq::CloudOpticalProps{FT},
                  clouds_ice::CloudOpticalProps{FT},
                  optical_props::AbstractOpticalPropsArry{FT}) where {FT<:AbstractFloat}

Compute single-scattering properties in

 - `optical_props` optical properties, see [`AbstractOpticalPropsArry`](@ref)

given

 - `this` cloud optics, see [`AbstractOpticalProps`](@ref)
 - `clouds_liq` cloud properties for liquid, see [`CloudOpticalProps`](@ref)
 - `clouds_ice` cloud properties for ice   , see [`CloudOpticalProps`](@ref)
"""
function cloud_optics!(
    this::CloudOpticsPade{FT},
    clouds_liq::CloudOpticalProps{FT},
    clouds_ice::CloudOpticalProps{FT},
    optical_props::AbstractOpticalPropsArry{FT},
) where {FT <: AbstractFloat}

    nbnd = get_nband(this)

    validate_cloud_optics!(this, clouds_liq, clouds_ice, optical_props)

    liq = TwoStream(FT, size(clouds_liq.wp)..., nbnd)
    ice = TwoStream(FT, size(clouds_ice.wp)..., nbnd)

    #### Compute cloud optical properties.

    compute_all_from_pade!(nbnd, clouds_liq, this.liq, liq)              # Liquid
    compute_all_from_pade!(nbnd, clouds_ice, this.ice, ice, this.icergh) # Ice

    # Copy total cloud properties onto outputs
    combine_optical_props!(optical_props, liq, ice)

end
function cloud_optics!(
    this::CloudOpticsLUT{FT},
    clouds_liq::CloudOpticalProps{FT},
    clouds_ice::CloudOpticalProps{FT},
    optical_props::AbstractOpticalPropsArry{FT},
) where {FT}

    nbnd = get_nband(this)

    validate_cloud_optics!(this, clouds_liq, clouds_ice, optical_props)

    liq = TwoStream(FT, size(clouds_liq.wp)..., nbnd)
    ice = TwoStream(FT, size(clouds_ice.wp)..., nbnd)

    #### Compute cloud optical properties.

    compute_all_from_table!(nbnd, clouds_liq, this.liq, liq)              # Liquid
    compute_all_from_table!(nbnd, clouds_ice, this.ice, ice, this.icergh) # Ice

    # Copy total cloud properties onto outputs
    combine_optical_props!(optical_props, liq, ice)

end
function cloud_optics!(
    this::CloudOpticsPade{FT},
    clouds_liq::CloudOpticalPropsPGP{FT},
    clouds_ice::CloudOpticalPropsPGP{FT},
    optical_props::AbstractOpticalPropsPGP{FT},
) where {FT <: AbstractFloat}

    nbnd = get_nband(this)

    validate_cloud_optics!(this, clouds_liq, clouds_ice, optical_props)

    liq = TwoStreamPGP(FT, nbnd)
    ice = TwoStreamPGP(FT, nbnd)

    #### Compute cloud optical properties.

    compute_all_from_pade!(nbnd, clouds_liq, this.liq, liq)              # Liquid
    compute_all_from_pade!(nbnd, clouds_ice, this.ice, ice, this.icergh) # Ice

    # Copy total cloud properties onto outputs
    combine_optical_props!(optical_props, liq, ice)

end
function cloud_optics!(
    this::CloudOpticsLUT{FT},
    clouds_liq::CloudOpticalPropsPGP{FT},
    clouds_ice::CloudOpticalPropsPGP{FT},
    optical_props::AbstractOpticalPropsPGP{FT},
) where {FT}

    nbnd = get_nband(this)

    validate_cloud_optics!(this, clouds_liq, clouds_ice, optical_props)

    liq = TwoStreamPGP(FT, nbnd)
    ice = TwoStreamPGP(FT, nbnd)

    #### Compute cloud optical properties.

    compute_all_from_table!(nbnd, clouds_liq, this.liq, liq)              # Liquid
    compute_all_from_table!(nbnd, clouds_ice, this.ice, ice, this.icergh) # Ice

    # Copy total cloud properties onto outputs
    combine_optical_props!(optical_props, liq, ice)

end

"""
    compute_all_from_table!(nbnd::I,
                            clouds::CloudOpticalProps{FT},
                            lut::LookUpTable{FT,I},
                            op::TwoStream{FT,I},
                            icergh::Union{Nothing,I}=nothing) where {FT<:AbstractFloat,I<:Int}

Linearly interpolate values from a lookup table with "nsteps" evenly-spaced
elements starting at "offset." in

 - `op` optical properties, see [`TwoStream`](@ref)

given

 - `nbnd` number of bands
 - `clouds` cloud optical properties, see [`CloudOpticalProps`](@ref)
 - `lut` interpolation method (look-up table), see [`LookUpTable`](@ref)
 - `icergh` ice surface roughness category - needed for Yang (2013) ice optics parameterization

The table's second dimension is band. Returns 0 where `clouds.wp ≤ 0`.
"""
function compute_all_from_table!(
    nbnd::I,
    clouds::CloudOpticalProps{FT},
    lut::LookUpTable{FT, I},
    op::TwoStream{FT, I},
    icergh::Union{Nothing, I} = nothing,
) where {FT <: AbstractFloat, I <: Int}

    ncol, nlay = size(clouds.re)

    offset = lut.rad_lwr

    τ_table = icergh ≠ nothing ? lut.ext[:, :, icergh] : lut.ext
    ssa_table = icergh ≠ nothing ? lut.ssa[:, :, icergh] : lut.ssa
    asy_table = icergh ≠ nothing ? lut.asy[:, :, icergh] : lut.asy

    @inbounds for ibnd in 1:nbnd
        @inbounds for ilay in 1:nlay
            @inbounds for icol in 1:ncol
                if clouds.wp[icol, ilay] > FT(0)
                    index =
                        convert(Int, min(floor((clouds.re[icol, ilay] - offset) / lut.step_size) + 1, lut.nsteps - 1))
                    fint = (clouds.re[icol, ilay] - offset) / lut.step_size - (index - 1)
                    τ =
                        clouds.wp[icol, ilay] *
                        (τ_table[index, ibnd] + fint * (τ_table[index + 1, ibnd] - τ_table[index, ibnd]))
                    τs = τ * (ssa_table[index, ibnd] + fint * (ssa_table[index + 1, ibnd] - ssa_table[index, ibnd]))
                    op.g[icol, ilay, ibnd] =
                        τs * (asy_table[index, ibnd] + fint * (asy_table[index + 1, ibnd] - asy_table[index, ibnd]))
                    op.ssa[icol, ilay, ibnd] = τs
                    op.τ[icol, ilay, ibnd] = τ
                else
                    op.τ[icol, ilay, ibnd] = FT(0)
                    op.ssa[icol, ilay, ibnd] = FT(0)
                    op.g[icol, ilay, ibnd] = FT(0)
                end
            end
        end
    end
end
function compute_all_from_table!(
    nbnd::I,
    clouds::CloudOpticalPropsPGP{FT},
    lut::LookUpTable{FT, I},
    op::TwoStreamPGP{FT, I},
    icergh::Union{Nothing, I} = nothing,
) where {FT <: AbstractFloat, I <: Int}

    offset = lut.rad_lwr

    τ_table = icergh ≠ nothing ? lut.ext[:, :, icergh] : lut.ext
    ssa_table = icergh ≠ nothing ? lut.ssa[:, :, icergh] : lut.ssa
    asy_table = icergh ≠ nothing ? lut.asy[:, :, icergh] : lut.asy

    @inbounds for ibnd in 1:nbnd
        if clouds.wp > FT(0)
            index = convert(Int, min(floor((clouds.re - offset) / lut.step_size) + 1, lut.nsteps - 1))
            fint = (clouds.re - offset) / lut.step_size - (index - 1)
            τ = clouds.wp * (τ_table[index, ibnd] + fint * (τ_table[index + 1, ibnd] - τ_table[index, ibnd]))
            τs = τ * (ssa_table[index, ibnd] + fint * (ssa_table[index + 1, ibnd] - ssa_table[index, ibnd]))
            op.g[ibnd] = τs * (asy_table[index, ibnd] + fint * (asy_table[index + 1, ibnd] - asy_table[index, ibnd]))
            op.ssa[ibnd] = τs
            op.τ[ibnd] = τ
        else
            op.τ[ibnd] = FT(0)
            op.ssa[ibnd] = FT(0)
            op.g[ibnd] = FT(0)
        end
    end
end

#####
##### Pade functions
#####

"""
    compute_all_from_pade!(nbnd::I,
                           clouds::CloudOpticalProps{FT},
                           pm::PadeMethod{FT},
                           op::TwoStream{FT,I},
                           icergh::Union{Nothing,I}=nothing) where {FT<:AbstractFloat,I<:Int}

 - `nbnd` number of bands
 - `clouds` cloud optical properties, see [`CloudOpticalProps`](@ref)
 - `pm` interpolation method (Pade), see [`PadeMethod`](@ref)
 - `op` optical properties, see [`TwoStream`](@ref)
 - `icergh` ice surface roughness category - needed for Yang (2013) ice optics parameterization
"""
function compute_all_from_pade!(
    nbnd::I,
    clouds::CloudOpticalProps{FT},
    pm::PadeMethod{FT},
    op::TwoStream{FT, I},
    icergh::Union{Nothing, I} = nothing,
) where {FT <: AbstractFloat, I <: Int}

    nsizes = size(pm.ext, 2)
    # Cloud optical properties from Pade coefficient method
    #   Hard coded assumptions: order of approximants, three size regimes
    m_ext, n_ext = 2, 3
    m_ssa, n_ssa = 2, 2
    m_asy, n_asy = 2, 2

    coeffs_ext = icergh ≠ nothing ? pm.ext[:, :, :, icergh] : pm.ext
    coeffs_ssa = icergh ≠ nothing ? pm.ssa[:, :, :, icergh] : pm.ssa
    coeffs_asy = icergh ≠ nothing ? pm.asy[:, :, :, icergh] : pm.asy

    ncol, nlay = size(clouds.re)
    @inbounds for ibnd in 1:nbnd
        @inbounds for ilay in 1:nlay
            @inbounds for icol in 1:ncol
                if clouds.wp[icol, ilay] > FT(0)
                    re = clouds.re[icol, ilay]
                    #
                    # Finds index into size regime table
                    # This works only if there are precisely three size regimes (four bounds) and it's
                    #   previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
                    #
                    irad = convert(Int, min(floor((re - pm.sizreg_ext[2]) / pm.sizreg_ext[3]) + 2, 3))
                    τ = clouds.wp[icol, ilay] * pade_eval(ibnd, nbnd, nsizes, m_ext, n_ext, irad, re, coeffs_ext)

                    irad = convert(Int, min(floor((re - pm.sizreg_ssa[2]) / pm.sizreg_ssa[3]) + 2, 3))
                    # Pade approximants for co-albedo can sometimes be negative
                    τs = τ * (FT(1) - max(FT(0), pade_eval(ibnd, nbnd, nsizes, m_ssa, n_ssa, irad, re, coeffs_ssa)))
                    irad = convert(Int, min(floor((re - pm.sizreg_asy[2]) / pm.sizreg_asy[3]) + 2, 3))
                    op.g[icol, ilay, ibnd] = τs * pade_eval(ibnd, nbnd, nsizes, m_asy, n_asy, irad, re, coeffs_asy)
                    op.ssa[icol, ilay, ibnd] = τs
                    op.τ[icol, ilay, ibnd] = τ
                else
                    op.τ[icol, ilay, ibnd] = FT(0)
                    op.ssa[icol, ilay, ibnd] = FT(0)
                    op.g[icol, ilay, ibnd] = FT(0)
                end
            end
        end
    end
end
function compute_all_from_pade!(
    nbnd::I,
    clouds::CloudOpticalPropsPGP{FT},
    pm::PadeMethod{FT},
    op::TwoStreamPGP{FT, I},
    icergh::Union{Nothing, I} = nothing,
) where {FT <: AbstractFloat, I <: Int}

    nsizes = size(pm.ext, 2)
    # Cloud optical properties from Pade coefficient method
    #   Hard coded assumptions: order of approximants, three size regimes
    m_ext, n_ext = 2, 3
    m_ssa, n_ssa = 2, 2
    m_asy, n_asy = 2, 2

    coeffs_ext = icergh ≠ nothing ? pm.ext[:, :, :, icergh] : pm.ext
    coeffs_ssa = icergh ≠ nothing ? pm.ssa[:, :, :, icergh] : pm.ssa
    coeffs_asy = icergh ≠ nothing ? pm.asy[:, :, :, icergh] : pm.asy

    @inbounds for ibnd in 1:nbnd
        if clouds.wp > FT(0)
            re = clouds.re
            #
            # Finds index into size regime table
            # This works only if there are precisely three size regimes (four bounds) and it's
            #   previously guaranteed that size_bounds(1) <= size <= size_bounds(4)
            #
            irad = convert(Int, min(floor((re - pm.sizreg_ext[2]) / pm.sizreg_ext[3]) + 2, 3))
            τ = clouds.wp * pade_eval(ibnd, nbnd, nsizes, m_ext, n_ext, irad, re, coeffs_ext)

            irad = convert(Int, min(floor((re - pm.sizreg_ssa[2]) / pm.sizreg_ssa[3]) + 2, 3))
            # Pade approximants for co-albedo can sometimes be negative
            τs = τ * (FT(1) - max(FT(0), pade_eval(ibnd, nbnd, nsizes, m_ssa, n_ssa, irad, re, coeffs_ssa)))
            irad = convert(Int, min(floor((re - pm.sizreg_asy[2]) / pm.sizreg_asy[3]) + 2, 3))
            op.g[ibnd] = τs * pade_eval(ibnd, nbnd, nsizes, m_asy, n_asy, irad, re, coeffs_asy)
            op.ssa[ibnd] = τs
            op.τ[ibnd] = τ
        else
            op.τ[ibnd] = FT(0)
            op.ssa[ibnd] = FT(0)
            op.g[ibnd] = FT(0)
        end
    end
end

#####
##### Ancillary functions
#####

"""
    pade_eval(iband::I,
              nbnd::I,
              nrads::I,
              m::I,
              n::I,
              irad::I,
              re::FT,
              pade_coeffs::AbstractArray{FT}) where {FT<:AbstractFloat, I<:Int}

Evaluate Pade approximant of order [m/n], given

 - `pade_coeffs` pade coefficients [nbnd, nrads, 1:m+n+1]
"""
function pade_eval(
    iband::I,
    nbnd::I,
    nrads::I,
    m::I,
    n::I,
    irad::I,
    re::FT,
    pade_coeffs::AbstractArray{FT},
) where {FT <: AbstractFloat, I <: Int}

    denom = pade_coeffs[iband, irad, n + m + 1]
    @inbounds for i in (n - 1 + m):-1:(1 + m)
        denom = pade_coeffs[iband, irad, i + 1] + re * denom
    end
    denom = 1 + re * denom

    numer = pade_coeffs[iband, irad, m + 1]
    @inbounds for i in (m - 1):-1:1
        numer = pade_coeffs[iband, irad, i + 1] + re * numer
    end
    numer = pade_coeffs[iband, irad, 1] + re * numer

    pade_eval_1 = numer / denom
    return pade_eval_1
end

end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
VERSION

# %%
function G(n)
    for i1 in 0:n-1, i2 in 0:n-1, i3 in 0:n-1, i4 in 0:n-1, i5 in 0:n-1
    end
end

@code_native debuginfo=:none G(20)

# %%
function H(n)
    for i in 0:n^5-1
    end
end

@code_native debuginfo=:none H(20)

# %%
function f(n)
    s = 0.0
    for (i1, i2, i3, i4, i5) in Iterators.product(0:n-1, 0:n-1, 0:n-1, 0:n-1, 0:n-1)
        x = evalpoly(n, (i5, i4, i3, i2, i1)) + 1
        s += sin(x)/x
    end
    2s + 1
end

function g(n)
    s = 0.0
    for i1 in 0:n-1, i2 in 0:n-1, i3 in 0:n-1, i4 in 0:n-1, i5 in 0:n-1
        x = evalpoly(n, (i1, i2, i3, i4, i5)) + 1
        s += sin(x)/x
    end
    2s + 1
end

function h(n)
    s = 0.0
    for i in 0:n^5-1
        x = float(i) + 1
        s += sin(x)/x
    end
    2s + 1
end

using BenchmarkTools
@show VERSION
@show f(20)
@show g(20)
@show h(20)
@show f(20) ≈ g(20) ≈ h(20)
@btime f(20)
@btime g(20)
@btime h(20);

# %%
function f(n)
    s = 0.0
    for (i1, i2, i3, i4, i5) in Iterators.product(0:n-1, 0:n-1, 0:n-1, 0:n-1, 0:n-1)
        x = evalpoly(n, (i5, i4, i3, i2, i1)) + 1
        s += sin(x)/x
    end
    2s + 1
end

function g(n)
    s = 0.0
    for i1 in 0:n-1, i2 in 0:n-1, i3 in 0:n-1, i4 in 0:n-1, i5 in 0:n-1
        x = evalpoly(n, (i1, i2, i3, i4, i5)) + 1
        s += sin(x)/x
    end
    2s + 1
end

function h(n)
    s = 0.0
    for i in 0:n^5-1
        x = float(i) + 1
        s += sin(x)/x
    end
    2s + 1
end

using BenchmarkTools
@show VERSION
@show f(20)
@show g(20)
@show h(20)
@show f(20) ≈ g(20) ≈ h(20)
@btime f(20)
@btime g(20)
@btime h(20);

# %%

max_stars_repo_path: 0007/five-fold loop benchmark.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/five-fold loop benchmark.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/five-fold loop benchmark.jl

max_forks_repo_name: genkuroki/public


text:
# @osx_only begin
#     const qtlibdir = "/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/~/usr/lib/"
#     const QtCore = joinpath(qtlibdir,"QtCore.framework/")
#     const QtWidgets = joinpath(qtlibdir,"QtWidgets.framework/")

#     addHeaderDir(qtlibdir; isFramework = true, kind = C_System)

#     dlopen(joinpath(QtCore,"QtCore_debug"))
#     addHeaderDir(joinpath(QtCore,"Headers"), kind = C_System)
#     addHeaderDir("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore")

#     cxxinclude("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore/private/qcoreapplication_p.h")

#     dlopen(joinpath(QtWidgets,"QtWidgets"))
#     addHeaderDir(joinpath(QtWidgets,"Headers"), kind = C_System)
# end

using Cxx

@linux_only begin
  # TODO: probably allow for overrides from ENV

  const qtincdir = "/usr/include/qt5"
  addHeaderDir(qtincdir, kind = C_System)
  for dir in ("QtCore", "QtWidgets", "QtGui")
    addHeaderDir(joinpath(qtincdir, "$dir/"), kind = C_System)
  end

  const qtlibdir = "/usr/lib64"
  for lib in ("libQt5Core.so", "libQt5Gui.so", "libQt5Widgets.so")
    Libdl.dlopen(joinpath(qtlibdir, lib), Libdl.RTLD_GLOBAL)
  end
end

# basic includes
cxx"""
  #include <iostream>
  #include <QtCore>
  #include <QApplication>
  #include <QMessageBox>
  #include <QPushButton>
  // #include <QMainWindow>
  // #include <QWindow>
  // #include <QRasterWindow>
  #include <QBrush>
  #include <QPen>
  #include <QPixmap>
  #include <QWidget>
  #include <QPainter>
"""

type QApplication
  qapp
end

function QApplication()

  # create the args/nargs
  global const ac = [Int32(1)]
  global const a = "julia"
  x = Ptr{UInt8}[pointer(a),C_NULL]

  # create the QApplication
  global const app = QApplication(@cxxnew QApplication(*(pointer(ac)),pointer(x)))

  # quit on finalize  # TODO: make this work... not getting called
  # finalizer(app, quitapp)

  startEventLoop()
  app
end

function quitapp(app)
  @cxx (app.qapp)->quit()
end

# this will be called from the julia main loop
function update_loop(_timer)
    icxx"""
        $(app.qapp)->processEvents();
    """
end

# starts a timer which will call the QApplication's processEvents()
function startEventLoop()
  Base.Timer( update_loop, 0.1, 0.005 )
end


text:
push!(LOAD_PATH, "../src/")

using Documenter, IDSGraphs

makedocs(
    sitename="IDSGraphs.jl Documentation",
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    ),
    modules=[IDSGraphs],
    pages=[
        "Home" => "index.md",
        "Quick Start" => "guide_quickstart.md",
        "API Reference" => [
            "Structures" => "api_structures.md",
            "Component Dependency Graphs" => "api_depgraphs.md",
            "Internals" => "api_internals.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/tmthyln/IDSGraphs.jl.git",
    devbranch = "main",
    devurl="latest"
    )


text:
function fn(a)
    a[:] .= 9
end

function fn(addr, length)
    a = unsafe_wrap(Array{UInt32}, Ptr{UInt32}(addr), length)
    fn(a)
end

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
@symbol_func function P_BR(cur_reactor::AbstractReactor)
  cur_P = K_BR(cur_reactor)

  cur_P *= cur_reactor.R_0 ^ 3

  cur_P *= cur_reactor.n_bar ^ 2

  cur_P *= sqrt(cur_reactor.T_bar)

  cur_P
end


text:
using PyCall

function checkpip()
   # Import pip
   try
      pipver = `$(PyCall.python) -m pip --version`
      run(`$(pipver)`)
   catch
       # If it is not found, install it
       println("""I couldn't find `pip`. I will try to download and install it
                  automatically, but if this fails, please install
                  manually, then try to build `JuLIP` again.""")
       get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
       download("https://bootstrap.pypa.io/get-pip.py", get_pip)
       run(`$(PyCall.python) $get_pip --user`)
   end
end

function pip(pkgname)
   checkpip()
   pipcmd = `$(PyCall.python) -m pip install --upgrade --user $(pkgname)`
   run(`$(pipcmd)`)
end

println("Installing Dependencies of `ASE.jl`: `ase` and `matscipy`")

if is_unix()
   try
      @pyimport ase as _ase_
   catch
      println("""`ase` was not found, trying to install via pip. If this fails,
               please file an issue and try to install it manually, following
               the instructions at `https://wiki.fysik.dtu.dk/ase/install.html`""")
      pip("ase")
   end

   try
      @pyimport matscipy as _matscipy_
   catch
      println("""`matscipy` was not found, trying to install it. If this fails,
               please file an issue and try to install it manually, following
               the instructions at `https://github.com/libAtoms/matscipy`""")
      pip("matscipy")
   end
else
   println("""it looks like this is a windows machine? I don't dare try to
            automatically build the dependencies here -- sorry!
            If installing them by hand turns out non-trivial, please file an
            issue.""")
end

# This here is a gist with a slightly different approach to installing
# the dependencies, using @pyimport pip
# https://gist.github.com/Luthaf/368a23981c8ec095c3eb


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Interpolant-basis weighted sum-of-squares polynomial (of dimension `U`) epigraph
of `ℓ₁` norm (of dimension `R`) cone, parametrized by vector of matrices `Ps`
derived from interpolant basis and polynomial domain constraints.

    $(FUNCTIONNAME){T}(R::Int, U::Int, Ps::Vector{Matrix{T}}, use_dual::Bool = false)
"""
mutable struct WSOSInterpEpiNormOne{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    R::Int
    U::Int
    Ps::Vector{Matrix{T}}
    nu::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_prod_updated::Bool
    inv_hess_prod_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    mats::Vector{Vector{Matrix{T}}}
    matfact::Vector{Vector}
    ΛLi_Λ::Vector{Vector{Matrix{T}}}
    Λ11::Vector{Matrix{T}}
    tempΛ11::Vector{Matrix{T}}
    tempLL::Vector{Matrix{T}}
    tempLU::Vector{Matrix{T}}
    tempLU2::Vector{Matrix{T}}
    Λ11LiP::Vector{Matrix{T}}
    tempUU_vec::Vector{Matrix{T}} # reused in update_hess
    tempUU::Matrix{T}
    tempUU2::Matrix{T}
    tempURU::Matrix{T}
    tempURU2::Matrix{T}
    # for each (2, 2)-block pertaining to (lambda_1, lambda_i),
    # P * inv(Λ)[1, 1] * Ps = P * inv(Λ)i[2, 2] * Ps
    PΛiPs1::Vector{Vector{Matrix{T}}}
    # for each (2, 2)-block pertaining to (lambda_1, lambda_i),
    # P * inv(Λ)[2, 1] * Ps = P * inv(Λ)[1, 2]' * Ps
    PΛiPs2::Vector{Vector{Matrix{T}}}
    ΛLiPs11::Vector{Vector{Matrix{T}}}
    ΛLiPs12::Vector{Vector{Matrix{T}}}
    ΛLiP_dir11::Vector{Vector{Matrix{T}}}
    ΛLiP_dir12::Vector{Vector{Matrix{T}}}
    ΛLiP_dir21::Vector{Vector{Matrix{T}}}
    dder3_half::Vector{Vector{Matrix}}
    Λfact::Vector
    hess_edge_blocks::Vector{Matrix{T}}
    hess_diag_blocks::Vector{Matrix{T}}
    hess_diag_facts::Vector
    hess_diags::Vector{Matrix{T}}
    hess_schur_fact::Factorization{T}
    point_views::Vector
    Ps_times::Vector{Float64}
    Ps_order::Vector{Int}

    function WSOSInterpEpiNormOne{T}(
        R::Int,
        U::Int,
        Ps::Vector{Matrix{T}};
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert R >= 2
        for Pk in Ps
            @assert size(Pk, 1) == U
        end
        cone = new{T}()
        cone.use_dual_barrier = !use_dual # using dual barrier
        cone.dim = U * R
        cone.R = R
        cone.U = U
        cone.Ps = Ps
        cone.nu = R * sum(size(Pk, 2) for Pk in Ps)
        return cone
    end
end

function setup_extra_data!(cone::WSOSInterpEpiNormOne{T}) where {T <: Real}
    U = cone.U
    R = cone.R
    Ps = cone.Ps
    K = length(Ps)
    Ls = [size(Pk, 2) for Pk in cone.Ps]
    cone.mats = [[zeros(T, L, L) for _ in 1:(R - 1)] for L in Ls]
    # TODO preallocate better
    cone.matfact = [[cholesky(hcat([one(T)])) for _ in 1:R] for _ in cone.Ps]
    cone.hess_edge_blocks = [zeros(T, U, U) for _ in 1:(R - 1)]
    cone.hess_diag_blocks = [zeros(T, U, U) for _ in 1:R]
    # TODO preallocate better
    cone.hess_diag_facts = Any[cholesky(hcat([one(T)])) for _ in 1:(R - 1)]
    cone.hess_diags = [zeros(T, U, U) for _ in 1:R - 1]
    cone.ΛLi_Λ = [[zeros(T, L, L) for _ in 1:(R - 1)] for L in Ls]
    cone.Λ11 = [zeros(T, L, L) for L in Ls]
    cone.tempΛ11 = [zeros(T, L, L) for L in Ls]
    cone.tempLL = [zeros(T, L, L) for L in Ls]
    cone.tempLU = [zeros(T, L, U) for L in Ls]
    cone.tempLU2 = [zeros(T, L, U) for L in Ls]
    cone.Λ11LiP = [zeros(T, L, U) for L in Ls]
    cone.tempUU_vec = [zeros(T, U, U) for _ in Ps]
    cone.tempUU = zeros(T, U, U)
    cone.tempUU2 = zeros(T, U, U)
    cone.tempURU = zeros(T, U, U * (R - 1))
    cone.tempURU2 = zeros(T, U * (R - 1), U)
    cone.PΛiPs1 = [Vector{Matrix{T}}(undef, R) for Pk in Ps]
    cone.PΛiPs2 = [Vector{Matrix{T}}(undef, R) for Pk in Ps]
    cone.PΛiPs1 = [[zeros(T, U, U) for _ in 1:(R - 1)] for Pk in Ps]
    cone.PΛiPs2 = [[zeros(T, U, U) for _ in 1:(R - 1)] for Pk in Ps]
    cone.ΛLiPs11 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.ΛLiPs12 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.ΛLiP_dir11 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.ΛLiP_dir12 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.ΛLiP_dir21 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.dder3_half = [[zeros(T, 2 * L, 2 * U) for _ in 1:(R - 1)] for L in Ls]
    cone.Λfact = Vector{Any}(undef, K)
    cone.point_views = [view(cone.point, block_idxs(U, i)) for i in 1:R]
    cone.Ps_times = zeros(K)
    cone.Ps_order = collect(1:K)
    return cone
end

reset_data(cone::WSOSInterpEpiNormOne) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated =
    cone.hess_prod_updated = cone.inv_hess_prod_updated = false)

function set_initial_point!(arr::AbstractVector, cone::WSOSInterpEpiNormOne)
    @views arr[1:cone.U] .= 1
    @views arr[(cone.U + 1):end] .= 0
    return arr
end

function update_feas(cone::WSOSInterpEpiNormOne)
    @assert !cone.feas_updated
    U = cone.U
    R = cone.R
    Λfact = cone.Λfact
    point_views = cone.point_views

    # order the Ps by how long it takes to check feasibility, to improve efficiency
    sortperm!(cone.Ps_order, cone.Ps_times, initialized = true) # stochastic

    cone.is_feas = true
    for k in cone.Ps_order
        cone.Ps_times[k] = @elapsed @inbounds begin
            Pk = cone.Ps[k]
            Λ11j = cone.Λ11[k]
            tempΛ11j = cone.tempΛ11[k]
            LLk = cone.tempLL[k]
            LUk = cone.tempLU[k]
            ΛLi_Λ = cone.ΛLi_Λ[k]
            matsk = cone.mats[k]
            factk = cone.matfact[k]

            # first lambda
            @. LUk = Pk' * point_views[1]'
            mul!(tempΛ11j, LUk, Pk)
            copyto!(Λ11j, tempΛ11j)
            Λfact[k] = cholesky!(Symmetric(tempΛ11j, :U), check = false)
            if !isposdef(Λfact[k])
                cone.is_feas = false
                break
            end

            uo = U + 1
            @inbounds for r in 2:R
                r1 = r - 1
                matr = matsk[r1]
                factr = factk[r1]
                @. LUk = Pk' * point_views[r]'
                mul!(LLk, LUk, Pk)

                ldiv!(ΛLi_Λ[r1], Λfact[k].L, LLk)
                copyto!(matr, Λ11j)
                mul!(matr, ΛLi_Λ[r1]', ΛLi_Λ[r1], -1, 1)

                factk[r1] = cholesky!(Symmetric(matr, :U), check = false)
                if !isposdef(factk[r1])
                    cone.is_feas = false
                    cone.feas_updated = true
                    return cone.is_feas
                end
                uo += U
            end
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(cone::WSOSInterpEpiNormOne)
    @assert cone.is_feas
    U = cone.U
    R = cone.R
    Λfact = cone.Λfact
    matfact = cone.matfact
    grad = cone.grad
    grad .= 0

    @inbounds for k in eachindex(cone.Ps)
        Pk = cone.Ps[k]
        Λ11LiP = cone.Λ11LiP[k]
        ΛLi_Λ = cone.ΛLi_Λ[k]
        factk = cone.matfact[k]
        ΛLiPs11 = cone.ΛLiPs11[k]
        ΛLiPs12 = cone.ΛLiPs12[k]

        # prep ΛLiPs
        ldiv!(Λ11LiP, cone.Λfact[k].L, Pk')

        for r in 1:(R - 1)
            mul!(ΛLiPs12[r], ΛLi_Λ[r]', Λ11LiP, -1, false)
            ldiv!(factk[r].L, ΛLiPs12[r])
            ldiv!(ΛLiPs11[r], factk[r].L, Pk')
        end

        @views for u in 1:U
            grad[u] -= sum(abs2, Λ11LiP[:, u]) * R
            idx = u + U
            for r in 1:(R - 1)
                grad[u] -= 2 * sum(abs2, ΛLiPs12[r][:, u])
                grad[idx] -= 2 * dot(ΛLiPs12[r][:, u], ΛLiPs11[r][:, u])
                idx += U
            end
        end
    end

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::WSOSInterpEpiNormOne)
    cone.hess_prod_updated || update_hess_prod(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    U = cone.U
    H .= 0

    @. @views H[1:U, 1:U] = cone.hess_diag_blocks[1]
    @inbounds for r in 1:(cone.R - 1)
        idxs = block_idxs(U, r + 1)
        LinearAlgebra.copytri!(cone.hess_edge_blocks[r], 'U')
        @. @views H[1:U, idxs] = cone.hess_edge_blocks[r]
        @. @views H[idxs, idxs] = cone.hess_diag_blocks[r + 1]
    end

    cone.hess_updated = true
    return cone.hess
end

function update_hess_prod(cone::WSOSInterpEpiNormOne)
    @assert cone.grad_updated
    U = cone.U
    R = cone.R
    R2 = R - 2
    diag_blocks = cone.hess_diag_blocks
    edge_blocks = cone.hess_edge_blocks

    @inbounds for r in 1:(R - 1)
        diag_blocks[r] .= 0
        edge_blocks[r] .= 0
    end
    diag_blocks[R] .= 0

    @inbounds for k in eachindex(cone.Ps)
        Λ11LiP = cone.Λ11LiP[k]
        PΛ11iP = cone.tempUU_vec[k]
        PΛiPs1 = cone.PΛiPs1[k]
        PΛiPs2 = cone.PΛiPs2[k]
        ΛLiPs11 = cone.ΛLiPs11[k]
        ΛLiPs12 = cone.ΛLiPs12[k]

        # P * inv(Λ_11) * P' for (1, 1) hessian block and adding to PΛiPs[r][r]
        mul!(PΛ11iP, Λ11LiP', Λ11LiP)
        for r in 1:(R - 1)
            mul!(PΛiPs2[r], ΛLiPs12[r]', ΛLiPs11[r])
            mul!(PΛiPs1[r], ΛLiPs12[r]', ΛLiPs12[r])
            @. PΛiPs1[r] += PΛ11iP
            for j in 1:U, i in 1:j
                ij1 = PΛiPs1[r][i, j]
                ij2 = PΛiPs2[r][i, j]
                uu = 2 * (abs2(ij1) + abs2(ij2))
                diag_blocks[1][i, j] += uu
                diag_blocks[r + 1][i, j] += uu
                edge_blocks[r][i, j] += 4 * (ij1 * ij2)
            end
        end
        for j in 1:U, i in 1:j
            diag_blocks[1][i, j] -= abs2(PΛ11iP[i, j]) * R2
        end
    end

    cone.hess_prod_updated = true
    return prod
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::WSOSInterpEpiNormOne,
    )
    if !cone.hess_prod_updated
        update_hess_prod(cone)
    end
    U = cone.U
    R = cone.R
    prod .= 0
    diag_blocks = cone.hess_diag_blocks

    @views mul!(prod[1:U, :], Symmetric(diag_blocks[1], :U), arr[1:U, :])
    @inbounds @views for r in 1:(R - 1)
        idxs = block_idxs(U, r + 1)
        edge_r = Symmetric(cone.hess_edge_blocks[r], :U)
        arr_r = arr[idxs, :]
        mul!(prod[1:U, :], edge_r, arr_r, true, true)
        mul!(prod[idxs, :], edge_r, arr[1:U, :])
        mul!(prod[idxs, :], Symmetric(diag_blocks[r + 1], :U), arr_r, true, true)
    end

    return prod
end

function update_inv_hess_prod(cone::WSOSInterpEpiNormOne{T}) where T
    if !cone.hess_prod_updated
        update_hess_prod(cone)
    end
    U = cone.U
    R = cone.R
    schur = cone.tempUU2
    diag_facts = cone.hess_diag_facts
    diag_blocks = cone.hess_diag_blocks
    Diz = cone.tempURU2
    schur_backup = cone.tempUU

    copyto!(schur, diag_blocks[1])

    @inbounds for r in 2:R
        r1 = r - 1
        diag_facts[r1] = posdef_fact_copy!(Symmetric(cone.hess_diags[r1], :U),
            Symmetric(diag_blocks[r], :U), false)

        z = cone.hess_edge_blocks[r1]
        LinearAlgebra.copytri!(z, 'U')
        @views Dizi = Diz[block_idxs(U, r1), :]
        ldiv!(Dizi, diag_facts[r1], z)
        mul!(schur, z', Dizi, -1, true)
    end

    cone.hess_schur_fact = posdef_fact_copy!(Symmetric(schur_backup, :U),
        Symmetric(schur, :U), false)

    cone.inv_hess_prod_updated = true
    return
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::WSOSInterpEpiNormOne,
    )
    if !cone.inv_hess_prod_updated
        update_inv_hess_prod(cone)
    end
    U = cone.U
    R = cone.R
    edge = cone.tempURU
    Diz = cone.tempURU2
    s_fact = cone.hess_schur_fact

    @inbounds for r in 1:(R - 1)
        idxs = block_idxs(U, r)
        @views copyto!(edge[:, idxs], cone.hess_edge_blocks[r])
    end
    @inbounds for r in 2:R
        idxs = block_idxs(U, r)
        @views ldiv!(prod[idxs, :], cone.hess_diag_facts[r - 1], arr[idxs, :])
    end
    # prod += u * inv(schur) * u' * arr
    ldiv!(s_fact, edge)
    @inbounds @views for j in 1:size(arr, 2)
        Dix = prod[(U + 1):end, j]
        p1j = prod[1:U, j]
        ldiv!(p1j, s_fact, arr[1:U, j])
        mul!(p1j, edge, Dix, -1, true)
        mul!(Dix, Diz, p1j, -1, true)
    end

    return prod
end

function dder3(cone::WSOSInterpEpiNormOne, dir::AbstractVector)
    @assert cone.grad_updated
    dder3 = cone.dder3
    dder3 .= 0
    R = cone.R
    U = cone.U

    # TODO refactor with other WSOS cones
    @inbounds for k in eachindex(cone.Ps)
        Pk = cone.Ps[k]
        LUk = cone.tempLU[k]
        LLk = cone.tempLL[k]
        ΛFLPk = cone.Λ11LiP[k]

        @views mul!(LUk, ΛFLPk, Diagonal(dir[1:U]))
        mul!(LLk, LUk, ΛFLPk')
        mul!(LUk, Hermitian(LLk), ΛFLPk)
        @views for u in 1:U
            dder3[u] += sum(abs2, LUk[:, u])
        end
    end
    @. dder3[1:U] *= 2 - R

    @inbounds for k in eachindex(cone.Ps)
        L = size(cone.Ps[k], 2)
        ΛLiPs11 = cone.ΛLiPs11[k]
        ΛLiPs12 = cone.ΛLiPs12[k]
        Λ11LiP = cone.Λ11LiP[k]
        LUk = cone.tempLU[k]
        ΛLiP_dir11 = cone.ΛLiP_dir11[k]
        ΛLiP_dir12 = cone.ΛLiP_dir12[k]
        ΛLiP_dir21 = cone.ΛLiP_dir21[k]
        dder3_half = cone.dder3_half[k]
        LLk = cone.tempLL[k]

        @views ΛLiP_dir22 = mul!(LUk, Λ11LiP, Diagonal(dir[1:U]))
        @views for r in 2:R
            D = Diagonal(dir[block_idxs(U, r)])
            mul!(ΛLiP_dir11[r - 1], ΛLiPs11[r - 1], Diagonal(dir[1:U]))
            mul!(ΛLiP_dir11[r - 1], ΛLiPs12[r - 1], D, true, true)
            mul!(ΛLiP_dir12[r - 1], ΛLiPs11[r - 1], D)
            mul!(ΛLiP_dir12[r - 1], ΛLiPs12[r - 1], Diagonal(dir[1:U]), true, true)
            mul!(ΛLiP_dir21[r - 1], Λ11LiP, D)
        end

        ΛLiP_dir22_Λ11LiP = ΛLiP_dir22 * Λ11LiP'
        @views for s in 1:(R - 1)
            mul!(LLk, ΛLiP_dir11[s], ΛLiPs12[s]')
            mul!(dder3_half[s][1:L, 1:U], LLk, ΛLiPs12[s])
            mul!(dder3_half[s][1:L, (U + 1):(2 * U)], LLk, ΛLiPs11[s])

            mul!(LLk, ΛLiP_dir12[s], ΛLiPs12[s]')
            mul!(dder3_half[s][1:L, 1:U], LLk, ΛLiPs11[s], true, true)
            mul!(dder3_half[s][1:L, (U + 1):(2 * U)], LLk, ΛLiPs12[s], true, true)

            mul!(LLk, ΛLiP_dir21[s], ΛLiPs12[s]')
            mul!(dder3_half[s][(L + 1):(2 * L), 1:U], LLk, ΛLiPs12[s])
            mul!(dder3_half[s][(L + 1):(2 * L), (U + 1):(2 * U)], LLk, ΛLiPs11[s])

            mul!(LLk, ΛLiP_dir22, ΛLiPs12[s]')
            mul!(dder3_half[s][(L + 1):(2 * L), 1:U], LLk, ΛLiPs11[s], true, true)
            mul!(dder3_half[s][(L + 1):(2 * L), (U + 1):(2 * U)], LLk, ΛLiPs12[s],
                true, true)

            mul!(LLk, ΛLiP_dir11[s], Λ11LiP')
            mul!(dder3_half[s][1:L, 1:U], LLk, Λ11LiP, true, true)
            mul!(LLk, ΛLiP_dir12[s], Λ11LiP')
            mul!(dder3_half[s][1:L, (U + 1):(2 * U)], LLk, Λ11LiP, true, true)
            mul!(LLk, ΛLiP_dir21[s], Λ11LiP')
            mul!(dder3_half[s][(L + 1):(2 * L), 1:U], LLk, Λ11LiP, true, true)
        end
        mul!(LLk, ΛLiP_dir22, Λ11LiP')
        mul!(LUk, LLk, Λ11LiP)

        idx = U + 1
        @views for s in 1:(R - 1)
            @. dder3_half[s][(L + 1):(2 * L), (U + 1):(2 * U)] += LUk
            for u in 1:U
                dder3[u] += sum(abs2, dder3_half[s][:, u])
                dder3[u] += sum(abs2, dder3_half[s][:, U + u])
                dder3[idx] += 2 * dot(dder3_half[s][:, U + u], dder3_half[s][:, u])
                idx += 1
            end
        end
    end

    return dder3
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    f, ax1, pl = scatter(1:4)
    f, ax2, pl = lines(1:4)
    f = Figure()
    Axis(f[1,1])
    return
end


text:
function exprify(x::NamedTuple; rename=true)
    names = Symbol[]

    function f(t::Tuple, path)
        Expr(:tuple, t...)
    end

    function f(t::NamedTuple, path)
        kvs = zip(keys(t), values(t))
        args = []
        for (k,v) in kvs
            push!(args, Expr(:(=), k, v))
        end
        Expr(:tuple, args...)
    end

    function f(x, path)
        if rename
            name = gensym()
        else
            name = x
        end
        push!(names, name)
        return name
    end

    body = fold(f, x)

    return (names, body)
end


text:
@testset "Continuous Quadratic" begin
    mock = MOIU.MockOptimizer(Model{Float64}())
    config = MOIT.TestConfig()

    @testset "QP" begin
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [4/7, 3/7, 6/7]))
        MOIT.qp1test(mock, config)
        MOIT.qp2test(mock, config)
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/4, 3/4]),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.0]))
        MOIT.qp3test(mock, config)
    end
    @testset "QCP" begin
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 7/4], MOI.FeasiblePoint))
        MOIT.qcp1test(mock, config)
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint))
        MOIT.qcp2test(mock, config)
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint))
        MOIT.qcp3test(mock, config)
    end
    @testset "SOCP" begin
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 1/2, 1/√2]))
        MOIT.socp1test(mock, config)
    end
end


text:
module Cassette

using Core: CodeInfo, SlotNumber, NewvarNode, GotoNode, SSAValue

include("context.jl")
include("pass.jl")
include("tagging.jl")
include("overdub.jl")

end # module


text:
struct Model{A,L}
	epistasis::EpistasisMu{A,L}
	x::Vector{Float64}
	λ::Vector{Float64}
	A::Int
	L::Int
	T::Int
end

#creates a model for epistasis
function Model(epistasis::EpistasisMu{A,L}, T::Int, x = [], λ = []) where {A,L}
	if length(x) == 0
		x = zeros(fields_length(epistasis))
	end
	if length(λ) == 0
		λ = zeros(T-1)
	end
	return Model(epistasis, x, λ, A, L, T)
end

#returns the total number of parameters to learn
function Model_length(A::Int, L::Int, T)
	return A^2*L*(L-1)÷2 + A*L + 1 + T -1
end

export Model, Model_length


text:
module TidalFluxPlots

using RecipesBase, Plots, TidalFluxQuantities, ADCPDataProcessing, TidalFluxCalibrations

@recipe function f(adcp::ADCPData)
    p = adcp.p
    v = adcp.v
    t = adcp.t
    temp = adcp.temp
    pitch= adcp.pitch
    roll = adcp.roll
    heading = adcp.heading

    bs = bins(adcp.dep.adcp)
    NB = length(bs)
    LB = div(NB,3)
    bsx = 1:LB:40
    bsl = bs[bsx]

    a1t = isnull(adcp.a1)
    a2t = isnull(adcp.a2)

    legend := false
    layout := @layout [pressure
                       vx
                       vy
                       vz
                       temp
                       angles
                       analog]

    @series begin
        seriestype := :path
        subplot := 1
        xticks := false
        p
    end

    @series begin
        seriestype := :heatmap
        subplot := 2
        xticks := false
        v[:,:,1]
    end

    @series begin
        seriestype := :heatmap
        subplot := 3
        xticks := false
        v[:,:,2]
    end

    @series begin
        seriestype := :heatmap
        subplot := 4
        xticks := false
        v[:,:,3]
    end

    @series begin
        seriestype := :path
        subplot := 5
        xticks := false
        temp
    end

    @series begin
        seriestype := :path
        subplot := 6
        xticks := false
        ylims := (-180,360)
        yticks := [-180; 0; 180; 360]
        [roll pitch heading]
    end

    if !a1t && !a2t
        a1 = get(adcp.a1)
        a2 = get(adcp.a2)
        @series begin
            seriestype := :path
            subplot := 7
            t, [log1p.(a1)]
        end
    else
        @series begin
            seriestype := :path
            subplot := 7
            t, zeros(length(t))
        end
    end
end

@recipe f(dep::Deployment) = load_data(dep)

@recipe function f(cal::Calibration)
    @series begin
        from_quantity(cal)
    end
    @series begin
        to_quantity(cal)
    end
end

@userplot CalScatter

@recipe function f(cs::CalScatter)
    cals = cs.args[1]
    grid := false
    leg := false

    #=
    xlabel := F.name.name
    ylabel := T.name.name
    =#

    for cal in cals
        f = TidalFluxCalibrations.interpolatecal(cal)
        t = to_quantity(cal)
        @series begin
            seriestype := :scatter
            quantity(f),quantity(t)
        end
    end
end

@recipe function f(cm::CalibrationModel{T,F};calpoints=false,ribbon=false) where {T,F}

    a,b = extrema(quantity(from_quantity(cm.c)))
    x = linspace(a,b,50)

    color --> :black
    fillcolor --> :grey70
    markerstrokecolor --> :grey70

    leg --> false
    grid --> false

#    xlabel --> F.name.name
#    ylabel --> T.name.name

    if ribbon == false
        y = predict(cm,x)
    elseif typeof(ribbon) <: Interval
        y,i = predict(cm,x,ribbon)
        ribbon := i
        linealpha := 0.0
    else
        error("$ribbon not supported")
    end

    @series begin
        x,y
    end

    if calpoints
        cal = cm.c

        f = TidalFluxCalibrations.interpolatecal(cal)
        t = to_quantity(cal)

        # Reset ribbon
        ribbon := nothing

        @series begin
            seriestype := :scatter
            quantity(f),quantity(t)
        end
    end

    # Reset ribbon and linealpha
    ribbon := nothing
    linealpha := 1.0

    @series begin
        x,y
    end
end

@recipe f(cs::CrossSectionData) = (cs.x,cs.z)

mask_quantity(q::Quantity,m::Mask) = quantity(q).*[x?1.0:NaN for x in quantity(m)]

@recipe f(q::Q) where {Q<:Quantity} = (TidalFluxQuantities.times(q),quantity(q))

@recipe f(q::Q,m::Mask) where {Q<:Quantity} = (TidalFluxQuantities.times(q),mask_quantity(q,m))

@userplot Fingerprint

@recipe function f(fp::Fingerprint)
    H,Q = fp.args

    xlabel --> "Discharge (m³/s)"
    ylabel --> "Stage (m)"

    grid --> false
    leg --> false

    framestyle --> :origin

    quantity(Q),quantity(H)
end

include("creeks.jl")

end # module


text:
using DimensionalData, Test, Unitful
using DimensionalData.LookupArrays, DimensionalData.Dimensions
using Statistics: mean

dim = Coord([(1.0,1.0,1.0), (1.0,2.0,2.0), (3.0,4.0,4.0), (1.0,3.0,4.0)], (X(), Y(), Z()))
da = DimArray(0.1:0.1:0.4, dim)
da2 = DimArray((0.1:0.1:0.4) * (1:1:3)', (dim, Ti(1u"s":1u"s":3u"s")))
Coord()

@testset "regular indexing" begin
    @test da[Coord()] === da[Coord(:)] === da
    @test da[Coord([1, 2])] == [0.1, 0.2]
    @test da[Coord(4)] == 0.4
    @test da2[Coord(4), Ti(3)] ≈ 1.2
end

@testset "coord selector indexing" begin
    @test da[Coord(:, :, :)] == [0.1, 0.2, 0.3, 0.4]
    @test da[Coord(Between(1, 5), :, At(4.0))] == [0.3, 0.4]
    @test da[Coord(:, Between(1, 3), :)] == [0.1, 0.2]
    @test da2[Ti(At(1u"s")), Coord(:, Between(1, 3), :)] ≈ [0.1, 0.2]
end

@testset "coord dimension indexing" begin
    @test da[Coord(Z(At(1.0)), Y(Between(1, 3)))] == [0.1]
end

@test index(da[Coord(:, Between(1, 3), :)], Coord) == [(1.0,1.0,1.0), (1.0,2.0,2.0)]

@test bounds(da) == (((1.0, 3.0), (1.0, 4.0), (1.0, 4.0)),)

@testset "coord named reduction" begin
    m = mean(da2; dims = Coord)
    @test size(m) == (1,3)
    @test length(dims(m, Coord)) == 1
    @test dims(m, Coord).val == DimensionalData.NoLookup(Base.OneTo(1))
    pure_mean = mean(da2.data; dims = 1)
    @test vec(pure_mean) == vec(m.data)
    @test Dimensions._tozero((1.0, 1.0)) == (0.0, 0.0)
end

@testset "coord indexing with intervals" begin
    @test da[Coord(Z(At(1.0)), Y(1..3))] == [0.1]
end


text:
Base.eltype(fg::AFG) = eltype(typeof(fg))
Base.length(fg::AFG) = length(vertices(fg))
nvertices(fg::AFG) = length(fg)
nvertices{n, T}(::Type{Simplex{n, T}}) = n
nvertices{N,T}(::Type{HyperCube{N,T}}) = 2^N
nvertices{N,T}(::Type{HyperRectangle{N,T}}) = 2^N
nvertices(s::Union{AbstractGeometry, Simplex}) = nvertices(typeof(s))

spacedim(s) = length(eltype(s))
spacedim{N,T}(::HyperCube{N,T}) = N
numtype(s) = eltype(eltype(s))
numtype{N,T}(::HyperCube{N,T}) = T

Base.push!(fl::AFG, pt) = (push!(vertices(fl), pt); fl)
Base.deleteat!(c::AbstractFlexibleGeometry, i) = (deleteat!(vertices(c), i); c)

Base.copy{FG <: AFG}(fl::FG) = FG(copy(vertices(fl)))
push(fl::AFG, pt) = push!(copy(fl), pt)

vertices(s::Simplex) = s._
standard_cube_vertices(::Type{Val{1}}) = [Vec(0), Vec(1)]
_vcat(v1,v2) = Vec(Tuple(v1)..., Tuple(v2)...)
function _combine_vcat(arr1, arr2)
    T = typeof(_vcat(first(arr1), first(arr2)))
    ret = T[]
    for v in arr1, w in arr2
        push!(ret, _vcat(v,w))
    end
    ret
end

@generated function standard_cube_vertices{N}(::Type{Val{N}})
    @assert N::Int > 0
    quote
        vert_last = standard_cube_vertices($(Val{N-1}))
        vert_1 = standard_cube_vertices(Val{1})
        _combine_vcat(vert_1, vert_last)
    end
end

@generated function vertices{N,T}(r::HyperRectangle{N,T})
    ret_type = NTuple{(2^N), Vec{N,T}}
    quote
        o = origin(r)
        v = widths(r)
        f(sv) = o + sv .* v
        tuple(map(f, standard_cube_vertices(Val{N}))...)::$ret_type
    end
end

vertices(c::HyperCube) = vertices(convert(HyperRectangle, c))
vertices(s::AbstractConvexHull) = s._

vertexmat(s::Simplex) = Mat(map(Tuple, vertices(s)))
vertexmat(s::AbstractGeometry) = Mat(map(Tuple, vertices(s)))
function vertexmat(s::AbstractFlexibleGeometry)
    tuptup = tuple(map(Tuple, vertices(s))...)
    Mat(tuptup) :: Mat{spacedim(s), nvertices(s), numtype(s)}
end
vertexmatrix(s::AbstractConvexHull) = Matrix(vertexmat(s))::Matrix{numtype(s)}

convert{S <: Simplex}(::Type{S}, fs::FlexibleSimplex) = S(tuple(vertices(fs)...))
convert{F <: AFG}(::Type{F}, s::Simplex) = F(collect(vertices(s)))
convert{FS <: FlexibleSimplex}(::Type{FS}, f::FS) = f
convert{FG <: AFG, FS <: FlexibleSimplex}(::Type{FG}, f::FS) = FG(vertices(f))
convert{R <: HyperRectangle}(::Type{R}, c::HyperCube) = R(origin(c), widths(c))
convert{F <: FlexibleConvexHull}(::Type{F}, s::Simplex) = F(collect(vertices(s)))
convert{F <: FlexibleConvexHull}(::Type{F}, c) = F(collect(vertices(c)))

function Base.isapprox(s1::AbstractConvexHull, s2::AbstractConvexHull;kw...)
    isapprox(vertexmat(s1), vertexmat(s2); kw...)
end