text:
@testset "wrappers" begin

    @testset "ActionTransformedEnv" begin
        env = TigerProblemEnv(; rng=StableRNG(123))
        env′ = ActionTransformedEnv(
            env;
            action_space_mapping=x -> Base.OneTo(3),
            action_mapping=i -> action_space(env)[i],
        )

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)
    end

    @testset "DefaultStateStyleEnv" begin
        rng = StableRNG(123)
        env = TigerProblemEnv(; rng=rng)
        S = InternalState{Int}()
        env′ = DefaultStateStyleEnv{S}(env)
        @test DefaultStateStyle(env′) === S

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)

        # fix #129
        env = TicTacToeEnv()
        s_env = state(env)
        @test s_env ∈ state_space(env)

        E = DefaultStateStyleEnv{Observation{Int}()}(env)  # error probably somewhere here
        s = state(E)
        @test s isa Int
        @test s ∈ state_space(E)
    end

    @testset "MaxTimeoutEnv" begin
        rng = StableRNG(123)
        env = TigerProblemEnv(; rng=rng)
        n = 100
        env′ = MaxTimeoutEnv(env, n)

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)

        while !is_terminated(env′)
            env′(:listen)
            n -= 1
            @test n >= 0
        end

        reset!(env′)
        @test env′.current_t == 1
        @test is_terminated(env′) == false
    end

    @testset "RewardOverriddenEnv" begin
        rng = StableRNG(123)
        env = TigerProblemEnv(; rng=rng)
        env′ = RewardOverriddenEnv(env, x -> sign(x))

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)

        while !is_terminated(env′)
            env′(rand(rng, legal_action_space(env′)))
            @test reward(env′) ∈ (-1, 0, 1)
        end
    end

    @testset "StateCachedEnv" begin
        rng = StableRNG(123)
        env = CartPoleEnv(; rng=rng)
        env′ = StateCachedEnv(env)

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)

        while !is_terminated(env′)
            env′(rand(rng, legal_action_space(env′)))
            s1 = state(env)
            s2 = state(env)
            @test s1 === s2
        end
    end

    @testset "StateTransformedEnv" begin
        rng = StableRNG(123)
        env = TigerProblemEnv(; rng=rng)
        # S = (:door1, :door2, :door3, :none)
        # env′ = StateTransformedEnv(env, state_mapping=s -> s+1)
        # RLBase.state_space(env::typeof(env′), ::RLBase.AbstractStateStyle, ::Any) = S

        # RLBase.test_interfaces!(env′)
        # RLBase.test_runnable!(env′)
    end

    @testset "StochasticEnv" begin
        env = KuhnPokerEnv()
        rng = StableRNG(123)
        env′ = StochasticEnv(env; rng=rng)

        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)
    end

    @testset "SequentialEnv" begin
        env =  RockPaperScissorsEnv()
        env′ = SequentialEnv(env)
        RLBase.test_interfaces!(env′)
        RLBase.test_runnable!(env′)

        # https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/issues/393
        @test RLBase.current_player(env′) == 1
        @test action_space(env′) == action_space(env′, 1)

        env′('💎')

        @test RLBase.current_player(env′) == 2
        @test is_terminated(env′) == false

        reset!(env′)

        @test RLBase.current_player(env′) == 1
        @test reward(env′) == 0
        @test is_terminated(env′) == false

        env′('💎')

        @test reward(env′) == 0
        @test reward(env′, 1) == 0
        @test reward(env′, 2) == 0
        @test is_terminated(env′) == false

        env′('📃')

        @test reward(env′, 1) == -1
        @test reward(env′, 2) == 1
        @test is_terminated(env′) == true

        reset!(env′)
        @test reward(env′) == 0
        @test reward(env′, 1) == 0
        @test reward(env′, 2) == 0
        @test is_terminated(env′) == false
        @test RLBase.current_player(env′) == 1
    end
end


text:
# ---
# title: Max min filters
# cover: assets/max_min.gif
# author: Deeptendu Santra
# date: 2020-10-08
# ---

# In this tutorial we see how can we can effectively use max and min filter to distinguish
# between smooth and texture edges  in grayscale images.

# We will use the [`mapwindow`](@ref) function in `ImageFiltering.jl` which provides a general
# functionality to apply any function to the window around each pixel.

using ImageCore, ImageShow, ImageFiltering
using TestImages

img = Gray.(testimage("house");)      # Original Image

# We can use the `minimum` function to compute the minimum of the grayscale values in the given
# matrix or array. For example:
minimum([Gray(0.7),Gray(0.5),Gray(0.0)]) # Should return Gray(0.0) i.e black.
#
filter_size = (15, 15)
## Using the `mapwindow` function, we create an image of the local minimum.
## `mapwindow` maps the given function over a moving window of given size.
img_min = mapwindow(minimum, img, filter_size)
## Similarly for maximum
img_max = mapwindow(maximum, img, filter_size)
## The max(min) filter
img_max_min = mapwindow(maximum, img_min, filter_size)
## The min(max) filter
img_min_max = mapwindow(minimum, img_max, filter_size)
mosaicview(img_min, img_max, img_max_min, img_min_max; nrow=1)

# Now that we are done with the basic filtered images, we proceed to the next part
# which is edge detection using these filters.

# For edge detection, we need to define thresholds for our image. The threshold is an important tool
# to binarize a grayscale image. The threshold value for a given pixel practically
# decides if the pixel visible or not in the output image. However, appyling a global threshold might not consider the
# variation of colors/brightness within the image. Thus we consider an adaptive type theresholding method
# here using max/min results.

# The max(min) and min(max) filter effectively follows the smooth edges of the image. Therefore,
# their average also follows the smooth parts of the image. If we use this image as a threshold for the original
# image, the smooth parts of the original image will get filtered out, leaving only the texture and/or noise behind.
# So we can use the average of `img_max_min` and `img_min_max` as the texture threshold.

# The average of min and max filters also gives us the smooth edges, but it also includes the noise in the image.
# So using average of `img_min` and `img_max` as a threshold will yeild only texture.

img_texture_noise_threshold = (img_max_min + img_max_min) ./ 2
img_texture_threshold = (img_max + img_min) ./ 2
mosaicview(img_texture_noise_threshold, img_texture_threshold; nrow=1)
#
## The Dynamic Gist is obtained by subtracting the img_texture_threshold from the original image.
## The filtered image gives us the texture of the image.
img_dynamic_gist = img - img_texture_threshold
## The Texture Gist is obtained by subtracting img_texture_noise_threshold from the original image.
## The filtered image gives us the texture along with the noise of the image.
img_texture_gist = img - img_texture_noise_threshold
mosaicview(img_dynamic_gist, img_texture_gist; nrow=1)
#
## We extract out the smooth/ramp parts of the image.
## The darker section of the image consist of the ramp edges. The brighter pixels are mostly noise.
ramp = img_dynamic_gist - img_texture_gist
## Filtered-out edges
edge = img_max - img_min
## Smoothed-out version of edge
edge_smoothed = img_min_max - img_max_min
mosaicview(img, ramp, edge, edge_smoothed; nrow=2)

# # References

# Verbeek, P. W., Vrooman, H. A., & Van Vliet, L. J. (1988). [Low-level image processing by max-min filters](https://core.ac.uk/download/pdf/194053536.pdf). Signal Processing, 15(3), 249-258.

## save covers #src
using ImageMagick #src
mkpath("assets")  #src
ImageMagick.save("assets/max_min.gif", cat(img, edge, edge_smoothed, map(clamp01nan, ramp);dims=3);fps=1) #src

ext: jl


text:
module JPVT

using Roots

export R

const R= 8.314 # J/(mol.K) gas constant

end # module


text:
module List
export List!, Cons, Nil
using MLStyle
using MLStyle.Data
import Base: ^, map, foreach, iterate, isempty

@data List!{T} begin
    Nil{T}
    Cons{T}(head :: T, tail :: List!{T})
end

# not to overload convert,

Base.convert(::Type{List!{G}}, it :: Nil{T}) where {G, T <: G} = Nil{G}()

Base.convert(::Type{List!{G}}, it :: Cons{T}) where {G, T <: G} =
    Cons{G}(
        convert(G, it.head),
        convert(List!{G}, it.tail))

Cons(head :: T, tail :: List!{G}) where {G, T} =
    Cons{typejoin(T, G)}(head, tail)

function (^)(head :: T, tail :: List!{G}) where {G, T}
    Cons{typejoin(T, G)}(head, tail)
end

Nil() = Nil{Any}()

Base.iterate(lst :: List!{T}) where T =
    iterate(lst, lst)

function Base.iterate(::List!{T}, lst :: List!{T})::Union{Nothing, Tuple{T, List!{T}}} where T
    @match lst begin
        Nil{T}()      => nothing
        Cons{T}(a, b) => (a, b)
    end
end

function Base.map(fn :: Fun{T, R}, lst :: List!{T}) :: List!{R} where {T, R}
    ret = Nil{R}()
    for each in lst
       ret = Cons{R}(fn(each), ret)
    end
    ret
end

function Base.map(fn :: Function, lst :: List!{T}) :: List! where T

    ret = Nil()
    for each in lst
        ret = Cons(fn(each), ret)
    end
    ret

end

function Base.foreach(fn :: Fun{T, R}, lst :: List!{T}) :: Nothing where {T, R}
    for each in lst
        fn(each)
    end
end

function Base.foreach(fn :: Function, lst :: List!{T}) :: Nothing  where T
    for each in lst
        fn(each)
    end
end

