text:
# For each loss function we have 24 definitions!:
# CPU[4]: { 2arg, 3arg } x { sparse, dense }
# GPU[20]: { 2arg, 3arg } x { sparse, dense } x { Julia, JuliaMixed, CUDA, C32, C64 }

"""
```
loss = softloss(ypred,ygold)
softloss(ypred,ygold,ygrad)
```

Cross entropy loss to use after a softmax output (see the `soft`
activation function).  ypred has the normalized probabilities output
by the model.  ygold has normalized probabilities from the answer key.
For general arrays normalization is across the last dimension,
i.e. sum(y[:,...,:,i])==1.

For the math let us use x for the input to softmax, q for ypred, p for
ygold, J for loss.

The soft function computes: qi = exp(xi)/Σ exp(xj)
Then softloss(ypred,ygold) computes: J  = Σ -p log q

Formally softloss(ypred,ygold,ygrad) should compute: ∂J/∂q = 1-p/q
;; Derivation requires making normalization explicit:
;; J = Σi -pi log (qi/Σqj)
;; = Σi -pi log qi + pi log Σqj
;; = (Σ -pi log qi) + log Σ qj
Then softback should compute: ∂J/∂xk = Σi (∂J/∂qi)(∂qi/∂xk)
;; ∂qi/∂xk = [(i=k)(exp xi)(Σ exp xj) - (exp xi)(exp xk)] / (Σ exp xj)^2
;; = (i=k) qi - qi qk

But this is both expensive and numerically unstable.  In particular
1-p/q blows up when q is small, and trying to fix this by adding
epsilons effects the gradients and convergence in some problems.  (One
example is LeNet on MNIST scaled to [0,255], thanks to Saman Zia for
discovering this.)

It is much better to directly compute: ∂J/∂x = q-p
;; ∂J/∂xk = Σi (∂J/∂qi)(∂qi/∂xk)
;; = Σi (1-pi/qi)((i=k) qi - qi qk)
;; = Σi (qi-pi)((i=k)-qk)
;; = Σi ((i=k) qi - qi qk - (i=k) pi + pi qk)
;; = qk - qk - pk + qk = qk - pk

So have softloss(ypred,ygold,ygrad) pass back xgrad=ypred-ygold,
i.e. ∂J/∂x instead of ∂J/∂q.  Then softback just passes this back.
"""
function softloss(ypred,ygold,xgrad; mask=nothing, o...) # xgrad=(ypred-ygold)/ycols
    ycols = ccount(ypred)
    xgrad===ypred || copysync!(xgrad,ypred)
    axpy!(-1,ygold,xgrad)
    scale!(1/ycols,xgrad)
    mask!=nothing && domask(mask,xgrad)
    gpusync()
    return xgrad
end

### The two argument version is for loss calculation:

# TODO: see if we need both cpu and gpu versions or cpu only is fast enough.

function softloss(ypred::Array, ygold::Array; mask=nothing)
    @assert size(ypred)==size(ygold)
    (yrows,ycols) = size2(ygold)
    logp=zero(Float64)
    for i=1:length(ygold)
        if (mask==nothing || mask[1 + div((i-1),yrows)] != 0)
            ypredi = epsbump(ypred[i])
            logp += ygold[i]*log(ypredi)
        end
    end
    return -logp/ycols
end

epsbump(x)=(x > eps(typeof(x)) ? x : eps(typeof(x)))

softlosstemp = nothing

@gpu function softloss{T}(ypred::CudaArray{T}, ygold::CudaArray{T}; tmp=nothing, mask=nothing, o...)
    global softlosstemp
    (yrows,ycols) = size2(ygold)
    mask = (mask==nothing ? C_NULL : convert(CudaArray,mask))
    ly = (issimilar(tmp,ypred) ? tmp :
          issimilar(softlosstemp,ypred) ? softlosstemp :
          (softlosstemp = similar(ypred)))
    T <: Float32 ? ccall((:softloss32,libknet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat},Ptr{Cuchar},Ptr{Cfloat}),
                         yrows,ycols,ypred,ygold,mask,ly) :
    T <: Float64 ? ccall((:softloss64,libknet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cuchar},Ptr{Cdouble}),
                         yrows,ycols,ypred,ygold,mask,ly) : error()
    loss = CUBLAS.asum(ly)/ycols
    gpusync()
    return loss
end

@gpu softloss(ypred::CudaArray, ygold::Array; o...)=softloss(ypred, CudaArray(ygold); o...)

function softloss(ypred::Array, ygold::SparseMatrixCSC; mask=nothing, o...)
    (yrows,ycols) = size(ygold)
    logp=zero(Float64)
    col = 1             # Column i is in colptr[i]:(colptr[i+1]-1)
    for nz = 1:nnz(ygold)
        while nz > ygold.colptr[col+1]-1; col += 1; end
        mask!=nothing && mask[col]==0 && continue
        ygoldi = ygold.nzval[nz]
        row = ygold.rowval[nz]
        i = (col-1) * yrows + row
        ypredi = epsbump(ypred[i])
        logp += (ygoldi * log(ypredi))
    end
    return -logp/ycols
end

@gpu function softloss{T}(ypred::CudaArray{T}, ygold::CudaSparseMatrixCSC{T}; mask=nothing, tmp=nothing, o...)
    (yrows,ycols) = size(ygold)
    mask = (mask == nothing ? C_NULL : convert(CudaArray,mask))
    ly = (tmp == nothing ? similar(ygold.nzVal) : tmp) # TODO: get rid of alloc
    length(ly) >= nnz(ygold) || error("not enough temp space")
    T <: Float32 ? ccall((:softloss32csc,libknet),Void,(Cint,Cint,Ptr{Cfloat},Cint,Ptr{Cfloat},Ptr{Cint},Ptr{Cint},Ptr{Cuchar},Ptr{Cfloat}),
                         yrows,ycols,ypred,ygold.nnz,ygold.nzVal,ygold.rowVal,ygold.colPtr,mask,ly) :
    T <: Float64 ? ccall((:softloss64csc,libknet),Void,(Cint,Cint,Ptr{Cdouble},Cint,Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cuchar},Ptr{Cdouble}),
                         yrows,ycols,ypred,ygold.nnz,ygold.nzVal,ygold.rowVal,ygold.colPtr,mask,ly) : error()
    loss = CUBLAS.asum(nnz(ygold),ly,1)/ycols
    ly === tmp || free(ly)
    gpusync()
    return loss
end

@gpu softloss(ypred::CudaArray, ygold::SparseMatrixCSC; o...)=softloss(ypred, CudaSparseMatrixCSC(ygold); o...)

### QUADLOSS:

# Quadratic loss:
# y stores the model output.
# ygold is the desired output.
# Overwrites ygrad with the gradient of quadratic loss wrt y, i.e. y-ygold
# J = (1/2)*sum((y-ygold)^2)
# dJ/dy = y-ygold

function quadloss(ypred::BaseArray, ygold::BaseArray, ygrad::BaseArray; mask=nothing)
    length(ypred)==length(ygold)==length(ygrad) || throw(DimensionMismatch("$(map(size,(ypred,ygold,ygrad)))"))
    ycols = ccount(ypred)
    ygrad === ygold || copysync!(ygrad, ygold) # TODO: avoid copy if possible
    scale!(-1/ycols, ygrad)
    axpy!(1/ycols, ypred, ygrad)
    mask != nothing && domask(mask, ygrad)
    gpusync()
    return ygrad
end

# cpu is faster in 2-arg quadloss: 875ms->745ms
function quadloss{T}(ypred::Array{T}, ygold::Array{T}; mask=nothing)
    length(ypred)==length(ygold) || throw(DimensionMismatch("$(map(size,(ypred,ygold)))"))
    loss = 0.0
    if mask == nothing
        @inbounds for i=1:length(ypred)
            loss += (ypred[i]-ygold[i]).^2
        end
    else
        m = csize(ypred)
        @inbounds for i=1:length(ypred)
            j = 1+div(i-1,m)
            mask[j]!=0 && (loss += (ypred[i]-ygold[i]).^2)
        end
    end
    loss / (2*ccount(ypred))
end

quadloss(ypred,ygold; mask=nothing)=quadloss(convert(Array,ypred), convert(Array,ygold); mask=(mask==nothing ? mask : convert(Array,mask)))
# quadloss(ypred,ygold; o...)=quadloss(convert(Array,ypred), convert(Array,ygold); o...)

function quadloss0(ypred::BaseArray, ygold::BaseArray; mask=nothing)
    global quadlosstemp         # using this instead of alloc: 950ms->875ms
    issimilar(quadlosstemp, ypred) || (quadlosstemp = similar(ypred))
    copysync!(quadlosstemp, ygold)
    axpy!(-1, ypred, quadlosstemp)
    mask != nothing && domask(mask, quadlosstemp)
    qloss = vecnorm(quadlosstemp)^2/(2*ccount(ypred))
    gpusync()
    return qloss
end

quadlosstemp = nothing

# Should mask be bool instead of cuchar?  It turns out a bool vector
# in Julia is a char vector but only the last bit is used, other bits
# may contain garbage.  Bool type in c++ is a char vector that
# contains strictly 0 or 1.  It is safer to just use Cuchar.

@gpu function domask{T}(m::CudaVector{Cuchar},x::CudaArray{T})
    T <: Float32 ? ccall((:mask32,libknet),Void,(Cint,Cint,Ptr{Cuchar},Ptr{Cfloat}),size(x,1),size(x,2),m,x) :
    T <: Float64 ? ccall((:mask64,libknet),Void,(Cint,Cint,Ptr{Cuchar},Ptr{Cdouble}),size(x,1),size(x,2),m,x) :
    error("$T not supported")
    gpusync(); return x
end

@gpu domask{T}(m::Vector{Cuchar},x::CudaArray{T})=domask(CudaArray(m),x)

function domask(mask::Vector{Cuchar},x::Array)
    m = clength(x)
    @inbounds for i=1:length(x)
        j = 1+div(i-1,m)
        mask[j] == 0 && (x[i]=0)
    end
    return x
end

### ZERO-ONE LOSS

function zeroone(ypred::Array, ygold::Array; mask=nothing, o...)
    @assert size(ypred)==size(ygold)
    (yrows,ycols) = size2(ypred)
    cost = 0
    tmin = typemin(eltype(ypred))
    for j=1:ycols
        mask!=nothing && mask[j]==0 && continue
        (cz,cy,ymax,zmax) = (0,0,tmin,tmin)
        i1=(j-1)*yrows+1; i2=j*yrows
        for i=i1:i2
            ypred[i] > ymax && ((cy,ymax) = (i,ypred[i]))
            ygold[i] > zmax && ((cz,zmax) = (i,ygold[i]))
        end
        (cz != cy) && (cost += 1)
    end
    return cost/ycols
end

zeroone(ypred,ygold; o...)=zeroone(convert(Array,ypred),convert(Array,ygold); o...)

### DEAD CODE

# TODO: ALL THESE NEED TESTING

# Loss Layers
# TODO: get rid of l.y in documentation

# abstract Loss <: Op

# ninputs(::Loss)=1
# infersize(::Loss,dims)=(dims==nothing ? nothing : (dims,dims))
# overwrites(::Loss)=true
# back_reads_x(::Loss)=false
# back_reads_y(::Loss)=true

# # Loss has slightly different input/output behavior compared to regular layers:
# # forw only records the outgoing y.
# # back takes dy, the desired output, and returns the loss gradient wrt y
# # loss takes dy, the desired output, and returns a loss value

# for (ltype, lback, lloss, lname) in
#     ((:QuadLoss, :quadlossback, :quadlossloss, :quadloss),
#      (:SoftLoss, :softlossback, :softlossloss, :softloss),
#      (:LogpLoss, :logplossback, :logplossloss, :logploss),
#      (:XentLoss, :xentlossback, :xentlossloss, :xentloss),
#      (:PercLoss, :perclossback, :perclossloss, :percloss),
#      (:ScalLoss, :scallossback, :scallossloss, :scalloss))
#     @eval begin
#         type $ltype <: Loss; end

#         # $lname() = $ltype()

#         function forw(l::$ltype, x, y; o...)
#             size(x) == size(y) || error(map(summary,(x,y)))
#             (y===x ? y : copysync!(y,x)) # TODO: is this copy necessary?
#         end

#         function back(l::$ltype, dy, dx; y=nothing, o...)
#             size(y)==size(dy) || error(map(summary,(dy,y)))
#             dx == nothing && return
#             size(y)==size(dx) || error(map(summary,(dx,y)))
#             $lback(y,dy,dx; o...)
#         end

#         function loss(l::$ltype, dy, y; o...)
#             size(y)==size(dy) || error(map(summary,(y,dy)))
#             $lloss(y,dy; o...)
#         end

#     end
# end

# # function softlossloss(y::Array, dy::Array; o...)
# #     cost=zero(Float64)
# #     for i=1:length(dy)
# #         dy[i]>0 && (cost -= (dy[i]*log(y[i])))
# #     end
# #     return cost/ccount(dy)
# # end

# # function softlossback(y::Array, dy::Array, dx::Array; o...)
# #     nx=ccount(dx)
# #     for i=1:length(dx)
# #         dx[i] = ((y[i]-dy[i])/y[i])/nx
# #     end
# #     return dx
# # end