function _flatten_cons_marker(ast :: Expr, res :: Vector)
    if ast.head === :call && ast.args[1] === :^
        push!(res, ast.args[2])
        _flatten_cons_marker(ast.args[3], res)
    else
        push!(res, ast)
    end
end

function _flatten_cons_marker(ast :: Symbol, res :: Vector)
    push!(res, ast)
end

function _flatten_cons_marker(ast)
    if isa(ast, Expr)
        res = []
        _flatten_cons_marker(ast, res)
    else
        [ast]
    end
end

PatternDef.App(^) do args, guard, tag, mod

    check_ty = :($isa($tag, $Cons))

    l, r = args

    rs = _flatten_cons_marker(r)

    let l = pattern_match(l, nothing, :($tag.head), mod),
        r = map(rs) do r
                pattern_match(r, nothing, :($tag.tail), mod)
            end |>
            function (last)
                reduce((a, b) -> Expr(:&&, a, b), last)
            end
        :($check_ty && $l && $r)
    end |>

    function (last)
        if guard === nothing
            last
        else
            :($last && $guard)
        end
    end

end

# TODO:
# Base.reduce, Base.flatmap, Base.all, Base.any, Base.isempty, Base.copy
# head, tail, find, take, zip...

function isempty(lst :: Cons{T}) where T
    false
end

function isempty(lst :: Nil{T}) where T
    true
end

end


text:
# Guide.manual_color_key("Legend", ["Points", "Line"], ["green", "deepskyblue"])

DOYTICKS=true

function toggleYTicks()
  global DOYTICKS
  DOYTICKS = DOYTICKS ? false : true
  return DOYTICKS
end

function draw1D!(bd::BallTreeDensity, bins::Union{Array{Float64,1},LinSpace{Float64}}, e, c::ASCIIString="deepskyblue", myStyle::ASCIIString="";
  xlbl="X",legend=Union{})
  global DOYTICKS

  yV = evaluateDualTree(bd,bins)
  # clamp max y values
  yV[3.0 .< yV] = 3.0
  if e == Union{}
    ptArr = Any[]

    l1 = Gadfly.layer(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)))
    push!(ptArr, l1)
    push!(ptArr, Guide.xlabel(xlbl))
    if !DOYTICKS
      # e=Gadfly.plot(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)),Guide.xlabel(xlbl))
    # else
      push!(ptArr, Guide.ylabel(""))
      push!(ptArr,Guide.yticks(ticks=nothing))
      # e=Gadfly.plot(x=bins,y=yV,Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)),Guide.xlabel(xlbl),Guide.ylabel(""),Guide.yticks(ticks=nothing))
    end
    if legend != Union{}
      push!(ptArr, legend)
    end
    e = Gadfly.plot(ptArr...)
  else
    push!(e.layers, layer(x=bins, y=yV, Geom.line, Gadfly.Theme(default_color=parse(Colorant,c)))[1])
  end

  #mx = maximum(y)
  #wts = getWeights(bd)
  return e
end

function plotKDE(bd::BallTreeDensity;
          N=200, c::Array{ASCIIString,1}=["black"],rmax=-Inf,rmin=Inf,
          xlbl="X", legend=Union{})
    plotKDE([bd],N=N,c=c,rmax=rmax,rmin=rmin,xlbl=xlbl,legend=legend)
end

# function getKDERange(bd::BallTreeDensity; extend::Float64=0.1)
#   if (bd.bt.dims == 1)
#     pts = getPoints(bd)
#     rangeV = [minimum(pts),maximum(pts)]
#   else
#       return error("getKDERange(::BTD) -- multidimensional not implemented yet")
#   end
#
#   dr = extend*(rangeV[2]-rangeV[1])
#   rangeV[1] = rangeV[1] - dr;
#   rangeV[2] = rangeV[2] + dr;
#   return rangeV
# end
#
# function getKDERangeLinspace(bd::BallTreeDensity; extend::Float64=0.2, N::Int=201)
#   v = getKDERange(bd,extend=extend)
#   return linspace(v[1],v[2],N)
# end

function plotKDE(darr::Array{BallTreeDensity,1};
      N::Int=200, c::Array{ASCIIString,1}=["black"],rmax=-Inf,rmin=Inf,
      xlbl="X",legend=Union{})
    if (length(c)<2)
        c = repmat(c,length(darr))
    end
    lg = Union{}
    if legend != Union{}
      lg = Guide.manual_color_key("Legend", legend, c)
    end
    H = Union{}
    i = 0
    for bd in darr
        i+=1
        dim = 1:bd.bt.dims
        if (bd.bt.dims == 1)
          rangeV = getKDERange(bd)
          if rangeV[1] > rmin  rangeV[1] = rmin end
          if rmax > rangeV[2]  rangeV[2] = rmax end
          H=draw1D!(bd,linspace(rangeV[1],rangeV[2],N), H, c[i],xlbl=xlbl,legend=lg) #,argsPlot,argsKDE
        else
            error("plotKDE(::BTD) -- multidimensional plotting not implemented yet")
        end
    end
    return H
end

# function drawProdElement!(bd::BallTreeDensity, bins::Union{Array{Float64,1},LinSpace{Float64}}, H, offs, height, mcmc; c::ASCIIString="blue", myStyle::ASCIIString="")
#   a = 0.4
#   if offs==0.0
#     c = "red"
#     a = 1.0
#   end
#
#   ZZ = evaluateDualTree(bd,bins) .+ height
#   YY = zeros(length(ZZ)) .+ (offs - 0.0*height) .+ mcmc
#
#   H=PyPlot.plot3D(YY,bins,ZZ, color=c, alpha=a)
#
#   if offs == 0.0
#     offs -= 0.5
#   end
#   offs -= 1.0
#   return H, offs
# end
#
# function plotKDEProd!(bdArr::Array{BallTreeDensity,1},minmax; h=0.0, mcmc=0.0, N=200)
#   offs = 0.0
#   H = []
#   for bd in bdArr
#         dim = 1:bd.bt.dims
#         if (bd.bt.dims == 1)
#           pts = getPoints(bd)
#           rangeV = [minimum(pts),maximum(pts)]
#           rangeV[1] = rangeV[1] - .15*(rangeV[2]-rangeV[1]);
#           rangeV[2] = rangeV[2] + .15*(rangeV[2]-rangeV[1]);
#           bins = linspace(rangeV[1],rangeV[2],N)
#           if (offs==0)
#             minval = rangeV[1]
#             #if minval > 0.0
#             #  minval = 0
#             #end
#             bins = linspace(minval,rangeV[2],N)
#           end
#           H,offs =drawProdElement!(bd, bins, H, offs, h, mcmc) #,argsPlot,argsKDE
#         else
#             error("plotKDEProd(::BTD) -- multidimensional plotting not implemented yet")
#         end
#
#         if rangeV[1] < minmax[1]
#           minmax[1] = rangeV[1]
#         end
#         if minmax[2]<rangeV[2]
#           minmax[2] = rangeV[2]
#         end
#     end
# end

macro ifund(exp)
    local e = :($exp)
    isdefined(e.args[1]) ? :($(e.args[1])) : :($(esc(exp)))
end

function drawHorDens(pDens::Array{BallTreeDensity,1}; N::Int=200,
                    gt=Union{}, lbls=ASCIIString[], extend::Float64=0.1)
    len = length(pDens)
    h = Array{Gadfly.Plot,1}(len)
    r = Array{RemoteRef,1}(len)

    if gt!=Union{}
      for i in 1:length(pDens)
        g = kde!(gt[i][1,:],[gt[i][2,1]])
        v = getKDERange(pDens[i],extend=extend)
        p = plotKDE([pDens[i];g],c=["black";"red"],rmin=v[1],rmax=v[2],xlbl=lbls[i])
        h[i] =p;
      end
    else
      for i in 1:length(pDens)
         h[i]=plotKDE([pDens[i]],xlbl=lbls[i])
      end
    end
    #println("draw $(len) plots")
    # [r[i] = @spawn plotKDE(pDens[i], N=N) for i in 1:len]
    # [h[i] = fetch(r[i]) for i in 1:len]

    hstack(h)
end

function stackMarginals(P::Array{BallTreeDensity,1}, m::Int64)
  ret = Array{BallTreeDensity,1}(length(P))
  for i in 1:length(P)
    ret[i] =  marginal(P[i],[m])#push!()
  end
  return ret
end

function vstackedPlots(plots::Array{Gadfly.Compose.Context,1})
    evalstr = ""
    for i in 1:length(plots)
        evalstr = string(evalstr, ",plots[$(i)]")
    end
    evalstr = string("vstack(",evalstr[2:end],")")
    eval(parse(evalstr))
end

# function vArrPotentials(potens::Array{BallTreeDensity,1})
#   vv = Array{Compose.Context,1}(length(potens))
#   i = 0
#   for p in potens
#       i+=1
#       vv[i] = investigatePoseKDE(p)
#   end
#   return vv
# end


text:
Base.getindex(x::Value, inds...) = register(Base.getindex, x, inds...)

function gradient(::typeof(Base.getindex), grad, output, x::AbstractArray, inds...)
    grad_output = fill!(similar(x), 0)
    setindex!(grad_output, grad, inds...)
    (grad_output, )
end

gradient(::typeof(Base.getindex), grad, output, x::Number, ind::Int) = (grad, )

Base.selectdim(x::Value, d, i) = register(Base.selectdim, x, d, i)

function gradient(::typeof(Base.selectdim), grad, output, x::AbstractArray, d, i)
    grad_output = fill!(similar(x), 0)
    subgrad = Base.selectdim(grad_output, d, i)
    setindex!(subgrad, grad, :)
    (grad_output, )
end

Base.view(x::Value, inds...) = register(Base.view, x, inds...)

function gradient(::typeof(Base.view), grad, output, x::AbstractArray, inds...)
    grad_output = fill!(similar(x), 0)
    subgrad = view(grad_output, inds...)
    setindex!(subgrad, grad, :)
    (grad_output, )
end

Base.reshape(x::Value, dims...) = register(Base.reshape, x, dims...)
function gradient(::typeof(Base.reshape), grad, output, x, dims...)
    (reshape(grad, size(x)...), )
end


text:
module XGPaint

using Interpolations
using QuadGK
using Roots
using Cosmology
using Unitful, UnitfulAstro
using Parameters
using Random
using Healpix

import ThreadsX
import Distributions

include("./pixell.jl")
include("./model.jl")
include("./util.jl")
include("./cib.jl")
include("./radio.jl")

export get_cosmology, read_halo_catalog_hdf5
export Radio_Sehgal2009, CIB_Planck2013
export paint!, generate_sources

end # module

max_stars_repo_name: xzackli/XGPaint

max_issues_repo_name: xzackli/XGPaint

max_forks_repo_name: xzackli/XGPaint


text:
function trinary_to_decimal(str::AbstractString)
    result = 0
    try
        for (i,d) in enumerate(Iterators.reverse(str))
            p = parse(Int,d)
            p > 3 && return 0
            result += p*(3^(i-1))
        end
    catch(e)
        isa(e,ArgumentError) && return 0
        throw(e)
    end
    result
end

# Approx 2x quicker. From @edit it looks like Julia is using :ccall, which I think counts as cheating
function trinary_to_decimal_reference(str::AbstractString)
    try
        parse(Int,str,base=3)
    catch(e)
        return 0
    end
end


text:
@testset "Int → Eng" begin
    @test english(0) == "zero"
    @test english(1) == "one"
    @test english(2) == "two"
    @test english(9) == "nine"
    @test english(10) == "ten"
    @test english(11) == "eleven"
    @test english(19) == "nineteen"
    @test english(20) == "twenty"
    @test english(21) == "twenty-one"
    @test english(29) == "twenty-nine"
    @test english(30) == "thirty"
    @test english(99) == "ninety-nine"
    @test english(100) == "one hundred"
    @test english(101) == "one hundred one"
    @test english(111) == "one hundred eleven"
    @test english(199) == "one hundred ninety-nine"
    @test english(200) == "two hundred"
    @test english(999) == "nine hundred ninety-nine"
    @test english(1000) == "one thousand"
    @test english(1001) == "one thousand one"
    @test english(1100) == "one thousand one hundred"
    @test english(1999) == "one thousand nine hundred ninety-nine"
    @test english(2000) == "two thousand"
    @test english(10000) == "ten thousand"
    @test english(805097) == "eight hundred five thousand ninety-seven"
    @test english(999999) ==
        "nine hundred ninety-nine thousand nine hundred ninety-nine"
    @test english(1000000) == "one million"
    @test english(1000001) == "one million one"
    @test english(1001000) == "one million one thousand"
    @test english(1001001) == "one million one thousand one"
    @test english(800000005) == "eight hundred million five"
    @test english(1000000000) == "one billion"
    @test english(BigInt(10)^21) == "one sextillion"
    @test english(BigInt(10)^63) == "one vigintillion"
end

@testset "Int ← Eng" begin
    for T in [Int8, Int16, Int32, Int64, Int128,
              UInt8, UInt16, UInt32, UInt64, UInt128]
        @test unenglish(T, english(0))::T == 0
        @test unenglish(T, english(typemax(T)))::T == typemax(T)
    end

    for i in 0:10
        @test unenglish(Int, english(i)) == i
    end

    for i in 100:37:900
        @test unenglish(Int, english(i)) == i
    end

    for i in 3:8
        @test unenglish(Int, english(10^i)) == 10^i
        @test unenglish(Int, english(10^i + 85)) == 10^i + 85
    end

    for i in 9:63
        @test unenglish(BigInt, english(big(10)^i)) == big(10)^i
        @test unenglish(BigInt, english(big(10)^i - 85)) == big(10)^i - 85
    end
end


text:
export  MDF,
        lookahead,
        iterative_policy_evaluation,
        policy_evaluation,
        ValueFunctionPolicy,
        greedy,
        PolicyIteration,
        solve,
        backup,
        ValueIteration,
        GaussSeidelValueIteration,
        LinearProgramFormulation,
        tensorform,
        LinearQuadraticProblem

"""
    struct MDP
        γ # discount factor
        𝒮 # state space
        𝒜 # action space
        T # transition function
        R # reward function
        TR # sample transition and reward
    end
"""
struct MDP
    γ # discount factor
    𝒮 # state space
    𝒜 # action space
    T # transition function
    R # reward function
    TR # sample transition and reward
end

"""
    function lookahead(𝒫::MDP, U, s, a)

Calculate state-action value function Q(s,a) from U and valid state transition distribution T:

`Returns`
    Q(s,a) = R(s,a) + γ*sum(T(s,a,s′)*U(s′) for s′ in 𝒮)

"""
function lookahead(𝒫::MDP, U, s, a)
    𝒮, T, R, γ = 𝒫.𝒮, 𝒫.T, 𝒫.R, 𝒫.γ
    return R(s,a) + γ*sum(T(s,a,s′)*U(s′) for s′ in 𝒮)
end

"""
    function lookahead(𝒫::MDP, U::Vector, s, a)

Calculate state-action value function Q(s,a) from U and valid state transition distribution T:

    return R(s,a) + γ*sum(T(s,a,s′)*U[i] for (i,s′) in enumerate(𝒮))

Complexity: 𝒪(|𝒮|)
"""
function lookahead(𝒫::MDP, U::Vector, s, a)
    𝒮, T, R, γ = 𝒫.𝒮, 𝒫.T, 𝒫.R, 𝒫.γ
    return R(s,a) + γ*sum(T(s,a,s′)*U[i] for (i,s′) in enumerate(𝒮))
end

"""
    function iterative_policy_evaluation(𝒫::MDP, π, k_max)

 Approximate policy evaluation for policy `π`, using k_max iterations of 1-step lookahead, starting from U(s) = 0.

Complexity: 𝒪(|𝒮|² × k_max)
"""
function iterative_policy_evaluation(𝒫::MDP, π, k_max)
    𝒮, T, R, γ = 𝒫.𝒮, 𝒫.T, 𝒫.R, 𝒫.γ
    U = [0.0 for s in 𝒮]
    for k in 1:k_max
        U = [lookahead(𝒫, U, s, π(s)) for s in 𝒮]
    end
    return U
end

"""
    function policy_evaluation(𝒫::MDP, π)

Exact policy evaluation using matrix-based Moore-Penrose pseudoinverse equation.
Requires O(|𝒮|³) time.

`Returns:` [ U(s) ∀ s ∈ {𝒮} ]
"""
function policy_evaluation(𝒫::MDP, π)
    𝒮, R, T, γ = 𝒫.𝒮, 𝒫.R, 𝒫.T, 𝒫.γ

    # R(s, π(s)) ∀ s ∈ {𝒮}
    R′ = [R(s, π(s)) for s in 𝒮] # nx1 = |𝒮|x1
    T′ = [T(s, π(s), s′) for s in 𝒮, s′ in 𝒮] # nxm = |𝒮|x|𝒮|
    return (I - γ*T′)\R′
end

"""
    struct ValueFunctionPolicy
        𝒫 # problem
        U # utility function
    end
"""
struct ValueFunctionPolicy
    𝒫 # problem
    U # utility function
end

"""
    function greedy(𝒫::MDP, U, s)

Find the greedy action and its expected utility starting from state `s` and using estimated optimal (greedy) state value function `U`.  Calculates estimate as maximum over action space of 1-step lookahead:

    u, a = _findmax(a->lookahead(𝒫, U, s, a), 𝒫.𝒜)

Complexity: 𝒪(|𝒜|×|𝒮|)
"""
function greedy(𝒫::MDP, U, s)
    u, a = _findmax(a->lookahead(𝒫, U, s, a), 𝒫.𝒜)
    return (a=a, u=u)
end

"""
    function (π::ValueFunctionPolicy)(s)

Greedy policy function mapping states to greedy action.
"""
function (π::ValueFunctionPolicy)(s)
    return greedy(π.𝒫, π.U, s).a
end

# ---------------- Policy Iteration ----------------

"""
    struct PolicyIteration
        π # initial policy
        k_max # maximum number of iterations
    end
"""
struct PolicyIteration
    π # initial policy
    k_max # maximum number of iterations
end

"""
    function solve(M::PolicyIteration, 𝒫::MDP)

Iteratively improves an initial policy `π` to obtain an optimal policy for an MDP `𝒫` with discrete state and action spaces.  Returns improved policy.
"""
function solve(M::PolicyIteration, 𝒫::MDP)
    π, 𝒮 = M.π, 𝒫.𝒮
    converged = false

    for k = 1:M.k_max
        # exact solution for U_π(s)
        U = policy_evaluation(𝒫, π)
        # Use new utility function for finding greedy action.
        π′ = ValueFunctionPolicy(𝒫, U)
        # same action given for every state, => convergence (current π is best)
        if all(π(s) == π′(s) for s in 𝒮)
            break
        end
        # Not converged: update policy with new greedy policy.
        π = π′
    end
    return π
end

# -------------------------- Value Iteration --------------------------