# # Convenience op combining soft and softloss:

# # softmax()=quote
# #     x = input()
# #     y = soft(x)
# #     z = softloss(y)
# # end

# # function quadlossloss(y::Array, dy::Array; o...)
# #     cost=zero(Float64)
# #     for i=1:length(dy)
# #         cost += (y[i]-dy[i])^2
# #     end
# #     0.5*cost/ccount(dy)
# # end

# # @gpu function quadlossloss(y::CudaArray, dy::CudaArray; tmp=nothing, o...)
# #     tmp == nothing && (tmp = similar(y)) # t:87/472
# #     copysync!(tmp, y)                        # t:29/472
# #     axpy!(-1, dy, tmp)                   # t:24/472
# #     vecnorm(tmp)^2/(2*ccount(y))         # t:330/472
# # end

# # # quadlossback(y::Array, dy::Array, dx::Array=dy; o...)=(nx=ccount(dx); for i=1:length(dx); dx[i] = (y[i]-dy[i])/nx; end; dx)
# # # @gpu quadlossback(y::CudaArray, dy::CudaArray, dx::CudaArray=dy; o...)=(dx===dy||copysync!(dx,dy); cudnnTransformTensor(1/ccount(y), y, -1/ccount(y), dx); dx)  ## cudnnTransformTensor is buggy

# # quadlossback(y, dy, dx=dy; o...)=(dx===dy||copysync!(dx,dy); scale!(-1/ccount(y), dx); axpy!(1/ccount(y), y, dx); dx)

# ### LOGPLOSS:

# # Cross entropy loss to use after the Logp layer.
# # l.y should be normalized log probabilities output by the model.
# # p has normalized probabilities from the answer key.
# # Normalization is across the last dimension, i.e. sum(p[:,...,:,i])==1
# # Overwrites p with the gradient of the loss wrt y, i.e. exp(y)-p:
# #
# # dy = sum(exp(y))   ;; normalization constant (should be 1 here)
# # q = exp(y)/dy      ;; model probabilities
# # logq = y - logz   ;; model (normalized) log prob
# # dlogz/dy = q
# #
# # J = (1/N) Σ[nc] -p[nc]*logq[nc]  ;; n=1..N: instance, c=1..C: class
# #   = (1/N) Σ[nc] -p[nc]*(y[nc]-logz[n])
# #   = (1/N) ((Σ[n] logz[n]) - (Σ[nc] p[nc]*y[nc]))
# #   = (1/N) (Σ[nc] -p[nc]*y[nc])   ;; all logz are 0
# #
# # dJ/dy[md] = (1/N) (q[md] - p[md])

# logplossloss(y::Array, dy::Array)=(nx = ccount(dy); cost = zero(Float64); for i=1:length(dy); cost -= (dy[i]*y[i]); end; cost/nx)
# logplossback(y::Array, dy::Array, dx::Array=dy)=(nx = ccount(dx); for i=1:length(dx); dx[i] = (exp(y[i])-dy[i])/nx; end; dx)
# @gpu (logplossback(y::CudaArray{Float32}, dy::CudaArray{Float32}, dx::CudaArray{Float32}=dy)=
#         (ccall((:logplossback32,libknet),Void,(Cint,Cdouble,Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}),
#                length(dy),1/ccount(dy),y,dy,dx); dx))
# @gpu (logplossback(y::CudaArray{Float64}, dy::CudaArray{Float64}, dx::CudaArray{Float64}=dy)=
#         (ccall((:logplossback64,libknet),Void,(Cint,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
#                length(dy),1/ccount(dy),y,dy,dx); dx))

# ### XENTLOSS:

# # Cross entropy loss to use after an unnormalized layer.
# # l.y is treated as unnormalized log probabilities output by the model.
# # p has normalized probabilities from the answer key.
# # Normalization is across the last dimension, i.e. sum(p[:,...,:,i])==1
# # Overwrites p with the gradient of the loss wrt y, i.e. q-p:
# #
# # z = sum(exp(y))    ;; normalization constant
# # q = exp(y)/z       ;; model probabilities
# # logq = y - logz    ;; model (normalized) log prob
# # dlogz/dy = q
# #
# # J = (1/N) Σ[nc] -p[nc]*logq[nc]  ;; n=1..N: instance, c=1..C: class
# #   = (1/N) Σ[nc] -p[nc]*(y[nc]-logz[n])
# #   = (1/N) ((Σ[n] logz[n]) - (Σ[nc] p[nc]*y[nc]))
# #
# # dJ/dy[md] = (1/N) (q[md] - p[md])

function xentloss{T}(y::Array{T}, p::Array{T})
    cost = zero(Float64)
    (nd,nx) = size2(p)
    for j=1:nx
        i1=(j-1)*nd+1; i2=j*nd
        z = sumpy = zero(Float64)
        ymax = typemin(eltype(y))
        @inbounds for i=i1:i2; y[i] > ymax && (ymax = y[i]); end
        @inbounds for i=i1:i2; yi=y[i]-ymax; z += exp(yi); sumpy += p[i]*yi; end
        cost += (log(z) - sumpy)
    end
    return cost/nx
end

@gpu xentloss{T}(y::CudaArray{T}, p::CudaArray{T})=xentloss(to_host(y),to_host(p))

function xentloss{T}(y::Array{T}, p::Array{T}, dx::Array{T})
    (nd,nx) = size2(p)
    for j=1:nx
        i1=(j-1)*nd+1; i2=j*nd
        z = zero(Float64)
        ymax = typemin(eltype(y)) # subtract ymax for numerical stability
        @inbounds for i=i1:i2; y[i] > ymax && (ymax = y[i]); end
        @inbounds for i=i1:i2; z += exp(y[i]-ymax); end
        @inbounds for i=i1:i2; yi = exp(y[i]-ymax)/z; dx[i] = (yi - p[i])/nx; end
    end
    return dx
end

@gpu (xentloss(y::CudaArray{Float32}, p::CudaArray{Float32}, dx::CudaArray{Float32})=
        ((nd,nx)=size2(p);ccall((:xentlossback32,libknet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}),nd,nx,y,p,dx);dx))
@gpu (xentloss(y::CudaArray{Float64}, p::CudaArray{Float64}, dx::CudaArray{Float64})=
        ((nd,nx)=size2(p);ccall((:xentlossback64,libknet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),nd,nx,y,p,dx);dx))

# ### PERCLOSS

# # Perceptron loss function.

# # Going forward perceptron computes y=w*x and PercLoss simply records
# # the output y.  size(w)=(nc,nd) where nc is the number of classes and
# # nd is the number of x dimensions (i.e. features).  size(x)=(nd,nx)
# # where nd is the number of features and nx is the batch size.  This
# # gives us size(y)=(nc,nx) where the highest entry in each column of y
# # indicates the predicted class.

# # Going back we get a dy matrix with size(dy)=(nc,nx) where the correct
# # answer is marked with the maximum entry in each column.
# # For a given column with input x, if cz is the correct answer and cy
# # is the predicted answer, the multiclass perceptron update rule is:

# # w[cz,:] += x;  w[cy,:] -= x

# # Note that there is no update if cz==cy.

# # The mmul updates are:
# # dw = dy*x'
# # dx = w'*dy

# # So the perceptron update will be performed if we pass a dy matrix
# # back where in each column we have all zeros if the predicted answer
# # is correct, otherwise the correct answer is marked with -1 and the
# # predicted answer is marked with a +1.  The signs might be confusing,
# # this is the gradient of the loss, i.e. going in this direction will
# # increase the loss.  We will overwrite the dy matrix.

# # This update can be seen as the gradient of a perceptron loss
# # function Sum(-y[I]+y[J]) where I are the indices for the correct
# # answers, and J are the indices for predicted answers.

# function percloss{T}(ypred::Array{T}, ygold::Array{T})
#     (nc,nx) = size2(ypred)
#     cost = zero(Float64)
#     for j=1:nx
#         (cz,cy,ymax,zmax) = (0,0,typemin(T),typemin(T))
#         i1=(j-1)*nc+1; i2=j*nc
#         for i=i1:i2
#             ypred[i] > ymax && ((cy,ymax) = (i,ypred[i]))
#             ygold[i] > zmax && ((cz,zmax) = (i,ygold[i]))
#         end
#         (cz != cy) && (cost += ypred[cy]; cost -= ypred[cz])
#     end
#     return cost/nx
# end

# function perclossback{T}(ypred::Array{T}, ygold::Array{T}, dx::Array{T}=ygold)
#     (nc,nx) = size2(ypred)
#     for j=1:nx
#         (cz,cy,ymax,zmax) = (0,0,typemin(T),typemin(T))
#         i1=(j-1)*nc+1; i2=j*nc
#         for i=i1:i2
#             ypred[i] > ymax && ((cy,ymax) = (i,ypred[i]))
#             ygold[i] > zmax && ((cz,zmax) = (i,ygold[i]))
#             dx[i] = zero(T)
#         end
#         # TODO: these should be scaled 1/nx, why isn't our gradient check complaining?
#         (cz != cy) && (dx[cz] = -1; dx[cy] = 1)
#     end
#     return dx
# end

# @gpu (perclossback(ypred::CudaArray{Float32}, ygold::CudaArray{Float32}, dx::CudaArray{Float32}=ygold)=((nd,nx)=size2(ygold);ccall((:perclossback32,libknet),Void,(Cint,Cint,Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}),nd,nx,ypred,ygold,dx);dx))
# @gpu (perclossback(ypred::CudaArray{Float64}, ygold::CudaArray{Float64}, dx::CudaArray{Float64}=ygold)=((nd,nx)=size2(ygold);ccall((:perclossback64,libknet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),nd,nx,ypred,ygold,dx);dx))

# ### SCALLOSS
# #
# # When we do structured training, gradients rather than answers come
# # back.  We should scale them using training batch size so the
# # learning rate is independent of batch size.  TODO: find a better
# # interface.  can't pass back target probabilities because network
# # output is not locally normalized.  can we pass back anything so one
# # of the existing loss functions would work?

# scallossloss(y,dy)=error("Not implemented")
# scallossback(y,dy,dx=dy)=(dx===dy||copysync!(dx,dy);scale!(1/ccount(dx), dx))

# ### DEAD CODE:

#         # TODO: can we take these out and make them apply to Loss?
#         # $lloss(y::KUdense{Array}, dy::KUdense{Array}; o...)=$lloss(convert(Array,y), convert(Array,dy); o...)
#         # $lloss(y::KUdense{Array}, dy::Array; o...)=$lloss(convert(Array,y), convert(Array,dy); o...)
#         # $lloss(y::KUdense{Array}, dy::KUdense{Array}; o...)=$lloss(convert(Array,y), convert(Array,dy); o...)
#         # $lloss(y::KUdense{Array}, dy::SparseMatrixCSC; o...)=$lloss(convert(Array,y), dy; o...)
#         # @gpu $lloss{T}(y::KUdense{CudaArray,T},dy::Array{T}; o...)=$lloss(convert(CudaArray,y), convert(CudaArray,dy); o...)
#         # @gpu $lloss{T}(y::KUdense{CudaArray,T},dy::KUdense{Array,T}; o...)=$lloss(convert(CudaArray,y), convert(CudaArray,dy); o...)
#         # @gpu $lloss{T}(y::KUdense{CudaArray,T},dy::SparseMatrixCSC{T}; o...)=$lloss(convert(CudaArray,y), convert(CudaSparseMatrixCSC,dy); o...)

#         # $lback(y::KUdense{Array}, dy::KUdense{Array}, dx::KUdense{Array};o...)=($lback(convert(Array,y), convert(Array,dy), convert(Array, dx);o...); dx)
#         # $lback(y::KUdense{Array}, dy::SparseMatrixCSC, dx::KUdense{Array};o...)=($lback(convert(Array,y), dy, convert(Array, dx);o...); dx)
#         # @gpu $lback(y::KUdense{CudaArray}, dy::KUdense{CudaArray}, dx::KUdense{CudaArray};o...)=($lback(convert(CudaArray,y), convert(CudaArray,dy), convert(CudaArray, dx);o...); dx)
#         # @gpu $lback(y::KUdense{CudaArray}, dy::CudaSparseMatrixCSC, dx::KUdense{CudaArray};o...)=($lback(convert(CudaArray,y), dy, convert(CudaArray, dx);o...); dx)

#         # $lback(y::KUdense, dy::KUdense, dx::KUdense=dy)=($lback(y.arr, dy.arr, dx.arr); dx)
#         # $lloss(y,dy)=$lloss(convert(Array,y), convert(Array,dy))  # TODO: handle sparse arrays, implement gpu

# # params(::Loss)=Any[]
# # overwrites(::Loss)=true

# function softloss(ypred::Array, ygold::Array, ygrad::Array; mask=nothing, o...)
#     @assert size(ypred)==size(ygold)==size(ygrad)
#     (yrows,ycols) = size2(ygrad)
#     for i=1:length(ygrad)
#         if (mask != nothing && mask[1 + div((i-1),yrows)] == 0)
#             ygrad[i] = 0
#         else
#             ypredi = epsbump(ypred[i])
#             ygrad[i] = ((ypredi-ygold[i])/ypredi)/ycols
#         end
#     end
# end