"""
    function backup(𝒫::MDP, U, s)

Calculate maxₐ Q(s,a).
Complexity: O(|𝒮|×|𝒜|).

    return maximum(lookahead(𝒫, U, s, a) for a in 𝒫.𝒜)
"""
function backup(𝒫::MDP, U, s)
    return maximum(lookahead(𝒫, U, s, a) for a in 𝒫.𝒜)
end

"""
    struct ValueIteration
        k_max # maximum number of iterations
    end
"""
struct ValueIteration
    k_max # maximum number of iterations
end

"""
    function solve(M::ValueIteration, 𝒫::MDP)

Iteratively improves a value function `U` to obtain an optimal policy for an MDP `𝒫` with discrete state and action spaces.  Terminates after `k_max` iterations.

Order O(k_max × |𝒮| × |𝒜| × |𝒮|)
# Returns

Optimal policy `π`
"""
function solve(M::ValueIteration, 𝒫::MDP)
    U = [0.0 for s in 𝒫.𝒮]
    for k = 1:M.k_max
        U = [backup(𝒫, U, s) for s in 𝒫.𝒮]
    end
    return ValueFunctionPolicy(𝒫, U)
end

#
# @doc raw"""
#     function solve(M::ValueIteration, 𝒫::MDP, δ::Float64)
#
# Iteratively improves a value function `U` to obtain an optimal policy for an MDP `𝒫` with discrete state and action spaces.  Terminates after Bellman residual falls below threshold δ:
#
# $\|\|U\_{k+1} - U\_k\|\|\_{\infty} < \delta$
#
# # Returns
#
# Optimal policy `π`
# """
# function solve(M::ValueIteration, 𝒫::MDP, δ::Float64)
#     U = [0.0 for s in 𝒫.𝒮]
#     for k = 1:M.k_max
#         U = [backup(𝒫, U, s) for s in 𝒫.𝒮]
#     end
#     return ValueFunctionPolicy(𝒫, U)
# end

# ----------------- Asyncronous Value Iteration --------------------
#
# Similar to Value Iteration, but only subset of states updated every iteration.
#
# Common asyncronous value iteration method: Gauss-Seidel.
#   Sweeps through an ordering of states and applies Bellman update in place.
"""
    struct GaussSeidelValueIteration
        k_max # maximum number of iterations
    end
"""
struct GaussSeidelValueIteration
    k_max # maximum number of iterations
end

"""
    function solve(M::GaussSeidelValueIteration, 𝒫::MDP)

Same as Value Iteration but update U one state at a time.  Convergence speed depends on ordering of states in `𝒮`.
"""
function solve(M::GaussSeidelValueIteration, 𝒫::MDP)
    𝒮, 𝒜, T, R, γ = 𝒫.𝒮, 𝒫.𝒜, 𝒫.T, 𝒫.R, 𝒫.γ
    U = [0.0 for s in 𝒮]
    for k = 1:M.k_max
        for (i, s) in enumerate(𝒮)
            u = backup(𝒫, U, s)
            U[i] = u
        end
    end
    return ValueFunctionPolicy(𝒫, U)
end

# ------------- Linear Programming --------------

"""
    struct LinearProgramFormulation end
"""
struct LinearProgramFormulation end

"""
    function tensorform(𝒫::MDP)

Convert an MDP into its tensor form, where the states and actions consist of integer indicies, the reward function is a matrix, and the transition function is a three-dimensional tensor.
"""
function tensorform(𝒫::MDP)
    𝒮, 𝒜, R, T = 𝒫.𝒮, 𝒫.𝒜, 𝒫.R, 𝒫.T
    𝒮′ = eachindex(𝒮)
    𝒜′ = eachindex(𝒜)
    R′ = [R(s,a) for s in 𝒮, a in 𝒜]
    T′ = [T(s,a,s′) for s in 𝒮, a in 𝒜, s′ in 𝒮]
    return 𝒮′, 𝒜′, R′, T′
end

"""
    solve(𝒫::MDP) = solve(LinearProgramFormulation(), 𝒫)

Default MDP solver is the LinearProgramFormulation solver.
"""
solve(𝒫::MDP) = solve(LinearProgramFormulation(), 𝒫)

"""
    function solve(M::LinearProgramFormulation, 𝒫::MDP)
"""
function solve(M::LinearProgramFormulation, 𝒫::MDP)
    𝒮, 𝒜, R, T = tensorform(𝒫)
    model = Model(GLPK.Optimizer)
    @variable(model, U[𝒮])
    @objective(model, Min, sum(U))
    @constraint(model, [s=𝒮,a=𝒜], U[s] ≥ R[s,a] + 𝒫.γ*T[s,a,:]⋅U)
    optimize!(model)
    return ValueFunctionPolicy(𝒫, value.(U))
end

# ---------- Linear Quadratic Problem -----------

"""
    struct LinearQuadraticProblem
        Ts # transition matrix with respect to state
        Ta # transition matrix with respect to action
        Rs # reward matrix with respect to state (negative semidefinite)
        Ra # reward matrix with respect to action (negative definite)
        h_max # horizon
    end
"""
struct LinearQuadraticProblem
    Ts # transition matrix with respect to state
    Ta # transition matrix with respect to action
    Rs # reward matrix with respect to state (negative semidefinite)
    Ra # reward matrix with respect to action (negative definite)
    h_max # horizon
end

"""
    function solve(𝒫::LinearQuadraticProblem)
"""
function solve(𝒫::LinearQuadraticProblem)
    Ts, Ta, Rs, Ra, h_max = 𝒫.Ts, 𝒫.Ta, 𝒫.Rs, 𝒫.Ra, 𝒫.h_max
    V = zeros(size(Rs))
    πs = Any[s -> zeros(size(Ta, 2))]
    for h in 2:h_max
        V = Ts'*(V - V*Ta*((Ta'*V*Ta + Ra) \ Ta'*V))*Ts + Rs
        L = -(Ta'*V*Ta + Ra) \ Ta' * V * Ts
        push!(πs, s -> L*s)
    end
    return πs
end


text:
@testset "Jordan and Rational canonical form" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^6 + rand(Qx, 1:5, -3:3)
  M = companion_matrix(f)
  @test minpoly(Qx, M) == f
  @test charpoly(Qx, M) == f

  ff = factor(f)
  g = first(keys(ff.fac))
  g = divexact(g, leading_coefficient(g))
  J = Hecke.jordan_block(g, 2)
  @test minpoly(Qx, J) == g^2
  @test charpoly(Qx, J) == g^2

  M = identity_matrix(FlintQQ, 6)
  for i = 1:6
    for j = i+1:6
      M[i, j] = rand(FlintQQ, -2:2)
    end
  end
  M1 = transpose(M)
  @test Hecke.issimilar(M, M1)
  MT = Hecke.conjugating_matrix(M, M1)
  @test MT*M1*inv(MT) == M
  J, S = jordan_normal_form(M)
  @test S*M*inv(S) == J
  for i = 1:6
    for j = 1:i-1
      @test iszero(J[i, j])
    end
  end
  for i = 1:5
    @test isone(J[i, i+1]) || iszero(J[i, i+1])
  end
  for i = 1:6
    for j = i+2:6
      @test iszero(J[i, i+2])
    end
  end

  f = x^3+ rand(Qx, 1:2, -2:2)
  g = x^3 + rand(Qx, 1:2, -2:2)
  f1 = f
  f2 = f*g
  C = zero_matrix(FlintQQ, 9, 9)
  Hecke._copy_matrix_into_matrix(C, 1, 1, companion_matrix(f1))
  Hecke._copy_matrix_into_matrix(C, 4, 4, companion_matrix(f2))
  CF, TM = rational_canonical_form(C)
  @test CF == C
  @test TM * C * inv(TM) == CF

end


text:
struct Identity{M,N,S} <: TaskMap{M,N,S} end
struct IdentityT{M,N,S} <: TaskMapT{M,N,S}
    base_map::Identity{M,N,S}
end
IdentityT{M,N,S}() where {M,N,S} = IdentityT{M,N,S}(Identity{M,N,S}())

task_map_emb(::EmbRep, ::EmbRep, xm, task_map::Identity{ℝ{m},ℝ{m}}) where m = xm

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::Identity{𝕊{m},𝕊{m}}) where m = xme
domain_coord_rep(::Identity{𝕊{n},𝕊{n}}) where n = EmbRep()
codomain_coord_rep(::Identity{𝕊{n},𝕊{n}}) where n = EmbRep()


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module Common

export Vec, Vec2, IVec, IVec2, unsafe_bitflip!, flipbits!, discretize,
       LocalFields, AllButOne

import Base: iterate

const Vec  = Vector{Float64}
const Vec2 = Vector{Vec}
const IVec = Vector{Int}
const IVec2 = Vector{IVec}

@inline function unsafe_bitflip!(Bc::Array{UInt64}, i::Int)
    i1, i2 = Base.get_chunks_id(i)
    u = UInt64(1) << i2
    @inbounds begin
        #Bc[i1] ⊻= u
        Bc[i1] = Bc[i1] ⊻ u
    end
end
@inline unsafe_bitflip!(B::BitArray, i::Int) = unsafe_bitflip!(B.chunks, i)

flipbits!(B::BitArray) = B .= .!(B)

mutable struct LocalFields{ET}
    lfields::Vector{ET}
    lfields_last::Vector{ET}
    move_last::Int
    function LocalFields{ET}(N::Int) where {ET}
        lfields = zeros(ET, N)
        lfields_last = zeros(ET, N)
        return new{ET}(lfields, lfields_last, 0)
    end
end