# epsbump(x)=(x > eps(typeof(x)) ? x : eps(typeof(x)))

# @gpu function softloss{T}(ypred::CudaArray{T}, ygold::CudaArray{T}, ygrad::CudaArray{T}; mask=C_NULL, o...)
#     (yrows,ycols) = size2(ygrad)
#     mask != C_NULL && (mask = convert(CudaArray, mask))
#     T <: Float32 ? ccall((:softlossback32,libknet),Void,(Cint,Cint,Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cuchar}, Ptr{Cfloat}),
#                          yrows,ycols,ypred,ygold,mask,ygrad) :
#     T <: Float64 ? ccall((:softlossback64,libknet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cuchar},Ptr{Cdouble}),
#                          yrows,ycols,ypred,ygold,mask,ygrad) : error()
#     gpusync()
#     return ygrad
# end

# @gpu softloss(ypred::CudaArray, ygold::Array, ygrad::CudaArray; o...)=softloss(ypred, CudaArray(ygold), ygrad; o...)

# function softloss(ypred::Array, ygold::SparseMatrixCSC, ygrad::Array; mask=nothing, o...)
#     (yrows, ycols) = size(ygold)
#     col = 1             # Column i is in colptr[i]:(colptr[i+1]-1)
#     for nz = 1:nnz(ygold)
#         while nz > ygold.colptr[col+1]-1; col += 1; end
#         if mask != nothing && mask[col] == 0
#             ygrad[i] = 0
#         else
#             ygoldi = ygold.nzval[nz]
#             row = ygold.rowval[nz]
#             i = (col-1) * yrows + row
#             ypredi = epsbump(ypred[i])
#             ygrad[i] = ((ypredi-ygoldi)/ypredi)/ycols
#         end
#     end
#     return ygrad
# end

# @gpu function softloss{T}(ypred::CudaArray{T}, ygold::CudaSparseMatrixCSC{T}, ygrad::CudaArray{T}; mask=C_NULL, o...)
#     (yrows,ycols) = size2(ygrad)
#     mask != C_NULL && (mask = convert(CudaArray, mask))
#     T <: Float32 ? ccall((:softlossback32csc,libknet),Void,(Cint,Cint,Ptr{Cfloat}, Cint,Ptr{Cfloat}, Ptr{Cint},Ptr{Cint},Ptr{Cuchar},Ptr{Cfloat}),
#                          yrows,ycols,ypred,ygold.nnz,ygold.nzVal,ygold.rowVal,ygold.colPtr,mask,ygrad) :
#     T <: Float64 ? ccall((:softlossback64csc,libknet),Void,(Cint,Cint,Ptr{Cdouble},Cint,Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cuchar},Ptr{Cdouble}),
#                          yrows,ycols,ypred,ygold.nnz,ygold.nzVal,ygold.rowVal,ygold.colPtr,mask,ygrad) : error()
#     gpusync()
#     return ygrad
# end
# @gpu softloss(ypred::CudaArray, ygold::SparseMatrixCSC, ygrad::CudaArray; o...)=softloss(ypred, CudaSparseMatrixCSC(ygold), ygrad; o...)


text:
function Quote(ex::N) where {N <: Number}
    if ex isa Int || ex isa Float64 || ex isa Float32
        return ex
    else
        return convert(Expr, ex)
    end
end

# QuoteFn
function Quote(name::A,ex::Expr,vars=free_symbols(ex)) where {A <: AbstractString}
    return Expr(:function, Expr(:call, Symbol(name), map(Symbol,vars)...), Quote(ex))
end

# QuoteFnCSE
function Quote(name::A,cse_ex::Tuple{<:Array, Array{Sym,1}},vars::Array{Sym,1}) where {A <: AbstractString}
    quot = Expr(:block)
    quot.head = :function
    push!(quot.args, :($(Symbol(name))($(map(Symbol,vars)...) )))
    fnbody = Quote(cse_ex)
    push!(quot.args,fnbody)
    return quot
end

# QuoteCSE
function Quote(cse_ex::Tuple{T, Array{Sym,1}}) where T <: AbstractArray
    cse,ex = cse_ex
    ex=ex[1]
    fnbody=Expr(:block)
    for (k,v) in cse
        push!(fnbody.args, :($(Quote(k)) = $(Quote(v))) )
    end
    push!(fnbody.args,Quote(ex))
    return fnbody
end

#QuoteArrCSE
function Quote(cse_ex::Tuple{T,Array{U,1}}) where {T <: AbstractArray, U <: Array{Sym}}
    # assume we did simultaneous CSE on the array
    # cse is common list of substitutions
    # ex is array of cse'd expressions
    cse,ex = cse_ex
    fnbody=Expr(:block)
    for (k,v) in cse
        push!(fnbody.args, :($(Quote(k)) = $(Quote(v))) )
    end
    if length(ex) == 1
        ex = ex[1]
    end
    #push!(fnbody.args,:([$(Quote.(ex)...)]))
    push!(fnbody.args,Quote(ex).args[1])
    return fnbody
end

# QuoteArr
function Quote(ex::Array{Sym})
    fnbody=Expr(:block)
    sh = size(ex)
    push!(fnbody.args,:(reshape([$(Quote.(ex)...)],$sh)))
    # possibly improve this with alg at https://github.com/symengine/SymEngine.jl/pull/114
    return fnbody
end

# QuoteFnArr
function Quote(name::A,ex::Array{Sym}, vars=free_symbols(ex)) where {A <: AbstractString}
    Expr(:function, Expr(:call, Symbol(name), map(Symbol,vars)...), Quote(ex).args[1])
end

# QuoteFnArrCSE
function Quote(name::A, cse_ex::Tuple{<:AbstractArray,Array{<:Array{Sym},1}}, vars::Array{Sym,1}) where {A <: AbstractString}
    quot = Expr(:block)
    quot.head = :function
    push!(quot.args, :($(Symbol(name))($(map(Symbol,vars)...) )))
    fnbody = Quote(cse_ex)
    push!(quot.args,fnbody)
    return quot
end

ext: jl


text:
using HCubature, Cuba, SimpleLinearPiecewiseApproximation, TransferFunctions, Printf

using SpecialFunctions, Interpolations
include("integrals.jl")

include("GKGT.jl")
include("convkernel.jl")
include("convkernel_cuba.jl")


text:
#!jl # # EFTfitter.jl - Tutorial
#jl # EFTfitter.jl - Tutorial

#jl #~This tutorial introduces the basic functionalities of EFTfitter.jl using a generic example.
#jl #~More functionalities of EFTfitter.jl, like handling nuisance correlations
#jl #~or ranking measurements and uncertainties, are shown in the
#jl #~[advanced tutorial](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/advanced_tutorial/).
#md # Table of contents:
#md # ```@contents
#md # Pages = ["tutorial.md"]
#md # Depth = 3
#md # ```
#nb #~First, we need to setup EFTfitter, BAT and some other Julia packages:
#nb using EFTfitter
#nb using BAT           # for sampling
#nb using IntervalSets  # for specifying the prior
#nb using Distributions # for specifying the prior
#nb using Plots         # for plotting

#!jl include("lit_tutorial_inputs.jl")
#!jl # ## File "runTutorial.jl"
#~Here, we create the `EFTfitterModel` from our inputs and run the actual analysis.

#!nb #~First, we need to setup EFTfitter, BAT and some other Julia packages:
#!nb using EFTfitter
#!nb using BAT            # for sampling
#!nb using IntervalSets   # for specifying the prior
#!nb using Distributions  # for specifying the prior
#!nb using Plots # for plotting
#jl
#jl #~We include the definitions of observables, measurements,
#jl #~uncertainties and correlations from our `tutorial_inputs.jl` file:
#jl include("tutorial_inputs.jl")

#~We can then build the `EFTfitterModel` which combines all our inputs into
#~one object that is then used to perform the analysis on.
model = EFTfitterModel(parameters, measurements, correlations)

#~To sample the posterior distribution, we specify that our `EFTfitterModel`
#~should be used and then setup BAT.jl to sample the EFTfitter likelihood.
posterior = PosteriorDensity(model)

algorithm = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)
samples = bat_sample(posterior, algorithm).result;
#~For further information on settings & algorithms when sampling with BAT.jl
#~see the BAT.jl [tutorial](https://bat.github.io/BAT.jl/dev/tutorial/#Parameter-Space-Exploration-via-MCMC)
#~and [documentation](https://bat.github.io/BAT.jl/dev/stable_api/#BAT.bat_sample).

#~We can then inspect the results of the sampling using BAT.jl's `SampledDensity`,
#~giving a summary of the sampling and the results of the model parameters.
sd = SampledDensity(posterior, samples)
display(sd)
#md ```
#md ```
#md BAT.jl - SampledDensity
#md ──────────────────────────────
#md
#md Sampling:
#md ─────────────────────────
#md total number of samples:      103343
#md effective number of samples: (C1 = 19897.99484950266, C2 = 18106.884734831812)
#md
#md
#md Parameter estimates:
#md ─────────────────────────
#md number of free parameters: 2
#md
#md Table with 5 columns and 2 rows:
#md      parameter  mean       std        global_mode  marginal_mode
#md    ┌────────────────────────────────────────────────────────────
#md  1 │ C1         0.864514   0.349673   0.821192     0.81
#md  2 │ C2         0.0129843  0.0210368  0.0155112    0.0165
#md
#md
#md Covariance matrix:
#md ─────────────────────────
#md 2×2 Named Array{Float64,2}
#md cov ╲  │          C1           C2
#md ───────┼─────────────────────────
#md C1     │    0.122271   -0.0070394
#md C2     │  -0.0070394  0.000442548

#~Information about the smallest 1d intervals containing p% proability can be
#~obtained using the `get_smallest_interval_edges` function:
intervals_1d_C1 = get_smallest_interval_edges(samples, :C1, 0.9, bins=200, atol=0.1)
println("lower interval edges: $(intervals_1d_C1.lower)")
println("upper interval edges: $(intervals_1d_C1.upper)")

#~The keyword `atol` controls the absolute tolerance for which intervals are joined
#~together when they are seperated less than this value. This is particularly useful
#~when a large number of bins is used.

#~Of course, plotting the resulting posterior distributions is also simple
#~using Plots.jl and the BAT.jl plotting recipes:
p = plot(samples)
#!nb savefig(p, "plot.pdf")
#md # ![example plot](plots/plot.png)

#~For information on how to customize plots of the samples, please see the BAT.jl
#~[plotting documentation](https://bat.github.io/BAT.jl/dev/plotting/) and
#~[examples](https://github.com/bat/BAT.jl/blob/master/examples/dev-internal/plotting_examples.jl).

p = plot(samples, 0.9)
#!nb savefig(p, "plot_1d.pdf")
#md # ![example plot](plots/interval_plot_1.png)

#~For customizing the plots of the 1D intervals, also see the EFTfitter
#~[plotting documentation](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/plotting/)
#~and [tutorial](https://github.com/tudo-physik-e4/EFTfitter.jl/blob/main/examples/tutorial/plotting.jl).


text:
function return_indices(ary::Vector)
    u = unique(ary)
    idx = enumerate(u)
    ticks = similar(u)
    labels = Array{String}(size(u))

    indices = Vector{Int64}(length(ary))
    for x in idx
        indices[ary.==x[2]] = x[1]
        ticks[x[1]] = x[1]
        labels[x[1]] = "$(x[2])"
    end

    rng = (minimum(u),maximum(u))

    return indices, ticks, labels, rng
end

## Code below is not fully implemented

function getsample(n::Int64,p::Int64,smooth::Symbol)::Array{Float64}
    x = rand(n,p)

    for i = 1:p
        x(:,i) = rank(x(:,i))
    end
    if smooth == :on
        x = x - rand(size(x))
    else
        x = x - 0.5
    end
    return x/n
end

function score(x::Int64,crit::Symbol)
    if size(x,1) < 2
        # score is meaningless with just one point
        return 0.
    end
    if crit == :correlation
        error("The '$crit' function has not been implemented yet")
    elseif crit == :maximin
        error("The '$crit' function has not been implemented yet")
    else
        error("argument not recognized")
    end
end

function lhs_rank(x::Vector{Float64})::Vector{Float64}
    z = sort(x)
    r = similar(x)
    for i = 1:length(z)
        r[z[i]] = i
    end
    return r
end

function load_network(nn_prefix::String, epoch::Int,sq_fname::String)
    # the default one doesn't work for some reason, so I'll do it by hand
    arch1, arg_params1, aux_params1 = mx.load_checkpoint(string(nn_prefix,"_net1"),epoch)
    arch2, arg_params2, aux_params2 = mx.load_checkpoint(string(nn_prefix,"_net2"),epoch)

    net1 = mx.FeedForward(arch1)
    net1.arg_params = arg_params1
    net1.aux_params = aux_params1

    net2 = mx.FeedForward(arch2)
    net2.arg_params = arg_params2
    net2.aux_params = aux_params2

    sq_data = JLD.load(sq_fname)
    input_sch = sq_data["input_sch"]
    output_sch = sq_data["output_sch"]
    range = sq_data["range"]

    SQ = SQ_model(input_sch,output_sch,net1,net2,range)
    return SQ
end
function make_label_from_keys(d::Dict)
    z = ""
    for x in collect(keys(d))
        z = string(z,x)
    end
    return z