function discretize(x::T, LEV::Tuple{S,Vararg{S}}) where {T<:Real,S<:Real}
    d = LEV[1]
    r = x - d
    for l = 2:length(LEV)
        d1 = LEV[l]
        r1 = x - d1
        if abs(r1) < abs(r)
            d, r = d1, r1
        end
    end
    return d, r
end

function discretize(cvec::Vector{T}, LEV::Tuple{S,Vararg{S}}) where {T<:Real,S<:Real}
    N = length(cvec)
    fields = Array{S}(undef, N)
    rfields = Array{T}(undef, N)
    for (i,x) in enumerate(cvec)
        d, r = discretize(x, LEV)
        fields[i] = d
        rfields[i] = r
    end
    return fields, rfields
end

function discretize(cvec::NTuple{N,T}, LEV::Tuple{S,Vararg{S}}) where {N,T<:Real,S<:Real}
    fields = Array{S}(undef, N)
    rfields = Array{T}(undef, N)
    for (i,x) in enumerate(cvec)
        d, r = discretize(x, LEV)
        fields[i] = d
        rfields[i] = r
    end
    return tuple(fields...), tuple(rfields...)
end

# iterate over 1:N except i, i.e.
#   1, 2, 3, ..., i-1, i+1, ..., N-1, N
# useful for fully-connected models
mutable struct AllButOne
    N::Int
    i::Int
    function AllButOne(N::Integer, i::Integer)
        N ≥ 1 || throw(ArgumentError("N must me larger than 1, given: $N"))
        1 ≤ i ≤ N || throw(ArgumentError("expected 1 ≤ i ≤ N, given i=$i N=$N"))
        return new(N, i)
    end
end

function iterate(n::AllButOne, j = 1 + (n.i == 1))
    j == n.N + 1 && return nothing
    return j, j + 1 + (j == n.i - 1)
end
Base.length(n::AllButOne) = n.N - 1

end # module


text:
_savehtml(pl::py.SyncPlot,filename,js) = py.savefig(pl,filename,js=js)
function _saverendered(pl::py.SyncPlot,filename,delay)
    display(pl)
    sleep(delay) # the problem is we don't know how long
    py.savefig(pl, filename)
    close(pl)
end

function saveplot(pl::py.SyncPlot,format::Symbol,outPath::AbstractString; delay=5, js=:remote)
    filename = string(outPath,'.',format)
    if format==:html
        _savehtml(pl,filename,js)
    elseif format in [:png, :pdf, :svg]
        _saverendered(pl,filename,delay)
    else
        error("Unknown image format $format")
    end
    filename
end

saveplot(pl::py.SyncPlot,formats::AbstractArray{Symbol},outPath::AbstractString;kwargs...) = map(fmt->saveplot(pl,fmt,outPath;kwargs...), formats)


text:
module Loader

    function test()
        println("test!")
    end

    export test
end


text:
@add_benchmark function LT_constraints(handle_problem)
    p = satisfy()
    x = [Variable() for _ = 1:1000]
    for (i, xi) in enumerate(x)
        push!(p.constraints, xi <= 1.0 * i)
    end
    return handle_problem(p)
end

@add_benchmark function LT_constraint(handle_problem)
    p = satisfy()
    x = Variable(1000)
    push!(p.constraints, x <= collect(1.0:1000.0))
    return handle_problem(p)
end

@add_benchmark function GT_constraints(handle_problem)
    p = satisfy()
    x = [Variable() for _ = 1:1000]
    for (i, xi) in enumerate(x)
        push!(p.constraints, xi >= 1.0 * i)
    end
    return handle_problem(p)
end

@add_benchmark function GT_constraint(handle_problem)
    p = satisfy()
    x = Variable(1000)
    push!(p.constraints, x >= collect(1.0:1000.0))
    return handle_problem(p)
end

@add_benchmark function equality_constraints(handle_problem)
    p = satisfy()
    x = [Variable() for _ = 1:1000]
    for (i, xi) in enumerate(x)
        push!(p.constraints, xi == 1.0 * i)
    end
    return handle_problem(p)
end

@add_benchmark function equality_constraint(handle_problem)
    p = satisfy()
    x = Variable(1000)
    push!(p.constraints, x == collect(1.0:1000.0))
    return handle_problem(p)
end

@add_benchmark function SDP_constraint(handle_problem)
    p = satisfy()
    x = Variable(44, 44) # 990 vectorized entries
    push!(p.constraints, x ⪰ 0)
    return handle_problem(p)
end

@add_benchmark function SDP_constraints(handle_problem)
    p = satisfy()
    x = [Variable(4, 4) for _ = 1:100] # 1000 total vectorized entries
    for v in x
        push!(p.constraints, v ⪰ 0)
    end
    return handle_problem(p)
end


text:
using Test
using Cleaner: CleanTable, drop_missing, size, cols

@testset "Drop_missing functions working as expected" begin
    testCT = CleanTable([:A, :B, :C], [[1, 2, 3], [4, missing, 5], ["x", "y", ""]])

    @test size(drop_missing(testCT)) == (2, 3)
    @test size(drop_missing(testCT; missing_values=[""])) == (1, 3)

    if Threads.nthreads() > 1
        testRM = CleanTable([:A, :B],
        [append!(Union{Missing, Int64}[missing], collect(1:1_000_000)),
        append!(Union{Missing, Int64}[missing], collect(1:1_000_000))])

        testRM1 = drop_missing(testRM; missing_values=[1, 3, 10])

        @test size(testRM1) == (999_997, 2)
        @test !in(10, cols(testRM1)[1]) && !in(10, cols(testRM1)[2])
    end
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
#=
computeKernels:
- Julia version:
- Author: qd4314
- Date: 2019-03-21
=#

function computeScatteringKernelConvolution(KS::KineticSolver)
    Nq = KS.nquadpoints
    weightMatrix = zeros(Nq, Nq)
    for q = 1:Nq
        Omega = KS.Q.pointsxyz[q, :]
        for l = 1:Nq
            Omegaprime = KS.Q.pointsxyz[l, :]
            #eps = KS.convolutionwidth/Nq^2 # do we need a second parameter to adjust how isotropic the scattering is?
            eps = KS.convolutionwidth / Nq
            #x = getdist(Omega,Omegaprime,"sphere")
            x = (1 - dot(Omega, Omegaprime))
            val = (KS.Q.weights[l] * 1.0 / eps * exp(-x^2 / (eps^2)))
            if (val < 1e-8)
                weightMatrix[q, l] = 0.0
            else
                weightMatrix[q, l] = val
            end
        end
        #weightMatrix[q, q] -= sum(weightMatrix[q, :])
        #weightMatrix[q,:] = weightMatrix[q,:]./KS.Q.weights[:]
        weightMatrix[q, :] = weightMatrix[q, :] ./ (sum(weightMatrix[q, :])) # normalize
        weightMatrix[q, :] = weightMatrix[q, :] ./ KS.Q.weights[:] # ensure that weightMatrix*KS.Q.weights[:] = ones(nq)
    end
    weightMatrix = sparse(weightMatrix)

    KS.OutscatteringVector = zeros(KS.nquadpoints)
    for k = 1:KS.nquadpoints
        KS.OutscatteringVector[k] = sum(weightMatrix[k, :])
    end

    return weightMatrix
end

function numericnnz(A, treshold = 1e-16)
    # Returns the number of nonzero entries of 'A'
    # A nonzero entry is any entry that has an
    # absolute value above the specified 'treshold'.
    count = 0
    for i = 1:prod(size(A))
        count += (abs(A[i]) < treshold)
    end
    return count
end

function analyze_nnz_for_convwidth(KS)
    N = zeros(0)
    I = zeros(0)
    nelements = 0
    for i = 10:10:1000
        KS.convolutionwidth = i
        S = computeScatteringKernelConvolution(KS)
        nelements = prod(size(S))
        n = nnz(S, 1e-8)
        append!(N, n)
        append!(I, i)
        println("$i, $n")
    end
    lineplot(I, N / nelements, xlabel = "convolutionwidth", ylabel = "ratio of nnz")
    return I, N
end

function computeScatteringKernelHenyeyGreenstein(KS::KineticSolver, g = 0.00000001)
    Nq = KS.nquadpoints
    HGScattering = zeros(Nq, Nq)
    hg(g, mu) = 1 / 4 / pi * (1 - g^2) ./ (1 + g^2 - 2 * g * mu) .^ (3 / 2)
    for q = 1:Nq
        for l = 1:Nq
            Omega = KS.Q.pointsxyz[q, :]
            Omegaprime = KS.Q.pointsxyz[l, :]
            mu = dot(Omega, Omegaprime)

            val = hg(g, mu)
            HGScattering[q, l] = ifelse(val < 1e-12, 0.0, val) * KS.Q.weights[l]
        end
        #	HGScattering[q, :] /= sum(HGScattering[q, :])##
        #	HGScattering[q,:] .* = KS.Q.weights[:] / sum(KS.Q.weights[:])
        #HGScattering[q,q] -= 1
    end
    #print((HGScattering*ones(Nq))
    v = HGScattering * ones(Nq)
    for q = 1:Nq
        for qprime = 1:Nq
            HGScattering[q, qprime] /= v[q]
        end
    end

    return HGScattering
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using SchumakerSpline
tol = 10*eps()

x = collect(range(1, stop=6, length=1000))
y = log.(x) + sqrt.(x)
analytical_first_derivative(e) = 1/e + 0.5 * e^(-0.5)
spline = Schumaker(x,y)
for i in 1:length(x)
    abs(evaluate(spline, x[i]) - y[2]) < tol
end

# Testing First derivative.
first_derivatives = evaluate.(spline, x, 1)
maximum(abs.(first_derivatives .- analytical_first_derivative.(x))) < 0.002

# Testing second derivatives
second_derivatives = evaluate.(spline, x, 2)
maximum(second_derivatives) < tol

# Test higher derivative
higher_derivatives = evaluate.(spline, x, 3)
maximum(second_derivatives) < tol

# Testing Integrals
analytic_integral(lhs,rhs) = rhs*log(rhs) - rhs + (2/3) * rhs^(3/2) - ( lhs*log(lhs) - lhs + (2/3) * lhs^(3/2) )
lhs = 2.0
rhs = 2.5
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.01

lhs = 1.2
rhs = 4.3
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.01

lhs = 0.8
rhs = 4.0
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.03

# Testing creation of a spline with gradient information.
first_derivs = analytical_first_derivative.(x)
spline = Schumaker(x,y; gradients = first_derivs)
first_derivatives = evaluate.(spline, x, 1)
maximum(abs.(first_derivatives .- analytical_first_derivative.(x))) < tol
# Testing creation of a spline with only the gradients on the edges.
first_derivs = analytical_first_derivative.(x)
spline = Schumaker(x,y; left_gradient = first_derivs[1], right_gradient = first_derivs[length(first_derivs)])
first_derivatives = evaluate.(spline, x, 1)
gaps = abs.(first_derivatives .- analytical_first_derivative.(x))
gaps[1] < tol
gaps[length(gaps)] < tol
minimum(gaps[2:(length(gaps)-1)]) > 10* tol

# Testing the other syntax for evaluation.
abs(spline(1.4) - evaluate(spline, 1.4)) < eps()
abs(spline(1.5) - evaluate(spline, 1.5)) < eps()
abs(spline(1.6) - evaluate(spline, 1.6)) < eps()

#=
# should be two random roots and one optima.
x = collect(range(-10, stop=10, length=1000))
function random_function(a::Int)
    vertex = (mod(a * 97, 89)- 45)/5
    y = -(x .- vertex).^2 .+ 1
    sp = Schumaker(x,y)
    return sp, vertex
end
using Optim
sp, vertex = random_function(2)
@time optimafinder = find_optima(sp)
@time optimize(x -> evaluate(sp,x[1]), -5.0, 5.0 )
=#

#= Testing AAD. Not in the full test batch to avoid another dependency
x = collect(range(-10, stop=10, length=1000))
function random_function(a::Int)
    vertex = (mod(a * 97, 89)- 45)/5
    y = -(x .- vertex).^2 .+ 1
    sp = Schumaker(x,y)
    return sp, vertex
end
sp, vertex = random_function(2)
function spl(x::Array{<:Real,1})
    return sum(evaluate.(Ref(sp), x))
end

using ForwardDiff
ForwardDiff.gradient(spl, [2,3])
=#


text:
# Ball Tree

NO_CHILD = -1 # Maybe do int64(uint32(-1)) ??
FIELD_NAMES = ["D", "N", "centers", "ranges", "weights",
            "lower", "upper", "leftch", "rightch", "perm"]
nfields = 10;

type BallTree
  dims::Int64                     # dimension of data
  num_points::Int64               # of points
  centers::Array{Float64,1}       # ball centers, dims numbers per ball
  ranges::Array{Float64,1}        # bounding box ranges, dims per ball, dist from center to one side
  weights::Array{Float64,1}       # total weight in each ball

  left_child::Array{Int64,1}
  right_child::Array{Int64,1}     # left, right children; no parent indices
  lowest_leaf::Array{Int64,1}
  highest_leaf::Array{Int64,1}    # lower & upper leaf indices for each ball
  permutation::Array{Int64,1}     # point's position in the original data

  next::Int64                     # internal var for placing the non-leaf nodes

  swapHandle::Function
  calcStatsHandle::Function
  data
end

root() = 1
Ndim(bt::BallTree) = bt.dims
Npts(bt::BallTree) = bt.num_points
Npts(bt::BallTree, i::Int64) = bt.highest_leaf[i]-bt.lowest_leaf[i]+1

## todo make sure these two are working properly -- this is the problem, we are always assigning new memory for each return
#function center(bt::BallTree, i::Int64)
#    return bt.centers[((i-1)*bt.dims+1):end]
#end
center(bt::BallTree, i::Int64) = bt.centers[((i-1)*bt.dims+1):end]
center(bt::BallTree, i::Int64, k::Int) = bt.centers[((i-1)*bt.dims+k)]

#function rangeB(bt::BallTree, i::Int64)
#    return bt.ranges[((i-1)*bt.dims+1):end]
#end
rangeB(bt::BallTree, i::Int64) = bt.ranges[((i-1)*bt.dims+1):end]
rangeB(bt::BallTree, i::Int64, k::Int) = bt.ranges[((i-1)*bt.dims+k)]

#function weight(bt::BallTree, i::Int64)
#    return bt.weights[i]
#end
weight(bt::BallTree, i::Int64) = bt.weights[i]

#function isLeaf(bt::BallTree, ind::Int64)
#    return ind >= bt.num_points
#end
isLeaf(bt::BallTree, ind::Int64) = ind >= bt.num_points

validIndex(bt::BallTree, ind::Int64) = ((0<ind) && (ind <= 2*bt.num_points))

left(bt::BallTree, i::Int64) = bt.left_child[i]

right(bt::BallTree, i::Int64) = bt.right_child[i]

leafFirst(bt::BallTree, i::Int64) = bt.lowest_leaf[i]

leafLast(bt::BallTree, i::Int64) = bt.highest_leaf[i]

# Convert a BallTree::index to the numeric index in the original data
getIndexOf(bt::BallTree, i::Int64) = bt.permutation[i]

function swap!(data, _i::Int64, _j::Int64)
  return data.swapHandle(data, _i, _j)
end
function calcStats!(data, root::Int64)
  #@show "Fancy calcStats"
  return data.calcStatsHandle(data, root)
end

# Swap the ith leaf with the jth leaf.  Actually, only swap the
# weights, permutation, and centers, so only for swapping
# leaves. Will not swap ranges correctly and will not swap children
# correctly.
function swapBall!(bt::BallTree, _i::Int64, _j::Int64)
  #println("swapBall! -- (i,j)=$((_i,_j))")
  i = _i
  j = _j
  tmp = 0.
  if (i==j)
    return Union{}
  end

  # swap weights
  tmp = bt.weights[i]
  bt.weights[i] = bt.weights[j]
  bt.weights[j] = tmp

  # swap perm
  k = bt.permutation[i]
  bt.permutation[i] = bt.permutation[j]
  bt.permutation[j] = k;

  # swap centers
  i = (i-1)*bt.dims
  j = (j-1)*bt.dims;
  for k in 1:bt.dims
    i+=1
    j+=1
    tmp = bt.centers[i];
    bt.centers[i]  = bt.centers[j];
    bt.centers[j] = tmp;
  end
end

# Find the dimension along which the leaves between low and high
# inclusive have the greatest variance
function  most_spread_coord(bt::BallTree, low::Int64, high::Int64)
  #BallTree::index dimension, point, max_dim;
  #double mean, variance, max_variance;
  #println("most_spread_coord -- low, high = $((low, high))")
  max_variance = 0
  max_dim = 1

  for dimension = 1:bt.dims
    mean = 0
    for point = (bt.dims*(low-1) + dimension):bt.dims:(bt.dims*(high-1))
      mean += bt.centers[point]
    end
    mean /= (high - low)

    variance = 0
    for point in (bt.dims*(low-1) + dimension):bt.dims:(bt.dims*(high-1))
      variance += (bt.centers[point] - mean) * (bt.centers[point] - mean);
    end
    if (variance > max_variance)
      max_variance = variance;
      max_dim = dimension;
    end
  end

  #println("most_spread_coord -- max_var=$(max_variance), return=$(max_dim)")
  return max_dim;
end

# straight from CLR, the unrandomized partition algorithm for
# quicksort.  Partitions the leaves from low to high inclusive around
# a random pivot in the given dimension.  Does not affect non-leaf
# nodes, but does relabel the leaves from low to high.
function partition!(bt::BallTree, dimension::Int64, low::Int64, high::Int64)
  pivot = low;  # not randomized, could set pivot to a random element

  while (low < high)
    while(bt.centers[bt.dims*(high-1) + dimension] >= bt.centers[bt.dims*(pivot-1) + dimension])
      high-=1
    end
    while(bt.centers[bt.dims*(low-1) + dimension] < bt.centers[bt.dims*(pivot-1) + dimension])
      low+=1
    end

    bt.swapHandle(bt, low, high)
    pivot = high
  end

  return high;
end

# Function to partition the data into two (equal-sized or near as possible)
#   sets, one of which is uniformly greater than the other in the given
#   dimension.
function select!(bt::BallTree, dimension::Int64, position::Int64, low::Int64, high::Int64)
    m = 0
    r = 0
    i = 0
  while low < high
    r = (floor(Int64,(low + high)/2))
    swap!(bt.data, r, low)
    m = low;
    for i in (low):high
      if (bt.centers[dimension+bt.dims*(i-1)] < bt.centers[dimension+bt.dims*(low-1)])
        m+=1
        swap!(bt.data, m, i);
      end
    end
    swap!(bt.data, low, m);
    if (m <= position) low=m+1 end
    if (m >= position) high=m-1 end
  end
  return Union{}
end