end
function searchdir(path,key1)
    # list of files matching key1
    filt_list = filter(x->contains(x,key1),readdir(path))

    return filt_list
end
function searchdir(path,key1,key2)
    # list of files matching key1
    filt_list = searchdir(path,key1)

    # subset of that list that contains key2
    filt_list2 = filter(x->contains(x,key2),filt_list)

    return filt_list2
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
function _laplace_single_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    eidx    ::Int
) where T
    elem = elements[eidx]
    for ξidx in 1:length(Ξ)
        dst[ξidx] = laplacepot_single(Ξ[ξidx], elem)
    end
    nothing
end

function _laplace_double_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    eidx    ::Int
) where T
    elem = elements[eidx]
    for ξidx in 1:length(Ξ)
        dst[ξidx] = laplacepot_double(Ξ[ξidx], elem)
    end
    nothing
end

@testset "[sl] ξ is vertex" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([elem.v1, elem.v2, elem.v3])
        dst = CuArray{T}(undef, 3)

        @cuda _laplace_single_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(SingleLayer, elem.v1, elem, zero(T))
        @test res[2] ≈ Rjasanow.laplacepot(SingleLayer, elem.v2, elem, zero(T))
        @test res[3] ≈ Rjasanow.laplacepot(SingleLayer, elem.v3, elem, zero(T))
    end
end

@testset "[sl] ξ in plane" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[0, -1, 0], T[0, -1, -1]])
        dst = CuArray{T}(undef, 2)

        @cuda _laplace_single_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(SingleLayer, T[0, -1, 0], elem, zero(T))
        @test res[2] ≈ Rjasanow.laplacepot(SingleLayer, T[0, -1, -1], elem, zero(T))
    end
end

@testset "[sl] ξ on normal through vertex" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[1, 0, 0], T[1, 1, 0], T[1, 0, 1]])
        dst = CuArray{T}(undef, 3)

        @cuda _laplace_single_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(SingleLayer, elem.v1, elem, one(T))
        @test res[2] ≈ Rjasanow.laplacepot(SingleLayer, elem.v2, elem, one(T))
        @test res[3] ≈ Rjasanow.laplacepot(SingleLayer, elem.v3, elem, one(T))
    end
end

@testset "[sl] ξ somewhere in space" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[1, -1, 0], T[1, -1, -1]])
        dst = CuArray{T}(undef, 2)

        @cuda _laplace_single_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(SingleLayer, T[0, -1, 0], elem, one(T))
        @test res[2] ≈ Rjasanow.laplacepot(SingleLayer, T[0, -1, -1], elem, one(T))
    end
end

@testset "[dl] ξ is vertex" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([elem.v1, elem.v2, elem.v3])
        dst = CuArray{T}(undef, 3)

        @cuda _laplace_double_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v1, elem, zero(T))
        @test res[2] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v2, elem, zero(T))
        @test res[3] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v3, elem, zero(T))
    end
end

@testset "[dl] ξ in plane" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[0, -1, 0], T[0, -1, -1]])
        dst = CuArray{T}(undef, 2)

        @cuda _laplace_double_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(DoubleLayer, T[0, -1, 0], elem, zero(T))
        @test res[2] ≈ Rjasanow.laplacepot(DoubleLayer, T[0, -1, -1], elem, zero(T))
    end
end

@testset "[dl] ξ on normal through vertex" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[1, 0, 0], T[1, 1, 0], T[1, 0, 1]])
        dst = CuArray{T}(undef, 3)

        @cuda _laplace_double_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v1, elem, one(T))
        @test res[2] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v2, elem, one(T))
        @test res[3] ≈ Rjasanow.laplacepot(DoubleLayer, elem.v3, elem, one(T))
    end
end

@testset "[dl] ξ somewhere in space" begin
    for T in [Float32, Float64]
        elem = Triangle(T[0, 0, 0], T[0, 1, 0], T[0, 0, 1])
        elements = TriangleVector([elem])
        Ξ = PositionVector([T[1, -1, 0], T[1, -1, -1]])
        dst = CuArray{T}(undef, 2)

        @cuda _laplace_double_kernel!(dst, Ξ, elements, 1)
        res = Array(dst)
        @test res[1] ≈ Rjasanow.laplacepot(DoubleLayer, T[0, -1, 0], elem, one(T))
        @test res[2] ≈ Rjasanow.laplacepot(DoubleLayer, T[0, -1, -1], elem, one(T))
    end
end


text:
##########################
# Earth Orientation Data #
##########################

# Test internal data type
let
    @test typeof(EOP.data) == Dict{Int,Tuple{Float64,Float64,Float64}}
end

# Test constructors
let
    finals_2000 = EarthOrientationData("FINALS_2000")
    eop_c04     = EarthOrientationData("C04_14")
end

# Test Data valutes
let
    @test isapprox(UT1_UTC(58483), -0.0351948, atol=1e-5)
    xp, yp = POLE_LOCATOR(58483)
    @test isapprox(xp, 0.088504*AS2RAD, atol=1e-5)
    @test isapprox(yp, 0.270753*AS2RAD, atol=1e-5)
    @test isapprox(XP(58483), 0.088504*AS2RAD, atol=1e-5)
    @test isapprox(YP(58483), 0.270753*AS2RAD, atol=1e-5)
end

# Test Set values
let
    set_eop(58747, -0.2, 0.225, 0.3)

    @test UT1_UTC(58747)      == -0.2
    @test POLE_LOCATOR(58747) == (0.225*AS2RAD, 0.3*AS2RAD)
    @test XP(58747)           == 0.225*AS2RAD
    @test YP(58747)           == 0.3*AS2RAD

    set_eop(58747, -0.1897929, 0.230292, 0.332704)
    @test UT1_UTC(58747)      == -0.1897929
    @test POLE_LOCATOR(58747) == (0.230292*AS2RAD, 0.332704*AS2RAD)
    @test XP(58747)           == 0.230292*AS2RAD
    @test YP(58747)           == 0.332704*AS2RAD
end

# Test interpolation values
let
    ut1_utc = (UT1_UTC(58748) + UT1_UTC(58747))/2.0
    @test UT1_UTC(58747.5, interp=true) == ut1_utc

    x1, y1 = POLE_LOCATOR(58747)
    x2, y2 = POLE_LOCATOR(58748)
    pole_locator = ((x2+x1)/2.0, (y2+y1)/2.0)
    @test POLE_LOCATOR(58747.5, interp=true)[1] == pole_locator[1]
    @test POLE_LOCATOR(58747.5, interp=true)[2] == pole_locator[2]

    xp = (XP(58748) + XP(58747))/2.0
    @test XP(58747.5, interp=true) == xp

    yp = (YP(58748) + YP(58747))/2.0
    @test YP(58747.5, interp=true) == yp
end

# Test loading new products into project global vairable
let
    load_eop("C04_14")
    @test UT1_UTC(37665) == 0.0326338

    # Return to original state
    load_eop("FINALS_2000")
end

# Test Update of EOP
# let
#     update_eop("C04_14")
#     update_eop("C04_80")
#     update_eop("FINALS_2000")
# end

#################
# Gravity Model #
#################

let
    @test typeof(GRAVITY_MODEL.data) == Array{Float64, 2}

    @test GRAVITY_MODEL.name == "EGM2008"
    @test GRAV_COEF(0, 0)   == 1.0
    @test GRAV_COEF(90, 90) == 0.733188520723327e-09
end

let
    gmdl = GravModel(:EGM2008_90)
end

let
    load_gravity_model(:GGM01S)
end


text:
# ================================================
# Gearbox model
# See https://easychair.org/publications/paper/cwl
# ================================================
using Reachability, HybridSystems, MathematicalSystems, MathematicalPredicates,
      SparseArrays
using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra

"""
    gearbox()

Construct the gearbox model.
"""
function gearbox()
    # variables
    vx = 1  # x velocity
    vy = 2  # y velocity
    px = 3  # x position
    py = 4  # y position
    𝐼 = 5  # accumulated impulse (italic version to avoid a clash with 'I')
    t = 6   # time

    # number of variables
    n = 5 + 1

    # constants
    Fs = 70.0  # shifting force on the sleeve [N]
    Tf = 1.0  # resisting moments on the second gear [N⋅m]
    # TODO where do the following values come from?
    ms = 3.2  # mass of the sleeve
    Rs = 0.08  # radius of the sleeve
    Jg₂ = 0.7  # inertia of the second gear
    Δp = -0.003  # horizontal (px) distance
    θ = 0.628318530717959  # included angle of the second gear [°]
    mg₂ = 18.1  # mass of the second gear
    ζ = 0.9  # coefficient of restitution

    # additional bloating to get a non-flat guard intersection
    guard_bloating = sqrt(eps(Float64))

    # discrete structure (graph)
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 1, 1)
    add_transition!(automaton, 1, 1, 2)
    add_transition!(automaton, 1, 2, 3)

    # mode 1 ("free")
    A = zeros(n, n)
    b = zeros(n)
    A[px, vx] = 1.
    A[py, vy] = 1.
    b[vx] = Fs / ms
    b[vy] = - (Rs * Tf) / Jg₂
    b[t] = 1.
    invariant = HalfSpace(sparsevec([px], [1.], n), Δp)
    # TODO The SpaceEx model adds more constraints, possibly to help with the
    # guard intersection:
    invariant = HPolyhedron([invariant,
        HalfSpace(sparsevec([px, py], [tan(θ), 1.], n), 0.),    # py <= -px * tan(θ)
        HalfSpace(sparsevec([px, py], [tan(θ), -1.], n), 0.)])  # py >= px * tan(θ)
    m_1 = CACS(A, b, invariant)

    # mode 2 ("meshed")
    A = zeros(n, n)
    b = zeros(n)
    b[t] = 1.
    invariant = Universe(n)
    m_2 = CACS(A, b, invariant)

    # modes
    modes = [m_1, m_2]

    # common assignment matrix (requires individual modifications)
    A_template = zeros(n, n)
    for i in 3:6
        A_template[i, i] = 1.
    end
    denominator = ms * cos(θ)^2 + mg₂ * sin(θ)^2
    A_template[vx, vx] = (ms * cos(θ)^2 - mg₂ * ζ * sin(θ)^2) / denominator
    A_template[vx, vy] = (-(ζ + 1.) * mg₂ * sin(θ) * cos(θ)) / denominator
    A_template[vy, vx] = (-(ζ + 1.) * ms * sin(θ) * cos(θ)) / denominator
    A_template[vy, vy] = (mg₂ * sin(θ)^2 - ms * ζ * cos(θ)^2) / denominator
    A_template[𝐼, vx] = ((ζ + 1.) * ms * mg₂ * sin(θ)) / denominator
    A_template[𝐼, vy] = ((ζ + 1.) * ms * mg₂ * cos(θ)) / denominator

    # transition l1 -> l1
    # TODO what happened to the term '2nb' and the whole second constraint in the paper?
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), -1.], n), 0.),  # py >= -px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), -cos(θ)], n), 0.)  # vx * sin(θ) + vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    t1 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l1
    # TODO same remark as with the other guard
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), 1.], n), 0.),  # py <= px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), cos(θ)], n), 0.)  # vx * sin(θ) - vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    A[vx, vy] *= -1.
    A[vy, vx] *= -1.
    A[𝐼, vy] *= -1.
    t2 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l2
    guard = HalfSpace(sparsevec([px], [-1.], n), -Δp + guard_bloating)  # px >= Δp - ε
    A = copy(A_template)
    A[vx, vx] = 0.
    A[vx, vy] = 0.
    A[vy, vx] = 0.
    A[vy, vy] = 0.
    A[𝐼, vx] = A[𝐼, vy] = ms
    t3 = ConstrainedLinearMap(A, guard)
    # TODO The SpaceEx model uses four transitions for this one, possibly to
    # help with the signs; note that the assignments are different:
    # * guard = px >= Δp & vx >= 0 & vy >= 0
    #   assignment = I:=I+ms*vx+ms*vy & vx:=0 & vy:=0
    # * guard = px >= Δp & vx >= 0 & vy <= 0
    #   assignment = I:=I+ms*vx-ms*vy & vx:=0 & vy:=0
    # * guard = px >= Δp & vx <= 0 & vy >= 0
    #   assignment = I:=I-ms*vx+ms*vy & vx:=0 & vy:=0
    # * guard = px >= Δp & vx <= 0 & vy <= 0
    #   assignment = I:=I-ms*vx-ms*vy & vx:=0 & vy:=0

    # transition annotations
    resetmaps = [t1, t2, t3]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Hyperrectangle([0., 0., -0.0167, 0.003,  0., 0.],
                        [0., 0.,  0.0001, 0.0001, 0., 0.])
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    # safety property
    cond_free = is_contained_in(HalfSpace([0., 0., 0., 0., 0., 1.], 0.2))    # t <= 0.2
    cond_global = is_contained_in(HalfSpace([0., 0., 0., 0., 1., 0.], 20.))  # I <= 20
    property_free = Conjunction([cond_free, cond_global])
    property_meshed = cond_global
    property = Dict(1 => property_free, 2 => property_meshed)

    # default options
    options = Options(:T=>0.5, :property=>property)

    return (system, options)
end