# Calculate the statistics of level "root" based on the statistics of
#   its left and right children.
function calcStatsBall!(bt::BallTree, root::Int64)
  #println("calcStatsBall! -- root=$(root)")
  Ni = 0
  NiL = 0
  NiR = 0
  d = 0

  leftI = left(bt, root)
  rightI=right(bt, root);   # get children indices
  if (!(validIndex(bt, leftI)) || !(validIndex(bt, rightI)))
    return Union{} # nothing to do if this
  end           #   isn't a parent node

  # figure out the center and ranges of this ball based on it's children
  maxi = 0.
  mini = 0.
  for d=1:bt.dims
    a = center(bt, leftI, d) + rangeB(bt, leftI, d)
    b = center(bt, rightI, d) + rangeB(bt, rightI, d)
    #@show (d, leftI, rightI, a, b)
    if (a > b)
      maxi = center(bt, leftI, d) + rangeB(bt, leftI, d)
    else
      maxi = center(bt, rightI, d) + rangeB(bt, rightI, d)
    end

    if (center(bt, leftI, d) - rangeB(bt, leftI, d) < center(bt, rightI, d) - rangeB(bt, rightI, d))
      mini = center(bt, leftI, d) - rangeB(bt, leftI, d)
    else
      mini = center(bt, rightI, d) - rangeB(bt, rightI, d)
    end

    #@show (d, mini,maxi);
    bt.centers[(root-1)*bt.dims+d] = (maxi+mini) / 2.0;
    bt.ranges[(root-1)*bt.dims+d] = (maxi-mini) / 2.0;
  end

  # if the left ball is the same as the right ball (should only
  # happen when calling the function directly with the same argument
  # twice), don't count the weight twice
  if (leftI != rightI)
    bt.weights[root] = bt.weights[leftI] + bt.weights[rightI];
  else
    bt.weights[root] = bt.weights[leftI]
  end
  return Union{}
end

# Given the leaves, build the rest of the tree from the top down.
# Split the leaves along the most spread coordinate, build two balls
# out of those, and then build a ball around those two children.
function buildBall!(bt::BallTree, low::Int64, high::Int64, root::Int64)
  #println("buildBall! -- (low, high, root)=$((low, high, root))")
  # special case for N=1 trees
  if (low == high)
    bt.lowest_leaf[root] = low;
    bt.highest_leaf[root] = high;
    bt.left_child[root] = low;

    # point right child to the same as left for calc stats, and then
    # point it to the correct NO_CHILD afterwards.  kinda kludgey
    bt.right_child[root] = high;
    calcStats!(bt.data, root)
    bt.right_child[root] = NO_CHILD;
    return Union{}
  end

  #BallTree::index coord, split, left, right;
  coord = most_spread_coord(bt, low, high); # find dimension of widest spread

  # split the current leaves into two groups, to build balls on them.
  # Choose the most spread coordinate to split them on, and make sure
  # there are the same number of points in each (+-1 for round off
  # error).
  split = (floor(Int64,(low + high) / 2))
  #@show coord, split, low, high
  select!(bt, coord, split, low, high)

  # an alternative is to use partition, but that doesn't deal well
  # with repeated numbers and it doesn't split into balanced sets.
  # split = partition(coord, low, high);

  # if the left sub-tree is just one leaf, don't make a new non-leaf
  # node for it, just point left_idx directly to the leaf itself.
  if (split <= low)
    left = low
  else
    left = bt.next;
    bt.next+=1
  end

  # same for the right
  if (split+1 >= high)
    right = high
  else
    right = bt.next
    bt.next+=1
  end

  bt.lowest_leaf[root]  = low;
  bt.highest_leaf[root] = high;
  bt.left_child[root]   = left;
  bt.right_child[root]  = right;

  # build sub-trees if necessary
  if(left != low)
    buildBall!(bt, low, split, left)
  end
  if(right != high)
    buildBall!(bt, split+1, high, right)
  end

  calcStats!(bt.data, root);
  return Union{}
end

# Public method to build the tree, just calls the private method with
# the proper starting arguments.
function buildTree!(bt::BallTree)
  #println("buildTree!(::BallTree) -- is running")
  i=bt.num_points
  for j in 1:bt.num_points
    for k in 1:bt.dims
      bt.ranges[i*bt.dims+k] = 0
    end
    i+=1
    bt.lowest_leaf[i] = i
    bt.highest_leaf[i] = i
    bt.left_child[i] = i
    bt.right_child[i] = NO_CHILD #int64(uint32(NO_CHILD))
    bt.permutation[i] = j
  end
  bt.next = 2

  buildBall!(bt, bt.num_points+1, 2*bt.num_points, 1); # chgd for indexing 1
  return Union{}
end

# New BallTree
function makeBallTree(_pointsMatrix::Array{Float64,2}, _weights::Array{Float64,1}, suppressBuildTree=false)
  # get fields from input arguments
  Nd = size(_pointsMatrix,1);
  Np = size(_pointsMatrix,2);

  _points = reshape(_pointsMatrix, Nd*Np, 1)
  centers = zeros(Nd*2*Np)
  ranges = zeros(Nd*2*Np)
  weights = zeros(2*Np)
  centers[(Np*Nd+1):end] = _points[1:end]
  weights[(Np+1):end] = _weights[1:end]

  bt = BallTree(Nd, Np, centers, ranges, weights,
                  ones(Int64,2*Np), ones(Int64,2*Np), ones(Int64,2*Np),
                  ones(Int64,2*Np), zeros(Int64,2*Np),
                  0, swapBall!, calcStatsBall!, [])
  bt.data = bt

  if (Np > 0 && suppressBuildTree!=true)
    buildTree!(bt)
  end
  return bt
end

function printBallTree(bt::BallTree)
    @show round(bt.centers,15);
    @show round(bt.weights,15);
    @show round(bt.ranges,15);
    @show bt.left_child;
    @show bt.right_child;
    @show bt.highest_leaf;
    @show bt.lowest_leaf;
    @show bt.permutation;
    return Union{}
end

function test01()
  ## test point 01
  ## matlab results
  #centers=[30]=4.4, 3.6, 5, 3.2, 0.18, 2, 5.5, 6.5, 7.5, 2.6, 0.18, 2.5, 0, 0, 0, 1.9, 1.2, 3, 3.3, -0.81, 2, 4.6, -0.56, 1, 4, 5, 6, 7, 8, 9,
  #weights=[10]=1, 0.6, 0.4, 0.4, 0, 0.2, 0.2, 0.2, 0.2, 0.2,
  #ranges=[30]=2.6, 4.4, 4, 1.4, 0.98, 1, 1.5, 1.5, 1.5, 0.7, 0.98, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  #highest_leaf=[10]=9, 7, 9, 6, 0, 5, 6, 7, 8, 9,
  #lowest_leaf=[10]=5, 5, 8, 5, 0, 5, 6, 7, 8, 9,
  #permutation=[10]=0, 0, 0, 0, 0, 2, 1, 0, 3, 4,
  mus = [[4.6173,    3.2641,    1.8729, 4, 7 ]',
   [-0.5592,   -0.8088,    1.1610, 5, 8]',
   [1.,2.,3, 6, 9]']

  bt = makeBallTree(mus, 0.2*ones(5))
  printBallTree(bt);
end


text:
function tune!(
    krg::Kriging,
    kernbounds,
    meanbounds,
    opt::Optimizer...;
    verbose::Bool = false,
)
    ngp = length(krg.gps)
    for o in opt
        for i in 1:ngp
            fx, x, ret = optimize!(krg.gps[i], kernbounds, meanbounds, o)
            if verbose
                @info "gp[$i] optimized with $(o.algo), fx = $fx"
            end
        end
    end
end

function tune!(
    krg::Kriging,
    kernbounds,
    meanbounds;
    global_eval::Integer = 0,
    global_population::Integer = 0,
    local_eval::Integer = 1000,
    verbose::Bool = false,
)
    if global_eval > 0
        tune!(
            krg,
            kernbounds,
            meanbounds,
            ISRESOptimizer(
                population = global_population,
                maxeval = global_eval,
            ),
            LBFGSOptimizer(maxeval = local_eval),
            verbose = verbose,
        )
    else
        tune!(
            krg,
            kernbounds,
            meanbounds,
            LBFGSOptimizer(maxeval = local_eval),
            verbose = verbose,
        )
    end
end

function optimize!(gp::GP.GPE, kernbounds, meanbounds, opt::Optimizer)
    params_kwargs = GP.get_params_kwargs(
        gp;
        domean = true,
        kern = true,
        noise = false,
        lik = false,
    )
    fobj = if need_gradient(opt)
        (hyp::AbstractVector, grad::AbstractVector) -> begin
            prev = GP.get_params(gp; params_kwargs...)
            try
                GP.set_params!(gp, hyp; params_kwargs...)
                GP.update_target_and_dtarget!(gp; params_kwargs...)
                grad[:] = gp.dtarget
                return gp.target
            catch err
                # reset parameters to remove any NaNs
                GP.set_params!(gp, prev; params_kwargs...)
                if !all(isfinite.(hyp))
                    println(err)
                    return -Inf
                elseif isa(err, ArgumentError)
                    println(err)
                    return -Inf
                elseif isa(err, LinearAlgebra.PosDefException)
                    println(err)
                    return -Inf
                else
                    throw(err)
                end
            end
        end
    else
        (hyp::AbstractVector, grad::AbstractVector) -> begin
            prev = GP.get_params(gp; params_kwargs...)
            try
                GP.set_params!(gp, hyp; params_kwargs...)
                GP.update_target!(gp)
                return gp.target
            catch err
                # reset parameters to remove any NaNs
                GP.set_params!(gp, prev; params_kwargs...)

                if !all(isfinite.(hyp))
                    println(err)
                    return -Inf
                elseif isa(err, ArgumentError)
                    println(err)
                    return -Inf
                elseif isa(err, LinearAlgebra.PosDefException)
                    println(err)
                    return -Inf
                else
                    throw(err)
                end
            end
        end
    end
    init = GP.get_params(gp; params_kwargs...)
    lb, ub = GP.bounds(
        gp,
        nothing,
        meanbounds,
        kernbounds,
        nothing;
        domean = true,
        kern = true,
        noise = false,
        lik = false,
    )
    fx, x, ret = maximize(opt, fobj, lb, ub, init)
    GP.set_params!(gp, x; params_kwargs...)
    GP.update_target_and_dtarget!(gp; params_kwargs...)
    fx, x, ret
end


text:
using Dash, DashHtmlComponents, DashCoreComponents, DataFrames

df8 = DataFrame(a = [1, 2, 3],
               b = [4, 1, 4],
               c = ["x", "y", "z"])
s
app = dash()

app.layout = html_div() do
    dcc_dropdown(id="dropdown-7",
                 options = [
                    (label = i, value = i) for i in df8[:, "c"]
                 ],
                 value="x"
    ),
    html_div(id="output-7")
end

callback!(
    app,
    Output("output-7", "children"),
    Input("dropdown-7", "value"),
) do value
    df8f = df8[df8.c .== value, :]
    return length(df8f[1, :])
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
"""

    conv4(w, x; kwargs...)

Execute convolutions or cross-correlations using filters specified
with `w` over tensor `x`.

Currently KnetArray{Float32/64,4/5} and Array{Float32/64,4} are
supported as `w` and `x`.  If `w` has dimensions `(W1,W2,...,I,O)` and
`x` has dimensions `(X1,X2,...,I,N)`, the result `y` will have
dimensions `(Y1,Y2,...,O,N)` where

    Yi=1+floor((Xi+2*padding[i]-Wi)/stride[i])

Here `I` is the number of input channels, `O` is the number of output
channels, `N` is the number of instances, and `Wi,Xi,Yi` are spatial
dimensions.  `padding` and `stride` are keyword arguments that can be
specified as a single number (in which case they apply to all
dimensions), or an array/tuple with entries for each spatial
dimension.

# Keywords

* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `upscale=1`: upscale factor for each dimension.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `algo=0`: specifies which convolution algorithm shoud be used to compute the results. See the CUDNN User Guide for details.
* `workSpace=C_NULL`: data pointer to GPU memory to a workspace needed to able to execute the specified algorithm.
* `workSpaceSizeInBytes=0`: the size in bytes of the provided workSpace. Default=0.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.

"""
function conv4{T}(w::KnetArray{T},x::KnetArray{T};
                  handle=cudnnhandle(), algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, alpha=1,
                  o...) # padding=0, stride=1, upscale=1, mode=0
    y = similar(x, cdims(w,x;o...))
    beta=0 # nonzero beta does not make sense when we create y
    @cuda(cudnn, cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,   UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(T(alpha)),TD(x),x,FD(w),w,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(T(beta)),TD(y),y)
    return y
end

function conv4x{T}(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T};
                   handle=cudnnhandle(), alpha=1, algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0,
                   o...) # padding=0, stride=1, upscale=1, mode=0
    beta = 0
    dx = similar(x)
    if cudnnVersion >= 4000
        @cuda(cudnn,cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(T(beta)),TD(dx),dx)
    elseif cudnnVersion >= 3000
        @cuda(cudnn,cudnnConvolutionBackwardData_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(T(beta)),TD(dx),dx)
    else
        @cuda(cudnn,cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,       Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),Ref(T(beta)),TD(dx),dx)
    end
    return dx