function run_gearbox(system, options)
    opC = BFFPSV18(:δ => 0.001)
    opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => false)
    options[:verbosity] = "info"
    options[:mode] = "check"
    options[:plot_vars] = [3, 4]
    options[:project_reachset] = true

    solve(system, options, opC, opD)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function incident_infections(p1::InfParams, p2::InfParams, S::Int, I1::Vector{Int}, I2::Vector{Int}, R::Int)
    iszero(S) && return 0,0
    infSum1 = 0
    infSum2 = 0
    N = S + sum(I1) + sum(I2) + R
    susceptible_prop = S/N

    for (i, inf) in enumerate(I1)
        d = p1.infectiousness[i]
        k,θ = Distributions.params(d)
        d′ = Gamma(k*susceptible_prop, θ) # conversion from R₀ to Rₜ
        infSum1 += rand(RVsum(d′, inf))
    end

    for (i, inf) in enumerate(I2)
        d = p2.infectiousness[i]
        k,θ = Distributions.params(d)
        d′ = Gamma(k*susceptible_prop, θ) # conversion from R₀ to Rₜ
        infSum2 += rand(RVsum(d′, inf))
    end

    infSum1 = floor(Int, infSum1)
    infSum2 = floor(Int, infSum2)
    inf_sum = infSum1 + infSum2
    overage = inf_sum - S

    if overage > 0
        half_over = div(overage, 2)
        remove1 = min(infSum1, half_over)
        infSum1 -= remove1
        overage -= remove1
        infSum2 -= overage
    end

    return infSum1, infSum2
end

function symptomatic_isolation(params::InfParams, I::Vector{Int})
    isolating = zero(I)
    for (i, inf) in enumerate(I)
        symptomatic_prob = params.symptom_probs[i]
        isolation_prob = symptomatic_prob*(1-params.asymptomatic_prob)*params.symptomatic_isolation_prob
        isolating[i] = rand(Binomial(inf,isolation_prob))
    end
    return isolating
end

function positive_tests(params::InfParams, I::Vector{Int}, T::Matrix{Int}, a::CovidAction)
    pos_tests = zeros(Int, length(params.pos_test_probs))

    for (i, inf) in enumerate(I)
        num_already_tested = sum(@view T[:,i])
        num_tested = floor(Int,(inf-num_already_tested)*a.testing_prop)
        pos_tests[i] = rand(Binomial(num_tested,params.pos_test_probs[i]))
    end
    return pos_tests
end

# Only record number that have taken the test, the number that return postive is
# Binomial dist, such that s' is stochastic on s.
function update_isolations(params::InfParams, I, R, T, a::CovidAction)

    sympt = symptomatic_isolation(params, I) # Number of people isolating due to symptoms
    pos_tests = positive_tests(params, I, T, a)

    sympt_prop = sympt ./ I # Symptomatic Isolation Proportion
    replace!(sympt_prop, NaN=>0.0)

    R += sum(sympt)
    I .-= sympt

    T[end,:] .= pos_tests

    @. T = floor(Int, (1 - sympt_prop)' * T)

    R += sum(@view T[1,:])
    @views I .-= T[1,:]

    @assert all(≥(0), I)

    CovidPOMDPs.shift_test!(T)

    return I, R, T, pos_tests
end

function CovidPOMDPs.sim_step(pomdp::DoubleCovidPOMDP, state::DoubleCovidState, a::CovidAction)
    (;S, I1, I2, R, T1, T2, params, prev_action) = state
    I1 = copy(I1)
    I2 = copy(I2)
    T1 = copy(T1)
    T2 = copy(T2)

    # Update symptomatic and testing-based isolations
    I1, R, T1, pos_tests1 = update_isolations(pomdp.params, I1, R, T1, a)
    I2, R, T2, pos_tests2 = update_isolations(params, I2, R, T2, a)

    # Incident Infections
    R += I1[end]
    R += I2[end]
    I1 = circshift(I1, 1) # shift_inf!(I1)
    I2 = circshift(I2, 1) # shift_inf!(I2)

    inf1, inf2 = incident_infections(pomdp.params, params, S, I1, I2, R)
    I1[1] = inf1
    I2[1] = inf2

    S -= (inf1 + inf2)
    sp = DoubleCovidState(S, I1, I2, R, T1, T2, params, a)

    return sp, inf1+inf2, pos_tests1 .+ pos_tests2
end


text:
import ..update!, ..select_action, ..reset!

using Random
using StaticArrays

using Ramnet
using Ramnet.Models.AltDiscriminators: FunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

struct DiscreteActionPolicyLoss <: Loss.AbstractLoss end

struct DiscreteActionPolicy{OS,AS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    opt::FunctionalOptimizer{DiscreteActionPolicyLoss}
    f::FunctionalDiscriminator{AS,T,E}
    rng::MersenneTwister
end

function DiscreteActionPolicy{AS}(::Type{O}, n::Int, encoder::E; start_learning_rate::Float64 = 0.1, end_learning_rate::Float64 = 1e5, learning_rate_decay::Int = 1000, epochs::Int = 1, partitioner::Symbol = :uniform_random, seed::Union{Nothing,UInt} = nothing) where {AS,OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    DiscreteActionPolicy{OS,AS,T,O,E}(
        FunctionalOptimizer(DiscreteActionPolicyLoss(); start_learning_rate, end_learning_rate, learning_rate_decay, epochs),
        FunctionalDiscriminator{AS}(OS, n, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
end

function Loss.grad(::DiscreteActionPolicyLoss, p::DiscreteActionPolicy{OS,AS,T,O,E}, observation::O, action::Int, G) where {OS,AS,T,O,E}
    probs = probabilities(p, observation)
    probs[action+1] -= 1
    return probs * G
end

function probabilities(p::DiscreteActionPolicy{OS,AS,T,O,E}, observation::O) where {OS,AS,T,O,E}
    logits = predict(p.f, observation)
    return exp.(logits) / sum(exp.(logits))
end

function probability(p::DiscreteActionPolicy{OS,AS,T,O,E}, action::Int, observation::O) where {OS,AS,T,O,E}
    return probabilities(p, observation)[action+1]
end

function mode(p::DiscreteActionPolicy{OS,AS,T,O,E}, observation::O) where {OS,AS,T,O,E}
    return argmax(predict(p.f, observation)) - 1
end

function update!(p::DiscreteActionPolicy{OS,AS,T,O,E}, observation::O, action::Int, G) where {OS,AS,T,O,E}
    gradient = Optimizers.grad(p.opt, p, observation, action, G)

    train!(p.f, observation, learning_rate(p.opt) * gradient)

    return nothing
end

function select_action(p::DiscreteActionPolicy{OS,AS,T,O,E}, observation::O; deterministic::Bool = false) where {OS,AS,T,O,E}
    if deterministic
        return mode(p, observation)
    end

    cdf = cumsum(probabilities(p, observation))
    s = rand(p.rng)

    # Linear search. Might be too expensive if number of possible actions is big
    for i in eachindex(cdf)
        if s < cdf[i]
            return i - 1
        end
    end
end

function reset!(p::DiscreteActionPolicy; seed::Union{Nothing,UInt} = nothing)
    Ramnet.reset!(p.opt)
    Ramnet.reset!(p.f; seed)
    Random.seed!(p.rng, seed)

    return nothing
end


text:
using PowerSystems
using PowerSimulations
using JuMP

base_dir = string(dirname(dirname(pathof(PowerSystems))))
println(joinpath(base_dir,"data/data_5bus.jl"))
include(joinpath(base_dir,"data/data_5bus.jl"))

battery = [GenericBattery(name = "Bat",
                status = true,
                bus = nodes5[1],
                activepower = 10.0,
                energy = 5.0,
                capacity = (min = 0.0, max = 0.0),
                inputactivepowerlimits = (min = 0.0, max = 50.0),
                outputactivepowerlimits = (min = 0.0, max = 50.0),
                efficiency = (in = 0.90, out = 0.80),
                )];

generators_hg = [
    HydroFix("HydroFix",true,nodes5[2],
        TechHydro(60.0, 15.0, (min = 0.0, max = 60.0), nothing, nothing, nothing, nothing),
        TimeSeries.TimeArray(DayAhead,solar_ts_DA)
    ),
    HydroCurtailment("HydroCurtailment",true,nodes5[3],
        TechHydro(60.0, 10.0, (min = 0.0, max = 60.0), nothing, nothing, (up = 10.0, down = 10.0), nothing),
        100.0,TimeSeries.TimeArray(DayAhead,wind_ts_DA) )
]

sys5b = PowerSystem(nodes5, append!(generators5, generators_hg), loads5_DA, branches5, battery,  100.0)

m = Model()
devices_netinjection =  PowerSimulations.JumpAffineExpressionArray(length(sys5b.buses), sys5b.time_periods)

Pin, Pout, inyection_array = PowerSimulations.powerstoragevariables(m, devices_netinjection,  sys5b.storage, sys5b.time_periods)
Es = PowerSimulations.energystoragevariables(m, sys5b.storage, sys5b.time_periods);
PowerSimulations.powerconstraints(m, sys5b.storage, sys5b.time_periods)
PowerSimulations.energyconstraints(m, sys5b.storage, sys5b.time_periods)
println("defining energy bookkeeping")
PowerSimulations.energybookkeeping(m , sys5b.storage, sys5b.time_periods)

true


text:
module MDAL

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("MDAL was not build properly. Please run Pkg.build(\"MDAL\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

using CEnum

include(joinpath(@__DIR__, "..", "gen", "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libmdal_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libmdal_api.jl"))

# export everything
foreach(names(@__MODULE__, all=true)) do s
   if startswith(string(s), "SOME_PREFIX")
       @eval export $s
   end
end

include("utils.jl")
include("geometry.jl")
include("data.jl")

export MDALMesh
export load

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -17.39516631142554 <= q <= 12.050938789121952)

                     add_NL_constraint(m, :(log(1 + exp(0.1892165831624366 + -0.6797303891528506*log(1 + exp(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3]))) + -0.5563706493225857*log(1 + exp(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3]))) + -0.5881166064998418*log(1 + exp(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3]))) + -0.5324749908347566*log(1 + exp(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))) + log(1 + exp(-0.7087227741366906 + -0.3339546766337098*log(1 + exp(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3]))) + -0.8987724666577948*log(1 + exp(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3]))) + 0.3907505378456695*log(1 + exp(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3]))) + -0.7983638213649877*log(1 + exp(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))) + log(1 + exp(0.16171979073463483 + 0.31171990946083294*log(1 + exp(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3]))) + 0.2757921127232579*log(1 + exp(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3]))) + -0.732639957443523*log(1 + exp(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3]))) + -0.947011017876227*log(1 + exp(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))) + log(1 + exp(-0.17479526548317725 + -0.7121643152938684*log(1 + exp(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3]))) + -0.20573796119512977*log(1 + exp(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3]))) + 0.920351043144791*log(1 + exp(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3]))) + -0.8770524767395065*log(1 + exp(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Decart
    using HorizonSideRobots

    export init, move!, get_coord

    X_COORDINATE=0
    Y_COORDINATE=0

    function init(x=0,y=0)
        global X_COORD, Y_COORD
        X_COORD=0
        Y_COORD=0
    end

    function move!(r,side)
        global X_COORD, Y_COORD
        if side==Nord
            Y_COORD+=1
        elseif side==Sud
            Y_COORD-=1
        elseif side==Ost
            X_COORD+=1
        else
            X_COORD-=1
        end
        HorizonSideRobots.move!(r,side)
    end

    det_coord()=( X_COORD, Y_COORD)
end
using .Decart
function mark_chess(r::Robot,n::Int)

    Decart.init()
    side=Ost
    mark_row(r,side,n)
    while isborder(r,Nord)==false
        Decart.move!(r,Nord)
        side = inverse(side)
        mark_row(r,side,n)
    end
end

function mark_row(r::Robot,size::HorizonSide, n::Integer)
    putmarker_chess!(r,n)
    while isborder!(r,side)==false
        Decart.move!(r,side)
        putmarker_chess!(r,n)
    end
end

function putmarker_chess!(r::Robot,n::Integer)
    x,y = Decart.get_coord()
    x = mod(x,2n)
    y = mod(y,2n)
    if (x in 0:n && y in 0:n)  || (x in n+1:2n && y in n+1:2n)
        putmarker!(r)
    end
end

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


text:
function floydtriangle(rows)
    r = collect(1:div(rows *(rows + 1), 2))
    spacing = Int(ceil(log10(r[end] + 1))) + 1
    for i in 1:rows
        for _ in 1:i
            print(lpad(popfirst!(r), spacing))
        end
        println()
    end
end

floydtriangle(5); println(); floydtriangle(14)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
mutable struct ResponseDecoder
    name::String
    decode::Function
    τ::AbstractArray
end

using ScikitLearn
@sk_import naive_bayes: GaussianNB
using ScikitLearn.CrossValidation: cross_val_score, StratifiedKFold
#@sk_import model_selection: StratifiedKFold
#import StatsBase.score
#function score(model_type, X, y)
#    model = model_type()
#    ScikitLearn.fit!(model, X, y)
#    accuracy = sum(ScikitLearn.predict(model, X) .== y) / length(y)
#    return accuracy
#end

function evaluate(model, X, y)
    return mean(cross_val_score(model, X, y, cv=CrossValidation.StratifiedKFold(y, shuffle=false, random_state=0, n_folds=5)))
end

function pca_score(activity_matrix_split, before_after_identifiers)
    pca = MultivariateStats.fit(PCA, Float64.(activity_matrix_split'))
    X = activity_matrix_split * pca.proj
    y = before_after_identifiers
    return evaluate(GaussianNB(priors=[0.5,0.5]), X, y)
end

function spike_count_score(activity_matrix_split, before_after_identifiers)
    X = sum(activity_matrix_split,dims=2)
    y = before_after_identifiers

    return evaluate(GaussianNB(priors=[0.5,0.5]), X, y)
end

function wavelet_decoder_score2(activity_matrix_split, before_after_identifiers; opts=WIOpts())
    X = activity_matrix_split
    y = before_after_identifiers
    clf = WaveletClassifier(nsurr=opts.nsurr, nscales=opts.nscales,
                            maxcoefs=opts.maxwvcoefs, mincoefs=opts.minwvcoefs,
                            percentile=opts.percentile)
    return evaluate(clf, X, y)
end

function wavelet_decoder_score(activity_matrix_split, before_after_identifiers,
                               opts=WIopts())
    decomposed_matrix = decompose(activity_matrix_split)
    feature_select,(mi,umi)  = WaveletInformation.mi_select(decomposed_matrix,
                                                            before_after_identifiers,opts)
    wavelet_parameters = Dict(
                              "chosen_features" => feature_select,
                              "mutual_information" => mi,
                              "unbiased_mutual_information" => umi
                             )
    X = decomposed_matrix[:,feature_select]
    y = before_after_identifiers

    return evaluate(GaussianNB(priors=[0.5,0.5]), X, y)
end

max_level(X) = floor(Int,log2(size(X,2)))

function decompose(X::Array{T,2}, args...; kwargs...) where T<:Integer
    decompose(Float64.(X), args...; kwargs...)
end

function decompose(X::Array{T,2}, W = wavelet(WT.haar);
                   level = max_level(X)
                  ) where T<:AbstractFloat
    # intialise array
    wX = zeros(size(X)...)
    for i in 1:size(wX,1)
        wX[i,:] = dwt(X[i,:], W, level)
    end
    return wX::Array{T,2}
end


text:
using SCF
using Test


text:
import YAML, TOML, JSON, JSONSchema, Glob
using Harmonie, Test

config_exp  = TOML.parsefile("config/config_exp.toml")

for harmonie_configdir in filter(x -> isdir(joinpath(Harmonie.testbeddir, x)), readdir(Harmonie.testbeddir))
    harmonie_config = TOML.parsefile("harmonie_configurations/$harmonie_configdir.toml")
    # harmonie_config can not be merged because it is missing categories
    #  merged1  = merge(merge, config_exp, harmonie_config)
    @testset "$(basename(harmonie_configdir))" begin
       for testbed_config in readdir(joinpath(Harmonie.testbeddir,harmonie_configdir))
           # merged2 = TOML.parsefile("harmonie_configurations/$harmonie_configdir.toml")
           @test 1==1 # dummy test for now until we can merge with categories
       end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

bat_default(::typeof(bat_initval), ::Val{:rng}, ::Any) = bat_rng()

bat_default(::typeof(bat_initval), ::Val{:algorithm}, ::AnyDensityLike) = InitFromTarget()


text:
function fx(t, x, fx)
    fx .= x
end

function fq(t, q, p, fq)
    fq .= q
end

function fp(t, q, p, fp)
    fp .= q.^2
end

function fϕ(t, x, fϕ)
    fϕ .= 0
end

function gx(t, x, λ, fx)
    fx .= x
end

function gq(t, q, p, λ, fλ)
    fλ .= q
end

function gp(t, q, p, λ, fλ)
    fλ .= q.^2
end

function gϕ(t, q, p, gϕ)
    gϕ .= p - q.^2
end


text:
using Cytof5
import Random
using Flux
include("../../../sims/cb/PreProcess.jl")

using JLD2, FileIO
function loadSingleObj(objPath)
  data = load(objPath)
  return data[collect(keys(data))[1]]
end

Random.seed!(2)

include("../VB.jl")

tau = .005 # FIXME: why can't I do .001?
use_stickbreak = false
noisy_var = 10.0

# SIMULATE_DATA = false
SIMULATE_DATA = true
if SIMULATE_DATA
  L = Dict(0=>3, 1=>3)
  I = 3
  J = 20
  K = 5

  println("Simulate Data")
  N = [8, 1, 2] * 500
  # @time dat = Cytof5.Model.genData(J, N, K, Dict{Int,Int}(L))
  mus=Dict(0 => -[1.0, 2.3, 3.5],
           1 => +[1.0, 2.0, 3.0])
  a_W=rand(K)*10
  a_eta=Dict(z => rand(L[z])*10 for z in 0:1)
  Z = Cytof5.Model.genSimpleZ(J, K)
  @time dat = Cytof5.Model.genData(J=J, N=N, K=K, L=L, Z=Z,
                                   beta=[-9.2, -2.3],
                                   sig2=[0.2, 0.1, 0.3],
                                   mus=mus,
                                   a_W=a_W,
                                   a_eta=a_eta,
                                   sortLambda=false, propMissingScale=0.7,
                                   eps=fill(.005, I))

  I = length(N)
  K_MCMC = 10
  L_MCMC = Dict(false=>5, true=>3)
  priors = VB.Priors(K=K_MCMC, L=L_MCMC, use_stickbreak=use_stickbreak)
  mc = Cytof5.Model.defaultConstants(Cytof5.Model.Data(dat[:y]),
                                     K_MCMC, Dict{Int64,Int64}(L_MCMC),
                                     yQuantiles=[0.0, 0.25, 0.5], pBounds=[.05, .8, .05])
  beta = [mc.beta[:, i] for i in 1:I]
  c = VB.Constants(I, N, J, K_MCMC, L_MCMC, tau, beta, use_stickbreak, noisy_var, priors)
  y = dat[:y]
else
  cbDataPath = "../../sims/cb/data/cytof_cb_with_nan.jld2"
  y = loadSingleObj(cbDataPath)
  K_MCMC=30
  L = Dict(false=>5, true=>3)
  goodColumns, _ = PreProcess.preprocess!(y, maxNanOrNegProp=.9, maxPosProp=.9,
                                          subsample=0.0, rowThresh=-6.0)
  cbData = Cytof5.Model.Data(y)
  mc = Cytof5.Model.defaultConstants(cbData,
                                     K_MCMC, Dict{Int64,Int64}(L),
                                     yQuantiles=[0.0, 0.25, 0.5], pBounds=[.05, .8, .05])
  beta = [mc.beta[:, i] for i in 1:cbData.I]
  priors = VB.Priors(K=mc.K, L=L, use_stickbreak=use_stickbreak)
  c = VB.Constants(cbData.I, cbData.N, cbData.J, mc.K, L,
                   tau, beta, use_stickbreak, noisy_var, priors)
end

if SIMULATE_DATA
  Random.seed!(1)
else
  Random.seed!(0)
end

println("fit model...")
out = VB.fit(y=y, niters=10000, batchsize=2000, c=c, opt=ADAM(1e-2))

# POST PROCESS
using RCall
@rlibrary rcommon
@rlibrary cytof3

# Plot ELBO trace
R"plot"(metrics[:elbo]/sum(c.N), xlab="iter", ylab="elbo", typ="l")

# Number of MC samples
B = 100

println("test rsample of state")
# @time realp, tranp, yout, log_qy = VB.rsample(state, y, c);
@time trace = [VB.rsample(s, y, c)[2] for s in state_hist[1:20:end]];
sig2_trace = hcat([s.sig2.data for s in trace]...)
R"plot"(sig2_trace[1,:], typ="l", xlab="sig2", ylab="trace", col=2, lwd=2,
        ylim=[minimum(sig2_trace), maximum(sig2_trace)]);
R"lines"(sig2_trace[2,:], col=3, lwd=2);
R"lines"(sig2_trace[3,:], col=4, lwd=2);

# Draw MC samples from variational distributions
@time samples = [VB.rsample(state)[2] for i in 1:B];

v = hcat([s.v for s in samples]...).data
v = reshape(v, 1, c.K, length(samples))
H = cat([s.H for s in samples]..., dims=3).data
Z = Int.(v .- H .> 0)
my_image(reshape(mean(Z, dims=3), c.J, c.K))

sig2 = hcat([s.sig2 for s in samples]...).data
if SIMULATE_DATA
  ymin = minimum(minimum.([dat[:sig2], sig2]))
  ymax = maximum(maximum.([dat[:sig2], sig2]))
  R"boxplot"(sig2', typ="l", xlab="", ylab="", ylim=[ymin, ymax]);
  R"abline"(h=dat[:sig2], col="grey", lty=2);
else
  R"boxplot"(sig2', typ="l", xlab="", ylab="");
end

# Plot mu
delta0 = hcat([s.delta0 for s in samples]...).data
mu0 = -cumsum(delta0, dims=1)
delta1 = hcat([s.delta1 for s in samples]...).data
mu1 = cumsum(delta1, dims=1)
mu = [mu0; mu1]

R"boxplot"(mu', typ="l", xlab="", ylab="");
if SIMULATE_DATA
  R"abline"(h=dat[:mus][0], col="grey", lty=2);
  R"abline"(h=dat[:mus][1], col="grey", lty=2);
end
R"abline"(v=c.L[0] + .5, h=0, col="grey");


text:
using Test

@timed_testset "integration/" begin
    include("integration/sequential.jl")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license
"""
    AbstractSparseArray{Tv,Ti,N}

Supertype for `N`-dimensional sparse arrays (or array-like types) with elements
of type `Tv` and index type `Ti`. [`SparseMatrixCSC`](@ref), [`SparseVector`](@ref)
and `SuiteSparse.CHOLMOD.Sparse` are subtypes of this.
"""
abstract type AbstractSparseArray{Tv,Ti,N} <: AbstractArray{Tv,N} end

"""
    AbstractSparseVector{Tv,Ti}

Supertype for one-dimensional sparse arrays (or array-like types) with elements
of type `Tv` and index type `Ti`. Alias for `AbstractSparseArray{Tv,Ti,1}``.
"""
const AbstractSparseVector{Tv,Ti} = AbstractSparseArray{Tv,Ti,1}
"""
    AbstractSparseMatrix{Tv,Ti}

Supertype for two-dimensional sparse arrays (or array-like types) with elements
of type `Tv` and index type `Ti`. Alias for `AbstractSparseArray{Tv,Ti,2}`.
"""
const AbstractSparseMatrix{Tv,Ti} = AbstractSparseArray{Tv,Ti,2}

"""
    issparse(S)

Returns `true` if `S` is sparse, and `false` otherwise.

# Examples
```jldoctest
julia> sv = sparsevec([1, 4], [2.3, 2.2], 10)
10-element SparseVector{Float64,Int64} with 2 stored entries:
  [1 ]  =  2.3
  [4 ]  =  2.2

julia> issparse(sv)
true

julia> issparse(Array(sv))
false
```
"""
issparse(A::AbstractArray) = false
issparse(S::AbstractSparseArray) = true

issparse(S::LinearAlgebra.Symmetric{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.Hermitian{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.LowerTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UnitLowerTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UpperTriangular{<:Any,<:AbstractSparseMatrix}) = true
issparse(S::LinearAlgebra.UnitUpperTriangular{<:Any,<:AbstractSparseMatrix}) = true

indtype(S::AbstractSparseArray{<:Any,Ti}) where {Ti} = Ti

function Base.reinterpret(::Type, A::AbstractSparseArray)
    error("""
          `reinterpret` on sparse arrays is discontinued.
          Try reinterpreting the value itself instead.
          """)
end

# The following two methods should be overloaded by concrete types to avoid
# allocating the I = findall(...)
_sparse_findnextnz(v::AbstractSparseArray, i::Integer) = (I = findall(!iszero, v); n = searchsortedfirst(I, i); n<=length(I) ? I[n] : nothing)
_sparse_findprevnz(v::AbstractSparseArray, i::Integer) = (I = findall(!iszero, v); n = searchsortedlast(I, i);  !iszero(n)   ? I[n] : nothing)

function findnext(f::typeof(!iszero), v::AbstractSparseArray, i::Integer)
    j = _sparse_findnextnz(v, i)
    while j !== nothing && !f(v[j])
        j = _sparse_findnextnz(v, j+1)
    end
    return j
end

function findprev(f::typeof(!iszero), v::AbstractSparseArray, i::Integer)
    j = _sparse_findprevnz(v, i)
    while j !== nothing && !f(v[j])
        j = _sparse_findprevnz(v, j-1)
    end
    return j
end

"""
    findnz(A)

Return a tuple `(I, J, V)` where `I` and `J` are the row and column indices of the stored
("structurally non-zero") values in sparse matrix `A`, and `V` is a vector of the values.

# Examples
```jldoctest
julia> A = sparse([1 2 0; 0 0 3; 0 4 0])
3×3 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1, 1]  =  1
  [1, 2]  =  2
  [3, 2]  =  4
  [2, 3]  =  3

julia> findnz(A)
([1, 1, 3, 2], [1, 2, 2, 3], [1, 2, 4, 3])
```
"""
function findnz end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
using Test

import Tar
import Pkg.GitTools: tree_hash

function gen_file(file::String, size::Int)
    open(file, write=true) do io
        for i = 1:size; write(io, i % UInt8); end
    end
end

function make_test_tarball(tar_create::Function = Tar.create)
    root = mktempdir()
    let i = 0, paths = String[]
        a_lengths = [0, 10, 154, 155, 156, 255]
        @static if Sys.iswindows()
            filter!(a -> 2 + a + length(root) < 260, a_lengths)
        end
        for a in a_lengths
            dir = joinpath(root, "d"^a)
            push!(paths, normpath("..", dir))
            a > 0 && mkdir(dir)
            b_lengths = [10, 99, 100, 101, 255]
            # don't attempt to create file paths > 260 on windows
            @static if Sys.iswindows()
                filter!(b -> 4 + b + a + length(root) < 260, b_lengths)
            end
            for b in b_lengths
                for s in [0, 511, 512, 513, 1000]
                    f = rpad("$s-", b, "f")
                    x = rpad("$s-", b, "x")
                    file = joinpath(dir, f)
                    exec = joinpath(dir, x)
                    gen_file(file, s)
                    gen_file(exec, s)
                    chmod(file, 0o644)
                    chmod(exec, 0o755)
                    push!(paths, file)
                    push!(paths, exec)
                end
                dir′ = joinpath(dir, "s"^b)
                mkpath(dir′)
                push!(paths, dir′)
                link = joinpath(dir, "l"^b)
                target = relpath(paths[i += 1], link)
                symlink(target, link)
                push!(paths, link)
                broken = joinpath(dir, "b"^b)
                if target != "."
                    symlink(chop(target), broken)
                    push!(paths, broken)
                end
            end
        end
    end
    hash = tree_hash(root)
    tarball = tar_create(root)
    rm(root, force=true, recursive=true)
    return tarball, hash
end

@testset "empty tarball" begin
    dir = mktempdir()
    tarball = Tar.create(dir)
    rm(dir, recursive=true)
    @test Tar.list(tarball) == [Tar.Header(".", :directory, 0o755, 0, "")]
    dir = Tar.extract(tarball)
    @test isempty(readdir(dir))
    rm(dir, recursive=true)
    open(tarball, append=true) do io
        write(io, zeros(UInt8, 512))
    end
    dir = Tar.extract(tarball)
    @test isempty(readdir(dir))
    rm(dir, recursive=true)
end

function check_tree_hash(hash::Vector{UInt8}, root::AbstractString)
    @test tree_hash(root) == hash
    rm(root, recursive=true)
end

@testset "test tarball" begin
    tarball, hash = make_test_tarball()
    @testset "Tar.list & check properties" begin
        headers = Tar.list(tarball)
        @test issorted(headers, by = hdr -> hdr.path)
        for hdr in headers
            @test !isempty(hdr.path)
            @test hdr.type in (:file, :directory, :symlink)
            if hdr.type == :file
                @test hdr.path[end] != '/'
                @test hdr.mode in (0o644, 0o755)
                @test isempty(hdr.link)
            elseif hdr.type == :directory
                @test hdr.path[end] == '/'
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test isempty(hdr.link)
            elseif hdr.type == :symlink
                @test hdr.path[end] != '/'
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test !isempty(hdr.link)
            end
        end
        @testset "Tar.list from IO, process, pipeline" begin
            @test headers == open(Tar.list, tarball)
            @test headers == open(Tar.list, `cat $tarball`)
            @test headers == open(Tar.list, pipeline(`bzip2 -c -9 $tarball`, `bzcat`))
        end
    end
    # skip `tar` tests when it doesn't exist or when we're on windows
    if Sys.which("tar") != nothing && !Sys.iswindows()
        @testset "extract with `tar` command" begin
            root = mktempdir()
            run(`tar -C $root -xf $tarball`)
            check_tree_hash(hash, root)
        end
    end
    @testset "Tar.extract" begin
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    @testset "Tar.extract from IO, process, pipeline" begin
        root = open(Tar.extract, tarball)
        check_tree_hash(hash, root)
        root = open(Tar.extract, `cat $tarball`)
        check_tree_hash(hash, root)
        root = open(Tar.extract, pipeline(`bzip2 -c -9 $tarball`, `bzcat`))
        check_tree_hash(hash, root)
    end
    open(tarball, append=true) do io
        write(io, zeros(UInt8, 512))
    end
    @testset "Tar.extract with trailing zeros" begin
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    rm(tarball)
end

if Sys.which("gtar") != nothing && !Sys.iswindows()
    @testset "POSIX extensions" begin
        # make a test POSIX tarball with `gtar` instead of Tar.create
        tarball, hash = make_test_tarball() do root
            tarball, io = mktemp(); close(io)
            run(`gtar --format=posix -C $root -cf $tarball .`)
            return tarball
        end
        # TODO: check that extended headers contain `mtime` etc.
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    @testset "GNU extensions" begin
        # make a test GNU tarball with `gtar` instead of Tar.create
        tarball, hash = make_test_tarball() do root
            tarball, io = mktemp(); close(io)
            run(`gtar --format=gnu -C $root -cf $tarball .`)
            return tarball
        end
        hdrs = Tar.list(tarball, raw=true)
        # test that both long link and long name entries are created
        @test any(h.path == "././@LongLink" && h.type == :L for h in hdrs)
        @test any(h.path == "././@LongLink" && h.type == :K for h in hdrs)
        # test that Tar can extract these GNU entries correctly
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
end

@testset "symlink attacks" begin
    # not dangerous but still not allowed
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("dir", :directory, 0o755, 0, ""))
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "dir"))
    Tar.write_header(io, Tar.Header("link/target", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # attempt to write through relative link out of root
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "../target"))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # attempt to write through absolute link
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # same attack with some obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("./link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
    # same attack with different obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "/tmp"))
    Tar.write_header(io, Tar.Header("dir/../link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    rm(tarball)
end

!Sys.iswindows() &&
@testset "symlink overwrite" begin
    # allowable and should work
    @testset "allow overwriting a symlink" begin
        tarball₁, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :symlink, 0o755, 0, "/tmp"))
        Tar.write_header(io, Tar.Header("path", :file, 0o644, 0, ""))
        close(io)
        tree₁ = Tar.extract(tarball₁)
        tarball₂, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :file, 0o644, 0, ""))
        close(io)
        tree₂ = Tar.extract(tarball₂)
        @test tree_hash(tree₁) == tree_hash(tree₂)
        rm(tree₁, recursive=true)
        rm(tree₂, recursive=true)
        rm(tarball₁)
        rm(tarball₂)
    end
    @testset "allow write into directory overwriting a symlink" begin
        # make sure "path" is removed from links set
        tarball₁, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :symlink, 0o755, 0, "/tmp"))
        Tar.write_header(io, Tar.Header("path", :directory, 0o755, 0, ""))
        Tar.write_header(io, Tar.Header("path/file", :file, 0o644, 0, ""))
        close(io)
        tree₁ = Tar.extract(tarball₁)
        tarball₂, io = mktemp()
        Tar.write_header(io, Tar.Header("path/file", :file, 0o644, 0, ""))
        close(io)
        tree₂ = Tar.extract(tarball₂)
        @test tree_hash(tree₁) == tree_hash(tree₂)
        rm(tree₁, recursive=true)
        rm(tree₂, recursive=true)
        rm(tarball₁)
        rm(tarball₂)
    end
end

function make_test_dir(gen_skip::Bool=false)
    dir = mktempdir()
    touch(joinpath(dir, "file"))
    if gen_skip
        touch(joinpath(dir, "file.skip"))
    end
    mkdir(joinpath(dir, "dir"))
    touch(joinpath(dir, "dir", "file"))
    mkdir(joinpath(dir, "empty"))
    if gen_skip
        touch(joinpath(dir, "dir", "file.skip"))
    end
    if gen_skip
        mkdir(joinpath(dir, "dir.skip"))
        touch(joinpath(dir, "dir.skip", "file"))
    end
    if !Sys.iswindows()
        symlink("file", joinpath(dir, "link"))
        if gen_skip
            symlink("file", joinpath(dir, "link.skip"))
        end
    end
    return dir
end

const test_dir_paths = ["dir/file", "empty/", "file", "link"]
Sys.iswindows() && pop!(test_dir_paths)

@testset "API: create" begin
    dir = make_test_dir()
    @test !any(splitext(name)[2] == ".skip" for name in readdir(dir))
    # create(dir::String)
    tarball = Tar.create(dir)
    bytes = read(tarball)
    @test isfile(tarball)
    rm(tarball)
    # create(dir::String, tarball::String)
    tarball = tempname()
    Tar.create(dir, tarball)
    @test read(tarball) == bytes
    rm(tarball)
    # create(dir::String, tarball::IO)
    mktemp() do tarball, io
        Tar.create(dir, tarball)
        close(io)
        @test read(tarball) == bytes
    end
    rm(dir, recursive=true)
    # test predicate versions
    dir = make_test_dir(true)
    @test any(splitext(name)[2] == ".skip" for name in readdir(dir))
    predicate = path -> splitext(path)[2] != ".skip"
    # create(predicate::Function, dir::String)
    tarball = Tar.create(predicate, dir)
    @test read(tarball) == bytes
    rm(tarball)
    # create(predicate::Function, dir::String, tarball::String)
    tarball = tempname()
    Tar.create(predicate, dir, tarball)
    @test read(tarball) == bytes
    rm(tarball)
    # create(predicate::Function, dir::String, tarball::IO)
    mktemp() do tarball, io
        Tar.create(predicate, dir, tarball)
        close(io)
        @test read(tarball) == bytes
    end
    rm(dir, recursive=true)
end

@testset "API: list" begin
    dir = make_test_dir()
    tarball = Tar.create(dir)
    rm(dir, recursive=true)
    # list(tarball::String)
    headers = Tar.list(tarball)
    @test test_dir_paths == [hdr.path for hdr in headers]
    # list(tarball::IO)
    headers = open(Tar.list, tarball)
    @test test_dir_paths == [hdr.path for hdr in headers]
    # add a sketchy entry to tarball
    open(tarball, append=true) do io
        Tar.write_header(io, Tar.Header("/bad", :file, 0o644, 0, ""))
    end
    paths = push!(copy(test_dir_paths), "/bad")
    # list(tarball::String; strict=true|false)
    @test_throws ErrorException Tar.list(tarball)
    @test_throws ErrorException Tar.list(tarball, strict=true)
    headers = Tar.list(tarball, strict=false)
    @test paths == [hdr.path for hdr in headers]
    # list(tarball::IO; strict=true|false)
    @test_throws ErrorException open(Tar.list, tarball)
    @test_throws ErrorException open(tarball) do io
        Tar.list(io, strict=true)
    end
    headers = open(tarball) do io
        Tar.list(io, strict=false)
    end
    @test paths == [hdr.path for hdr in headers]
end

@testset "API: extract" begin
    dir = make_test_dir()
    hash = tree_hash(dir)
    tarball = Tar.create(dir)
    rm(dir, recursive=true)

    # extract(tarball::String)
    dir = Tar.extract(tarball)
    check_tree_hash(hash, dir)
    # extract(tarball::String, dir::String) — non-existent
    dir = tempname()
    Tar.extract(tarball, dir)
    check_tree_hash(hash, dir)
    # extract(tarball::String, dir::String) — existent, empty
    dir = mktempdir()
    Tar.extract(tarball, dir)
    check_tree_hash(hash, dir)
    # extract(tarball::String, dir::String) — non-directory (error)
    dir = tempname()
    touch(dir)
    @test_throws ErrorException Tar.extract(tarball, dir)
    rm(dir)
    # extract(tarball::String, dir::String) — non-empty directory (error)
    dir = mktempdir()
    touch(joinpath(dir, "file"))
    @test_throws ErrorException Tar.extract(tarball, dir)
    rm(dir, recursive=true)

    # extract(tarball::IO)
    dir = open(Tar.extract, tarball)
    check_tree_hash(hash, dir)
    # extract(tarball::IO, dir::String) — non-existent
    dir = tempname()
    open(tarball) do io
        Tar.extract(io, dir)
    end
    check_tree_hash(hash, dir)
    # extract(tarball::IO, dir::String) — existent, empty
    dir = mktempdir()
    open(tarball) do io
        Tar.extract(io, dir)
    end
    check_tree_hash(hash, dir)
    # extract(tarball::IO, dir::String) — non-directory (error)
    dir = tempname()
    touch(dir)
    @test_throws ErrorException open(tarball) do io
        Tar.extract(io, dir)
    end
    rm(dir)
    # extract(tarball::IO, dir::String) — non-empty directory (error)
    dir = mktempdir()
    touch(joinpath(dir, "file"))
    @test_throws ErrorException open(tarball) do io
        Tar.extract(io, dir)
    end
    rm(dir, recursive=true)
end


text:
#=
let E be a symmetric matrix sparsity pattern:
(1) find sparse PSD matrix with given sparsity pattern, "nearest" to A
    max_X tr(A, X) :
    tr(X) = 1
    X in intersect(S_+^n, S^n(E))
(2) find sparse PSD-completable matrix with given sparsity pattern, "nearest" to A
    max_X tr(A, X) :
    tr(X) = 1
    X in proj_E(S_+^n)

adapted from
"Decomposition methods for sparse matrix nearness problems"
(2015) by Sun & Vandenberghe
=#

using SparseArrays

struct NearestPSDJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    side::Int
    use_completable::Bool # solve problem (2) above, else solve problem (1)
    use_chordal_sparsity::Bool # use a chordal sparsity pattern
    use_sparsepsd::Bool # use sparse PSD cone formulation, else dense PSD cone
    use_cholmod_impl::Bool # use sparsePSD CHOLMOD implementation, else dense
end

function build(inst::NearestPSDJuMP{T}) where {T <: Float64}
    side = inst.side
    sparsity = min(3.0 / side, 1.0) # sparsity factor (before chordal extension)

    # generate random symmetric A (indefinite) with sparsity pattern E (nonchordal)
    A = tril!(sprandn(side, side, sparsity)) + Diagonal(randn(side))
    if inst.use_chordal_sparsity
        # compute a (heuristic) chordal extension of A using CHOLMOD functions
        copyto!(A, I)
        A = sparse(cholesky(Symmetric(A, :L)).L)
        (row_idxs, col_idxs, _) = findnz(A)
        A_vals = randn(length(row_idxs))
    else
        (row_idxs, col_idxs, A_vals) = findnz(A)
    end
    diag_idxs = findall(row_idxs .== col_idxs)

    model = JuMP.Model()

    if inst.use_sparsepsd || !inst.use_completable
        JuMP.@variable(model, X[1:length(row_idxs)])
        JuMP.@objective(model, Max, 2 * dot(A_vals, X) -
            sum(A_vals[k] * X[k] for k in diag_idxs)) # tr(A, X)
        JuMP.@constraint(model, sum(X[diag_idxs]) == 1) # tr(X) == 1

        if inst.use_sparsepsd
            rt2 = sqrt(2)
            X_scal = [X[k] * (row_idxs[k] == col_idxs[k] ? 1.0 : rt2)
                for k in eachindex(X)]
            impl = (inst.use_cholmod_impl ? Cones.PSDSparseCholmod :
                Cones.PSDSparseDense)
            cone = Hypatia.PosSemidefTriSparseCone{impl, T, T}
            JuMP.@constraint(model, X_scal in cone(side, row_idxs,
                col_idxs, inst.use_completable))
        else
            X_sparse = sparse(row_idxs, col_idxs, X)
            JuMP.@SDconstraint(model, Symmetric(Matrix(X_sparse), :L) >= 0)
        end
    else
        @assert inst.use_completable
        JuMP.@variable(model, X[1:side, 1:side], PSD)
        # tr(A, X):
        JuMP.@objective(model, Max, 2 * sum(X[row_idxs[k], col_idxs[k]] * A_vals[k]
            for k in eachindex(row_idxs)) - dot(A_vals[diag_idxs], diag(X)))
        JuMP.@constraint(model, tr(X) == 1) # tr(X) == 1
    end

    return model
end


text:
# Stream management

export
    CuStream, CuDefaultStream, CuStreamLegacy, CuStreamPerThread,
    priority, priority_range, synchronize

mutable struct CuStream
    handle::CUstream
    ctx::CuContext
end

Base.unsafe_convert(::Type{CUstream}, s::CuStream) = s.handle

Base.:(==)(a::CuStream, b::CuStream) = a.handle == b.handle
Base.hash(s::CuStream, h::UInt) = hash(s.handle, h)

@enum_without_prefix CUstream_flags_enum CU_

"""
    CuStream(; flags=STREAM_DEFAULT, priority=nothing)

Create a CUDA stream.
"""
function CuStream(; flags::CUstream_flags=STREAM_DEFAULT,
                    priority::Union{Nothing,Integer}=nothing)
    handle_ref = Ref{CUstream}()
    if priority === nothing
        cuStreamCreate(handle_ref, flags)
    else
        priority in priority_range() || throw(ArgumentError("Priority is out of range"))
        cuStreamCreateWithPriority(handle_ref, flags, priority)
    end

    ctx = CuCurrentContext()
    obj = CuStream(handle_ref[], ctx)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(s::CuStream)
    if isvalid(s.ctx)
        cuStreamDestroy_v2(s)
    end
end

function Base.show(io::IO, stream::CuStream)
    print(io, "CuStream(")
    @printf(io, "%p", stream.handle)
    print(io, ", ", stream.ctx, ")")
end

"""
    CuDefaultStream()

Return the default stream.

!!! note

    It is generally better to use `stream()` to get a stream object that's local to the
    current task. That way, operations scheduled in other tasks can overlap.
"""
@inline CuDefaultStream() = CuStream(convert(CUstream, C_NULL), CuContext(C_NULL))

"""
    CuStreamLegacy()

Return a special object to use use an implicit stream with legacy synchronization behavior.

You can use this stream to perform operations that should block on all streams (with the
exception of streams created with `STREAM_NON_BLOCKING`). This matches the old pre-CUDA 7
global stream behavior.
"""
@inline CuStreamLegacy() = CuStream(convert(CUstream, 1), CuContext(C_NULL))

"""
    CuStreamPerThread()

Return a special object to use an implicit stream with per-thread synchronization behavior.
This stream object is normally meant to be used with APIs that do not have per-thread
versions of their APIs (i.e. without a `ptsz` or `ptds` suffix).

!!! note

    It is generally not needed to use this type of stream. With CUDA.jl, each task already
    gets its own non-blocking stream, and multithreading in Julia is typically
    accomplished using tasks.
"""
@inline CuStreamPerThread() = CuStream(convert(CUstream, 2), CuContext(C_NULL))

"""
    query(s::CuStream)

Return `false` if a stream is busy (has task running or queued)
and `true` if that stream is free.
"""
function query(s::CuStream)
    res = unsafe_cuStreamQuery(s)
    if res == ERROR_NOT_READY
        return false
    elseif res == SUCCESS
        return true
    else
        throw_api_error(res)
    end
end

"""
    synchronize()
    synchronize([stream::CuStream]; blocking=true)

Wait until `stream` has finished executing, with `stream` defaulting to the stream
associated with the current Julia task. If `blocking` is true (the default), Julia will be
asked to yield to any other scheduled task, otherwise it will spin indefinitely.

See also: [`synchronize_all`](@ref)
"""
function synchronize(s::CuStream=stream(); blocking::Bool=true)
    # TODO: exponential back-off and sleep?
    while !query(s)
        if blocking
            yield()
        else
            ccall(:jl_cpu_pause, Cvoid, ())
            # Temporary solution before we have gc transition support in codegen.
            ccall(:jl_gc_safepoint, Cvoid, ())
        end
    end

    check_exceptions()
end

"""
    priority_range()

Return the valid range of stream priorities as a `StepRange` (with step size  1). The lower
bound of the range denotes the least priority (typically 0), with the upper bound
representing the greatest possible priority (typically -1).
"""
function priority_range()
    least_ref = Ref{Cint}()
    greatest_ref = Ref{Cint}()
    cuCtxGetStreamPriorityRange(least_ref, greatest_ref)
    step = least_ref[] < greatest_ref[] ? 1 : -1
    return least_ref[]:Cint(step):greatest_ref[]
end

"""
    priority_range(s::CuStream)

Return the priority of a stream `s`.
"""
function priority(s::CuStream)
    priority_ref = Ref{Cint}()
    cuStreamGetPriority(s, priority_ref)
    return priority_ref[]
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "coreutils"
version = v"9.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/coreutils/coreutils-$(version.major).$(version.minor).tar.xz", "ce30acdf4a41bc5bb30dd955e9eaa75fa216b4e3deb08889ed32433c7b3b97ce")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/coreutils-9.*

# Fix `configure: error: you should not run configure as root (set FORCE_UNSAFE_CONFIGURE=1 in environment to bypass this check)`
if [[ ${target} == x86_64-linux-musl* ]]; then
    export FORCE_UNSAFE_CONFIGURE=1
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=p -> !(Sys.islinux(p) | Sys.isfreebsd(p)))