end

function conv4w{T}(w::KnetArray{T},x::KnetArray{T},dy::KnetArray{T};
                   handle=cudnnhandle(), alpha=1, algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0,
                   o...) # padding=0, stride=1, upscale=1, mode=0
    beta = 0
    dw = similar(w)
    if cudnnVersion >= 4000
        @cuda(cudnn,cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(T(beta)),FD(dw),dw)
    elseif cudnnVersion >= 3000
        @cuda(cudnn,cudnnConvolutionBackwardFilter_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,workSpaceSizeInBytes,Ref(T(beta)),FD(dw),dw)
    else
        @cuda(cudnn,cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,       Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),Ref(T(beta)),FD(dw),dw)
    end
    return dw
end

@primitive conv4(w,x; o...),dy  conv4w(w,x,dy;o...)  conv4x(w,x,dy;o...)
@zerograd  conv4x(w,x,dy;o...)
@zerograd  conv4w(w,x,dy;o...)

"""

    pool(x; kwargs...)

Compute pooling of input values (i.e., the maximum or average of
several adjacent values) to produce an output with smaller height
and/or width.

Currently 4 or 5 dimensional KnetArrays with `Float32` or `Float64`
entries are supported.  If `x` has dimensions `(X1,X2,...,I,N)`, the
result `y` will have dimensions `(Y1,Y2,...,I,N)` where

    Yi=1+floor((Xi+2*padding[i]-window[i])/stride[i])

Here `I` is the number of input channels, `N` is the number of
instances, and `Xi,Yi` are spatial dimensions.  `window`, `padding`
and `stride` are keyword arguments that can be specified as a single
number (in which case they apply to all dimensions), or an array/tuple
with entries for each spatial dimension.

# Keywords:

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

"""
function pool{T}(x::KnetArray{T}; handle=cudnnhandle(), alpha=1,
                 o...) # window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    y = similar(x, pdims(x; o...))
    beta = 0
    @cuda(cudnn, cudnnPoolingForward,
          (Cptr, Cptr,      Ptr{T},    Cptr,Ptr{T},Ptr{T},   Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(T(alpha)),TD(x),x,    Ref(T(beta)),TD(y),y)
    return y
end

function poolx{T}(x::KnetArray{T},y::KnetArray{T},dy::KnetArray{T}; handle=cudnnhandle(), alpha=1, mode=0,
                  o...) # window=2, padding=0, stride=window, maxpoolingNanOpt=0
    if alpha!=1 && mode==0; error("Gradient of pool(alpha!=1,mode=0) broken in CUDNN"); end
    dx = similar(x)
    beta = 0
    @cuda(cudnn,cudnnPoolingBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,PD(x;mode=mode,o...),Ref(T(alpha)),TD(y),y,TD(dy),dy,TD(x),x,Ref(T(beta)),TD(dx),dx)
    return dx
end

@primitive pool(x;o...),dy,y  poolx(x,y,dy;o...)
@zerograd  poolx(x,y,dy;o...)

"""

Unpooling; `reverse` of pooling.

    x == pool(unpool(x;o...); o...)

"""
function unpool(x; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,x))
    y = similar(x,updims(x; window=window, o...))
    poolx(y,x,x.*w; o..., window=window, mode=1, alpha=1/alpha)
end

function unpoolx(dy; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,dy))
    pool(dy; o..., window=window, mode=1, alpha=1/alpha) * w
end

# @primitive unpool(x;o...),dy,y -pool(-dy;o...)
@primitive  unpool(x;o...),dy,y  unpoolx(dy;o...)

"""

Deconvolution; `reverse` of convolution.

"""
function deconv4(w,x; o...)
    y = similar(x,dcdims(w,x;o...))
    return conv4x(w,y,x;o...)
end

function deconv4w(w,x,dy; o...)
    return conv4w(w,dy,x;o...)
end

function deconv4x(w,x,dy; o...)
    return conv4(w,dy;o...)
end

@primitive deconv4(w,x; o...),dy,y  deconv4w(w,x,dy; o...)  deconv4x(w,x,dy; o...)
@zerograd deconv4w(w,x,dy; o...)
@zerograd deconv4x(w,x,dy; o...)

# cudnn descriptors

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

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

type CD; ptr
    function CD(w::KnetArray,x::KnetArray; padding=0, stride=1, upscale=1, mode=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        if cudnnVersion >= 4000
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        elseif cudnnVersion > 3000 # does not work when cudnnVersion==3000
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor_v3,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        else
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode)
        end
        cd = new(d[1])
        finalizer(cd, x->@cuda(cudnn,cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr))
        return cd
    end
end

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

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

# fill and reverse Cint array with padding etc. for cudnn calls
function cdsize(w, nd)
    if isa(w,Number)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch("$w $nd"))
    end
end

# convert padding etc. size to an Int array of the right dimension
function psize(p, x)
    nd = ndims(x)-2
    if isa(p,Number)
        fill(Int(p),nd)
    elseif length(p)==nd
        collect(Int,p)
    else
        throw(DimensionMismatch("psize: $p $nd"))
    end
end

DT(::KnetArray{Float32})=UInt32(0)
DT(::KnetArray{Float64})=UInt32(1)
DT(::KnetArray{Float16})=UInt32(2)

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

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

function dcdims(w,x; padding=0, stride=1, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            si*(size(x,i)-1) + size(w,i) - 2*pi
        elseif i == N-1
            size(w,N)
        else
            size(x,N)
        end
    end
end

function updims(x; window=2, padding=0, stride=window, o...)
    window = psize(window,x)
    stride = psize(stride,x)
    padding = psize(padding,x)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            (size(x,i)-1)*stride[i]+window[i]-2*padding[i]
        else
            size(x,i)
        end
    end
end

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

### CPU convolution using im2col from Mocha.jl

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

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

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

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

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

# Functions from conv.cpp:

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

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

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

    ### CPU pooling from Mocha.jl

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

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


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GDAL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GDAL")
JLLWrappers.@generate_main_file("GDAL", UUID("a7073274-a066-55f0-b90d-d619367d196c"))
end  # module GDAL_jll