# The products that we will ensure are always built
products = [
    ExecutableProduct("b2sum", :b2sum),
    ExecutableProduct("timeout", :timeout),
    ExecutableProduct("chown", :chown_bin),
    ExecutableProduct("factor", :factor),
    ExecutableProduct("sha512sum", :sha512sum),
    ExecutableProduct("pathchk", :pathchk),
    ExecutableProduct("fmt", :fmt),
    ExecutableProduct("sort", :sort_bin),
    ExecutableProduct("sha1sum", :sha1sum),
    ExecutableProduct("yes", :yes),
    ExecutableProduct("groups", :groups),
    ExecutableProduct("uname", :uname),
    ExecutableProduct("env", :env),
    ExecutableProduct("mv", :mv_bin),
    ExecutableProduct("sync", :sync),
    ExecutableProduct("base64", :base64),
    ExecutableProduct("truncate", :truncate_bin),
    ExecutableProduct("tsort", :tsort),
    ExecutableProduct("numfmt", :numfmt),
    ExecutableProduct("nice", :nice),
    ExecutableProduct("sleep", :sleep_bin),
    ExecutableProduct("stat", :stat_bin),
    ExecutableProduct("link", :link),
    ExecutableProduct("users", :users),
    ExecutableProduct("printenv", :printenv),
    ExecutableProduct("head", :head),
    ExecutableProduct("split", :split_bin),
    ExecutableProduct("sha224sum", :sha224sum),
    ExecutableProduct("shuf", :shuf),
    ExecutableProduct("chroot", :chroot),
    ExecutableProduct("csplit", :csplit),
    ExecutableProduct("ptx", :ptx),
    LibraryProduct("libstdbuf", :libstdbuf, "libexec/coreutils"),
    ExecutableProduct("dir", :dir),
    ExecutableProduct("chmod", :chmod_bin),
    ExecutableProduct("vdir", :vdir),
    ExecutableProduct("test", :test),
    ExecutableProduct("true", :true_bin),
    ExecutableProduct("sha256sum", :sha256sum),
    ExecutableProduct("expr", :expr),
    ExecutableProduct("basename", :basename_bin),
    ExecutableProduct("pwd", :pwd_bin),
    ExecutableProduct("unexpand", :unexpand),
    ExecutableProduct("cp", :cp_bin),
    ExecutableProduct("ln", :ln),
    ExecutableProduct("runcon", :runcon),
    ExecutableProduct("kill", :kill_bin),
    ExecutableProduct("mkfifo", :mkfifo),
    ExecutableProduct("tac", :tac),
    ExecutableProduct("echo", :echo),
    ExecutableProduct("od", :od),
    ExecutableProduct("logname", :logname),
    ExecutableProduct("rm", :rm_bin),
    ExecutableProduct("pinky", :pinky),
    ExecutableProduct("pr", :pr),
    ExecutableProduct("chcon", :chcon),
    ExecutableProduct("tty", :tty),
    ExecutableProduct("touch", :touch_bin),
    ExecutableProduct("df", :df),
    ExecutableProduct("comm", :comm),
    ExecutableProduct("nohup", :nohup),
    ExecutableProduct("date", :date),
    ExecutableProduct("seq", :seq),
    ExecutableProduct("hostid", :hostid),
    ExecutableProduct("sum", :sum_bin),
    ExecutableProduct("cut", :cut),
    ExecutableProduct("readlink", :readlink_bin),
    ExecutableProduct("realpath", :realpath_bin),
    ExecutableProduct("md5sum", :md5sum),
    ExecutableProduct("install", :install),
    ExecutableProduct("mkdir", :mkdir_bin),
    ExecutableProduct("cksum", :cksum),
    ExecutableProduct("tail", :tail_bin),
    ExecutableProduct("rmdir", :rmdir),
    ExecutableProduct("unlink", :unlink_bin),
    ExecutableProduct("whoami", :whoami),
    ExecutableProduct("false", :false_bin),
    ExecutableProduct("[", :left_bracket),
    ExecutableProduct("mknod", :mknod),
    ExecutableProduct("fold", :fold),
    ExecutableProduct("stty", :stty),
    ExecutableProduct("tr", :tr),
    ExecutableProduct("cat", :cat_bin),
    ExecutableProduct("chgrp", :chgrp),
    ExecutableProduct("ls", :ls),
    ExecutableProduct("shred", :shred),
    ExecutableProduct("who", :who),
    ExecutableProduct("tee", :tee),
    ExecutableProduct("wc", :wc),
    ExecutableProduct("base32", :base32),
    ExecutableProduct("printf", :printf),
    ExecutableProduct("join", :join_bin),
    ExecutableProduct("id", :id),
    ExecutableProduct("mktemp", :mktemp_bin),
    ExecutableProduct("uptime", :uptime),
    ExecutableProduct("stdbuf", :stdbuf),
    ExecutableProduct("uniq", :uniq),
    ExecutableProduct("dd", :dd),
    ExecutableProduct("nl", :nl),
    ExecutableProduct("dircolors", :dircolors),
    ExecutableProduct("basenc", :basenc),
    ExecutableProduct("du", :du),
    ExecutableProduct("sha384sum", :sha384sum),
    ExecutableProduct("nproc", :nproc),
    ExecutableProduct("dirname", :dirname_bin),
    ExecutableProduct("expand", :expand),
    ExecutableProduct("paste", :paste)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil