text:
import GLFW
using ModernGL

include("gltools.jl")
include("rgb.jl")

type MouseState
    x::Float64
    y::Float64
    pressed::Bool
end

function display(width, height, f; title="Julia")
    # Width and height are in screen coordinates, not pixels

    GLFW.Init()

    @osx_only begin
        # OS X-specific GLFW hints to specify the correct
        # version of OpenGL at context creation time
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
    end

    # Create a window with an OpenGL context
    window = GLFW.CreateWindow(width, height, title)

    # Make the window's context current
    GLFW.MakeContextCurrent(window)

    # Initialize a vertex array
    vao = glGenVertexArray()
    glBindVertexArray(vao)
    glCheckError("binding vertex array")

    # Initialize a buffer for vertices
    vbo = glGenBuffer()
    glBindBuffer(GL_ARRAY_BUFFER, vbo)
    glCheckError("binding buffer")

    # Bind the data to the VBO
    # Specifies a fullscreen rectangle in standard coordinates.
    # [top left, top right, bottom left, bottom right]
    data = GLfloat[1, -1, 1, 1, -1, -1, -1, 1]
    glBufferData(GL_ARRAY_BUFFER, sizeof(data), data, GL_STATIC_DRAW)
    glCheckError("setting buffer data")

    # Initialize a texture for color data
    tex = glGenTexture()
    glBindTexture(GL_TEXTURE_2D, tex)
    glTexStorage2D(GL_TEXTURE_2D, 1, GL_RGB8, width, height)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
    glCheckError("initializing texture")

    # Create and initialize shaders
    program = begin
        const vsh = """
            #version 330
            in vec2 position;

            void main() {
                gl_Position = vec4(position, 0.0, 1.0);
            }
        """

        const fsh = """
            #version 330
            out vec4 outColor;
            uniform vec2 screenDimsInPixels;
            uniform sampler2D tex;
            uniform float zoom;

            void main() {
                vec2 pos = gl_FragCoord.xy / screenDimsInPixels;
                vec2 zeroToCenter = vec2(0.5, 0.5);
                outColor = texture(tex, (pos - zeroToCenter) * zoom + zeroToCenter);
            }
        """

        vertexShader = createShader(vsh, GL_VERTEX_SHADER)
        fragmentShader = createShader(fsh, GL_FRAGMENT_SHADER)
        program = createShaderProgram(vertexShader, fragmentShader) do prog
            glBindFragDataLocation(prog, 0, "outColor")
        end

        glUseProgram(program)

        screenDimsInPixelsUniform = glGetUniformLocation(program, "screenDimsInPixels")
        @assert screenDimsInPixelsUniform > -1
        (pixelwidth::Float64, pixelheight::Float64) = GLFW.GetFramebufferSize(window)
        glUniform2f(screenDimsInPixelsUniform, pixelwidth, pixelheight)

        textureUniform = glGetUniformLocation(program, "tex")
        @assert textureUniform > -1

        zoomUniform = glGetUniformLocation(program, "zoom")
        @assert zoomUniform > -1
        glUniform1f(zoomUniform, 1.0)

        positionAttribute = glGetAttribLocation(program, "position")
        @assert positionAttribute > -1
        glEnableVertexAttribArray(positionAttribute)
        glVertexAttribPointer(positionAttribute, 2, GL_FLOAT, false, 0, 0)

        program
    end

    function drawtexture()
        # Activate and bind the texture
        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, tex)

        # The uniform value for a sampler refers to the texture unit (0 refers to GL_TEXTURE0)
        glUniform1i(textureUniform, 0)

        # Draw the textured rectangle
        glBindVertexArray(vao)
        glBindBuffer(GL_ARRAY_BUFFER, vbo)
        glDrawArrays(GL_TRIANGLE_STRIP, 0, 4)
    end

    regenerate(f, args...) = global generator = @task f(zeros(GLPixel, width, height), args...)
    regenerate(f)

    GLFW.SetKeyCallback(window, (key, scancode, action, mods) -> begin
        if action == GLFW.RELEASE
            key = uppercase(convert(Char, key))
            key == 'R' && regenerate(f)
            key == 'Q' && GLFW.SetWindowShouldClose(window, GL_TRUE)
        end
    end)

    # For now, disable garbage collection due to the inefficiency of the pixel-placing code
    gc_disable()

    zoom::GLfloat = 1.0

    # Event loop
    while !GLFW.WindowShouldClose(window)
        # Clear the screen
        glClearColor(0.0, 0.0, 0.0, 1.0)
        glClear(GL_COLOR_BUFFER_BIT)

        if !istaskdone(generator)
            (mx, my) = GLFW.GetCursorPos(window)
            mx = clamp(mx, 0, width) / width
            my = clamp(my, 0, height) / height
            pressed = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_LEFT) == GLFW.PRESS
            glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, consume(generator, MouseState(mx, my, pressed)))
        end

        glUniform1f(zoomUniform, zoom)

        # The pixel dimensions of the window may change here, e.g.
        # if the window is dragged from a high-dpi monitor to a low-dpi one.
        (pixelwidth, pixelheight) = GLFW.GetFramebufferSize(window)
        glUniform2f(screenDimsInPixelsUniform, pixelwidth, pixelheight)

        drawtexture()

        # Swap front and back buffers
        GLFW.SwapBuffers(window)

        # Poll for and process events
        GLFW.PollEvents()
    end

    gc_enable()
    gc()

    glDeleteTextures(1, [tex])
    GLFW.Terminate()
end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
export ProjMPO,
       LProj,
       RProj,
       product

mutable struct ProjMPO
  lpos::Int
  rpos::Int
  nsite::Int
  H::MPO
  LR::Vector{ITensor}
  ProjMPO(H::MPO) = new(0,length(H)+1,2,H,fill(ITensor(),length(H)))
end

nsite(pm::ProjMPO) = pm.nsite
length(pm::ProjMPO) = length(pm.H)

function LProj(pm::ProjMPO)::ITensor
  (pm.lpos <= 0) && return ITensor()
  return pm.LR[pm.lpos]
end

function RProj(pm::ProjMPO)::ITensor
  (pm.rpos >= length(pm)+1) && return ITensor()
  return pm.LR[pm.rpos]
end

function product(pm::ProjMPO,
                 v::ITensor)::ITensor
  Hv = v
  if isNull(LProj(pm))
    if !isNull(RProj(pm))
      Hv *= RProj(pm)
    end
    for j=pm.rpos-1:-1:pm.lpos+1
      Hv *= pm.H[j]
    end
  else #if LProj is not null
    Hv *= LProj(pm)
    for j=pm.lpos+1:pm.rpos-1
      Hv *= pm.H[j]
    end
    if !isNull(RProj(pm))
      Hv *= RProj(pm)
    end
  end
  return noprime(Hv)
end

(pm::ProjMPO)(v::ITensor) = product(pm,v)

function size(pm::ProjMPO)::Tuple{Int,Int}
  d = 1
  if !isNull(LProj(pm))
    for i in inds(LProj(pm))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  for j=pm.lpos+1:pm.rpos-1
    for i in inds(pm.H[j])
      plev(i) > 0 && (d *= dim(i))
    end
  end
  if !isNull(RProj(pm))
    for i in inds(RProj(pm))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  return (d,d)
end

function makeL!(pm::ProjMPO,
                psi::MPS,
                k::Int)
  while pm.lpos < k
    ll = pm.lpos
    if ll <= 0
      pm.LR[1] = psi[1]*pm.H[1]*dag(prime(psi[1]))
      pm.lpos = 1
    else
      pm.LR[ll+1] = pm.LR[ll]*psi[ll+1]*pm.H[ll+1]*dag(prime(psi[ll+1]))
      pm.lpos += 1
    end
  end
end

function makeR!(pm::ProjMPO,
                psi::MPS,
                k::Int)
  N = length(pm.H)
  while pm.rpos > k
    rl = pm.rpos
    if rl >= N+1
      pm.LR[N] = psi[N]*pm.H[N]*dag(prime(psi[N]))
      pm.rpos = N
    else
      pm.LR[rl-1] = pm.LR[rl]*psi[rl-1]*pm.H[rl-1]*dag(prime(psi[rl-1]))
      pm.rpos -= 1
    end
  end
end

function position!(pm::ProjMPO,
                   psi::MPS,
                   pos::Int)
  makeL!(pm,psi,pos-1)
  makeR!(pm,psi,pos+nsite(pm))

  #These next two lines are needed
  #when moving LProj and RProj backward
  pm.lpos = pos-1
  pm.rpos = pos+nsite(pm)
end

ext: jl


text:
using DifferentialEquations

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(lotka_volterra,u0,tspan,p)

sol = solve(prob)
using Plots
plot(sol)


text:
"""
    cubic(v, i, bc=:constant, constval=zero(eltype(v)))

Evaluate the vector `v` at the interpolated index `i` using cubic interpolation.
`bc` sets the boundary condition, which is used for extrapolation of the signal,
when `i` is outside the original support of `v`.

`bc` can be set to one of:
- :constant     - The function is `constval` outside the given indices
- :nearest      - The first and last values are extended
- :linear       - The first two and last two values are used to compute a slope
- :periodic     - The function is assumed to repeat infinitely
- :reflect_at   - The function is assumed to reflect at the first and last
                  points (so they don't repeat). This is equivalent to scipy's
                  `mirror` mode.
- :reflect_past - The function is assumed to reflect at a point 1/2 sample past
                  the first and last points, so they repeat. This is equivalent
                  to scipy's `reflect` mode.
"""
function cubic(v, i, bc=:linear, constval=zero(eltype(v)))
    # TODO: add :quadratic and :cubic boundary conditions that set a constant
    # 2nd and 3rd derivatives
    # TODO: for efficiency we also may want a `cubic_inbounds` that doesn't
    # check `i` and assumes we're within the signal. Definitely set up some
    # profiling and comparison to existing interpolation first.
    L = length(v)
    idxs = floor(Int, i) .+ (-1, 0, 1, 2)
    # generate y valuels, extrapolating if we're near the edge of the signal
    ys = if bc == :constant
        map(idx->1 <= idx <= L ? v[idx] : constval, idxs)
    elseif bc == :nearest
        map(x->v[clamp(x, 1, L)], idxs)
    elseif bc == :linear
        map(idxs) do idx
            if idx < 1
                v[1] + (v[2]-v[1]) * (idx-1)
            elseif idx > L
                v[L] + (v[L]-v[L-1]) * (idx-L)
            else
                v[idx]
            end
        end
    elseif bc == :periodic
        map(x->v[mod1(x, L)], idxs)
    elseif bc == :reflect_at # reflect at the edge
        # ABCDCB|ABCD|CBABCD
        map(idxs) do idx
            idx = mod1(idx, 2L-2)
            idx > L ? v[2L-idx] : v[idx]
        end
    elseif bc == :reflect_past # reflect 1/2 sample past the edge
        # DCBA|ABCD|DCBA
        map(idxs) do idx
            idx = mod1(idx, 2L)
            idx > L ? v[2L+1-idx] : v[idx]
        end
    else
        throw(ArgumentError("Unrecognized boundary condition $bc"))
    end

    # coeficients from http://paulbourke.net/miscellaneous/interpolation/
    # this is the Paul Breeuwsma / Catmull-Rom variant
    a0 = -0.5*ys[1] + 1.5*ys[2] - 1.5*ys[3] + 0.5*ys[4]
    a1 = ys[1] - 2.5*ys[2] + 2*ys[3] - 0.5*ys[4]
    a2 = -0.5*ys[1] + 0.5*ys[3]
    a3 = ys[2]

    i -= floor(i)
    a0*i^3 + a1*i^2 + a2*i + a3
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
environments = [
    :ant,
    :atlas,
    :cartpole,
    :halfcheetah,
    :hopper,
    :pendulum,
    :quadruped,
    :raiberthopper,
    :rexhopper,
    :walker,
    :block
]

for name in environments
    env = get_environment(name)
    @test true
end


text:
module DistributionsTest
using OnlineStats, Distributions, Base.Test
srand(02082016)

@testset "Distributions" begin
@testset "fitdistribution" begin
    y = rand(Beta(), 100)
    fitdistribution(Beta, y)

    y = rand(Categorical([.2, .2, .2, .4]), 100)
    fitdistribution(Categorical, y)

    y = rand(Cauchy(), 100)
    fitdistribution(Cauchy, y)

    y = rand(Gamma(2, 6), 100)
    fitdistribution(Gamma, y)

    y = randn(100)
    fitdistribution(Normal, y)

    x = rand(10)
    y = rand(Multinomial(5, x/sum(x)), 100)'
    fitdistribution(Multinomial, y)

    y = rand(MvNormal(zeros(4), diagm(ones(4))), 100)'
    fitdistribution(MvNormal, y)
end
@testset "Beta" begin
    y = rand(Beta(), 100)
    o = FitBeta(y)
    d = fit(Beta, y)
    @test mean(o)       ≈ mean(d)
    @test var(o)        ≈ var(d)
    @test params(o)[1]  ≈ params(d)[1]
    @test params(o)[2]  ≈ params(d)[2]
    @test nobs(o)       ≈ 100
end
@testset "Categorical" begin
    y = rand(Categorical([.2, .2, .2, .4]), 1000)
    o = FitCategorical(y)
    @test ncategories(o) == 4

    y = rand(Bool, 1000)
    o = FitCategorical(y)
    @test ncategories(o) == 2

    y = rand([:a, :b, :c, :d, :e, :f, :g], 1000)
    o = FitCategorical(y)
    @test ncategories(o) == 7
end
@testset "Cauchy" begin
    y = rand(Cauchy(), 10000)
    o = FitCauchy(y, LearningRate())
    fit!(o, y)
    @test_approx_eq_eps params(o)[1] 0.0 0.1
    @test_approx_eq_eps params(o)[2] 1.0 0.1
    @test nobs(o) == 2 * 10000
end
@testset "Gamma" begin
    y = rand(Gamma(2, 6), 100)
    o = FitGamma(y)
    @test mean(o) ≈ mean(y)
end
@testset "LogNormal" begin
    y = rand(LogNormal(), 100)
    o = FitLogNormal(y)
    @test_approx_eq_eps mean(o) mean(y) .1
end
@testset "Normal" begin
    y = randn(100)
    o = FitNormal(y)
    @test mean(o) ≈ mean(y)
    @test std(o) ≈ std(y)
    @test var(o) ≈ var(y)
end
@testset "Multinomial" begin
    x = rand(10)
    y = rand(Multinomial(5, x/sum(x)), 100)'
    o = FitMultinomial(y)
    @test mean(o) ≈ vec(mean(y, 1))
    @test nobs(o) == 100
end
@testset "MvNormal" begin
    y = rand(MvNormal(zeros(4), diagm(ones(4))), 100)'
    o = FitMvNormal(y)
    @test mean(o) ≈ vec(mean(y, 1))
    @test var(o) ≈ vec(var(y, 1))
    @test std(o) ≈ vec(std(y, 1))
    @test cov(o) ≈ cov(y)
    @test nobs(o) == 100
end
@testset "NormalMix" begin
    d = MixtureModel(Normal, [(0,1), (2,3), (4,5)])
    y = rand(d, 50_000)
    o = NormalMix(y, 3)

    fit!(o, y)
    fit!(o, y, 10)
    @test_approx_eq_eps mean(o) mean(y) .5
    @test_approx_eq_eps var(o)  var(y)  .5
    @test_approx_eq_eps std(o)  std(y)  .5
    @test length(componentwise_pdf(o, 0.5)) == 3
    @test ncomponents(o) == 3
    @test typeof(component(o, 1)) == Normal{Float64}
    @test length(probs(o)) == 3
    @test pdf(o, randn()) > 0
    @test 0 < cdf(o, randn()) < 1
    @test value(o) == o.value
    @test_approx_eq_eps quantile(o, [.25, .5, .75]) quantile(y, [.25, .5, .75]) .5
    quantile(o, collect(.01:.01:.99))

    fit!(o, y, 1)
    fit!(o, y, 2)
    fit!(o, y, 5)
    NormalMix(3, y)
end
end

end#module


text:
using FactCheck
include(Pkg.dir("PseudoGenomes", "src","pseudogenome-snp.jl"))
mm10path =joinpath(homedir(),"Dropbox/Data/GENOME_DATA/mm10full/mm10.fa")
vcf_path = joinpath( Pkg.dir(),"PseudoGenomes","testdata", "vcf-file1.vcf")
outpath = joinpath( tempdir(), "mm10.castsnps.fasta")
substitute_snps_in_fasta_file(mm10path,vcf_path,outpath)


text:
using DiffEqSensitivity,OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, DiffEqBase, ForwardDiff, Calculus, QuadGK,
      LinearAlgebra
using Test

fb = @ode_def_bare begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c
f = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c

p = [1.5,1.0,3.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
probb = ODEProblem(fb,[1.0;1.0],(0.0,10.0),p)
solb = solve(probb,Tsit5(),abstol=1e-14,reltol=1e-14)
sol_end = solve(probb,Tsit5(),abstol=1e-14,reltol=1e-14,
          save_everystep=false,save_start=false)

# Do a discrete adjoint problem
println("Calculate discrete adjoint sensitivities")
t = 0.0:0.5:10.0 # TODO: Add end point handling for callback
# g(t,u,i) = (1-u)^2/2, L2 away from 1
function dg(out,u,p,t,i)
  (out.=2.0.-u)
end

easy_res = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
easy_res2 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(quad=true,backsolve=false))
easy_res3 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(quad=false,backsolve=false))
easy_res4 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(backsolve=true))
easy_res5 = adjoint_sensitivities(sol,Kvaerno5(nlsolve=NLAnderson(), smooth_est=false),dg,t,abstol=1e-12,
                                 reltol=1e-10,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(backsolve=true))
easy_res6 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true,quad=true),
                                  checkpoints=sol.t[1:5:end])
easy_res7 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true,quad=false),
                                  checkpoints=sol.t[1:5:end])

adj_prob = ODEAdjointProblem(sol,dg,t)
adj_sol = solve(adj_prob,Tsit5(),abstol=1e-14,reltol=1e-14)
integrand = AdjointSensitivityIntegrand(sol,adj_sol)
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-12)

@test isapprox(res, easy_res, rtol = 1e-10)
@test isapprox(res, easy_res2, rtol = 1e-10)
@test isapprox(res, easy_res3, rtol = 1e-10)
@test isapprox(res, easy_res4, rtol = 1e-10)
@test isapprox(res, easy_res5, rtol = 1e-7)
@test isapprox(res, easy_res6, rtol = 1e-9)
@test isapprox(res, easy_res7, rtol = 1e-9)

println("Calculate adjoint sensitivities ")

easy_res8 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  save_everystep=false,save_start=false,
                                  sensealg=SensitivityAlg(backsolve=true))

@test isapprox(res, easy_res8, rtol = 1e-9)

end_only_res = adjoint_sensitivities(sol_end,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  save_everystep=false,save_start=false,
                                  sensealg=SensitivityAlg(backsolve=true))

@test isapprox(res, end_only_res, rtol = 1e-9)

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((2 .- A).^2)./2)
end
G([1.5,1.0,3.0])
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])

using Flux
res4 = Flux.Tracker.gradient(G,[1.5,1.0,3.0])[1]

using ReverseDiff
res5 = ReverseDiff.gradient(G,[1.5,1.0,3.0])

@test norm(res' .- res2) < 1e-8
@test norm(res' .- res3) < 1e-6
@test norm(res' .- res4) < 1e-6
@test norm(res' .- res5) < 1e-6

# Do a continuous adjoint problem

# Energy calculation
g(u,p,t) = (sum(u).^2) ./ 2
# Gradient of (u1 + u2)^2 / 2
function dg(out,u,p,t)
  out[1]= u[1] + u[2]
  out[2]= u[1] + u[2]
end

adj_prob = ODEAdjointProblem(sol,g,nothing,dg)
adj_sol = solve(adj_prob,Tsit5(),abstol=1e-14,reltol=1e-10)
integrand = AdjointSensitivityIntegrand(sol,adj_sol)
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-10)

println("Test the `adjoint_sensitivities` utility function")
easy_res = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
println("2")
easy_res2 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(quad=false))
println("3")
easy_res3 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
println("4")
easy_res4 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(autodiff=false))
println("5")
easy_res5 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true))
@test norm(easy_res .- res) < 1e-8
@test norm(easy_res2 .- res) < 1e-8
@test norm(easy_res3 .- res) < 1e-8
@test norm(easy_res4 .- res) < 1e-8
@test norm(easy_res5 .- res) < 1e-8

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=eltype(p).(prob.u0),p=p,
                    tspan=eltype(p).(prob.tspan))
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14)
  res,err = quadgk((t)-> (sum(sol(t)).^2)./2,0.0,10.0,atol=1e-14,rtol=1e-10)
  res
end
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])

@test norm(res' .- res2) < 1e-8
@test norm(res' .- res3) < 1e-6

# Buffer length test
f = (du, u, p, t) -> du .= 0
p = zeros(3); u = zeros(50)
prob = ODEProblem(f,u,(0.0,10.0),p)
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test_nowarn res = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

println("Adjoints of u0")

function dg(out,u,p,t,i)
  out .= 1 .- u
end

ū0,adj = adjoint_sensitivities_u0(sol,Tsit5(),dg,t,abstol=1e-14,
                         reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

adj2 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                        reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

ū02,adj3 = adjoint_sensitivities_u0(sol,Tsit5(),dg,t,abstol=1e-14,
                         sensealg=SensitivityAlg(backsolve=true),
                         reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

ū03,adj4 = adjoint_sensitivities_u0(sol,Tsit5(),dg,t,abstol=1e-14,
                        save_everystep=false, save_start=false,
                        sensealg=SensitivityAlg(backsolve=true),
                        reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

ū03,adj4 = adjoint_sensitivities_u0(sol,Tsit5(),dg,t,abstol=1e-14,
                        save_everystep=false, save_start=false,
                        sensealg=SensitivityAlg(backsolve=true),
                        reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

res = ForwardDiff.gradient(prob.u0) do u0
  tmp_prob = remake(prob,u0=u0)
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((1 .- A).^2)./2)
end

@test ū0 ≈ res rtol = 1e-10
@test ū02 ≈ res rtol = 1e-10
@test ū03 ≈ res rtol = 1e-10
@test adj ≈ adj2 rtol = 1e-10
@test adj ≈ adj3 rtol = 1e-10
@test adj ≈ adj4 rtol = 1e-10


text:
elapsed(i::Int, N::Int) = (i % N) == 0
elapsed(i::UnitRange, N::Int) = any([i...] .% N .== 0)

function TensorBoardLogger.log_value(logger::WBLogger, name::AbstractString, value::Real; step=0)
    info_dict = Dict(string(name) => value)
    wandb.log(info_dict, step=step)
end

@with_kw mutable struct LoggerParams
    dir::String = "log/"
    period::Int64 = 500
    use_wandb::Bool = false
    config::Union{Dict, Nothing} = nothing
    project::Union{AbstractString, Nothing} = nothing
    entity::Union{AbstractString, Nothing} = nothing
	notes::Union{AbstractString, Nothing} = nothing
    logger::Union{TBLogger, WBLogger} = use_wandb ? WBLogger(config=config, project=project, entity=entity, notes=notes) : TBLogger(dir, tb_increment)
    fns = Any[log_undiscounted_return(10), log_episode_averages([:r], period)]
    writeout::Dict{Int, Any} = Dict() # Other things to write to disk. Period => Function
    verbose::Bool = true
    sampler::Union{Sampler, Nothing, Vector} = nothing
end

Base.log(p::Nothing, i, data...; kwargs...)  = nothing

#Note that i can be an int or a unitrange
function Base.log(p::LoggerParams, i::Union{Int, UnitRange}, data...; 𝒮=nothing)

    # Write things to disc
    for (period, fn) in p.writeout
        elapsed(i, period) && fn(i=i[end], s=p.sampler, dir=p.dir, logger=p.logger)
    end

    # Save other run information
    !elapsed(i, p.period) && return
    i = i[end]
    p.verbose && print("Step: $i")
	dicts = Any[p.fns..., data...]

	if !isnothing(p.sampler)
    	π_explore = (p.sampler isa Vector) ?  first(p.sampler).agent.π_explore : p.sampler.agent.π_explore
		push!(dicts, log_exploration(π_explore))
	end

    for dict in dicts
        d = dict isa Function ? dict(s=p.sampler, i=i, 𝒮=𝒮) : dict
        for (k,v) in d
            p.verbose && print(", ", k, ": ", v)
            log_value(p.logger, string(k), v, step = i)
        end
    end
    p.verbose && println()
end
# @printf("%5d / %5d eps %0.3f |  avgR %1.3f | Loss %2.3e | Grad %2.3e | EvalR %1.3f \n",
                        #t, solver.max_steps, nt[1], avg100_reward, loss_val, grad_val, scores_eval)

function aggregate_info(infos)
    res = Dict()
    for k in unique(vcat([collect(keys(info)) for info in infos]...))
        res[k] = mean([info[k] for info in filter((x)->haskey(x, k), infos)])
    end
    res
end

# Built-in functions for logging common training things
log_performance(s::AbstractVector, name, fn; kwargs...) = Dict("$(name)/T$i" => fn(s[i]; kwargs...) for i=1:length(s))
log_performance(s::Sampler, name, fn; kwargs...) = Dict(name => fn(s; kwargs...))

log_discounted_return(Neps) = (;s, kwargs...) -> log_performance(s, "discounted_return", discounted_return, Neps=Neps)
log_undiscounted_return(Neps; name="undiscounted_return") = (;s, kwargs...) -> log_performance(s, name, undiscounted_return, Neps=Neps)
log_undiscounted_return(s, Neps; name="undiscounted_return") = (;kwargs...) -> log_performance(s, name, undiscounted_return, Neps=Neps)
log_failure(Neps) = (;s, kwargs...) -> log_performance(s, "failure_rate", failure, Neps=Neps)
log_metric_by_key(key, Neps) = (;s, kwargs...) -> log_performance(s, string(key), metric_by_key, key=key, Neps=Neps)

function log_metrics_by_key(keys, Neps; kwargs...)
    (;s,kwargs2...) -> begin
        vals = metrics_by_key(s, keys=keys, Neps=Neps, kwargs...)
        Dict(k=>v for (k,v) in zip(keys, vals))
    end
end

log_validation_error(loss, 𝒟_val; name="validation_error") = (;s, kwargs...) -> Dict(name => loss(s.π, 𝒟_val))

log_exploration(policy) = (;kwargs...) -> Dict()
log_exploration(policy::MixedPolicy; name = "eps") = (;i, kwargs...) -> Dict(name => policy.ϵ(i))
log_exploration(policy::GaussianNoiseExplorationPolicy; name = "noise_std") = (;i, kwargs...) -> Dict(name => policy.σ(i))
function log_exploration(policy::FirstExplorePolicy; name = "first_explore_on")
    (;i, kwargs...) -> begin
        d = Dict{String, Any}(name => i < policy.N)
        !isnothing(policy.after_policy) && merge!(d, log_exploration(policy.after_policy)(;i=1,kwargs...))
        d
    end
end

function log_episode_averages(keys, period)
    (;𝒮, kwargs...) -> begin
        d = Dict()
        if hasproperty(𝒮, :buffer)
            indices = get_last_N_indices(𝒮.buffer, period)
            for k in keys
                avg_val = sum(𝒮.buffer[k][indices]) / sum(𝒮.buffer[:episode_end][indices])
                d[Symbol(string("avg_", k))] = avg_val
            end
        end
        d
    end
end

function log_experience_sums(keys, period)
    (;𝒮, kwargs...) -> begin
        d = Dict()
        if hasproperty(𝒮, :buffer)
            indices = get_last_N_indices(𝒮.buffer, period)
            for k in keys
                avg_val = sum(𝒮.buffer[k][indices])
                d[Symbol(string("avg_", k))] = avg_val #TODO: Switch back to sum!
            end
        end
        d
    end
end

function save_gif(;base_name="demo", log_at_zero=false, log_wandb=true, append_rand=false, kwargs...)
    (;i, s, dir, logger) -> begin
        !log_at_zero && i==0 && return
        pref = append_rand ? replace(replace(string(now(), "_"), (":"=>"_")), "-"=>"_") : ""
        filename = string(dir,"/", pref, base_name, "_$i.gif")
        @info "writing gif to $filename"
        Crux.gif(s.mdp, s.agent.π, filename; kwargs...)

        if log_wandb && logger isa WBLogger
            vid = wandb.Video(filename, format="gif")
            wandb.log(Dict("gif"=>vid))
            rm(filename)
        end
    end
end


text:
# =====================================
# MNA 1
#
# system type: affine continuous system
# state dimension: 578
# =====================================
using ReachabilityBenchmarks, MathematicalSystems, MAT, SparseArrays

function mna1_model()
    file = matopen(@relpath "mna1.mat")

    # system matrix
    A = sparse(read(file, "A"))

    # affine term
    b = sparsevec(570:578, [fill(-0.1, 5); fill(-0.2, 4)], 578)

    # affine continuous system
    S = @system(x' = Ax + b)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using SafeTestsets

@safetestset "Simple" begin include("simple.jl") end


text:
using LinearAlgebra
using Test
using Plots

# Initial system state

# Time settings
Δt = 0.1
t_max = 100
t = Int(div(t_max, Δt))

# Masses of bodies
m = [0.999, 0.0005, 0.0005]
M = sum(m)
n = length(m)

# Gravitational constant
G = 1.0

# Matrices of shape [body, time, value]
R = zeros(n, t, 2)
V = zeros(n, t, 2)
A = zeros(n, t, 2)

# Initial values for position and velocity
R[:, 1, :] = [0 0; 1 0; 2 0]
V[:, 1, :] = [0.25 0; 0 1; 0 0.7]

"""
 Compute acceleration vector for timestep `t` and body `i`

# Arguments
- `R`: Matrix of shape `[body, time, position vector]`
- `G`: Gravitational constant
- `m`: Vector of masses of bodies
- `t`: Index of time step into `R`
- `i`: Index of body into `R`/`G`/`m`

# Returns
``a_i(t)`` Acceleration vector of body `i` and timestep `t`
"""
function acceleration_i(R, G, m, t, i)
    a_new = []

    for j = 1:length(R[:,1,1])
        if i == j
            continue
        end
        r_i = R[i, t, :]
        r_j = R[j, t, :]
        r_ij = r_j - r_i

        a_ij = r_ij * G * m[j] / norm(r_ij)^3
        push!(a_new, a_ij)
    end
    return sum(a_new)
end

"""
Compute complete acceleration matrix for timestep `t`

See `acceleration_i` for details on params

# Returns
``A(t)`` Acceleration matrix of timestep `t`
"""
function acceleration(R, G, m, t)
    return hcat(map(i->acceleration_i(R, G, m, t, i), 1:n)...)'
end

@testset "acceleration" begin
    time = 2
    m = [0.999, 0.0005, 0.0005]
    n = 3
    R = zeros(n, time, 2)
    R[:, 1, :] = [0 0; 1 0; 2 0]
    G = 1.0
    t = 1
    i = 1
    @test acceleration_i(R, G, m, t, i) == [0.000625, 0.0]
    i = 2
    @test acceleration_i(R, G, m, t, i) == [-0.9985, 0.0]
    i = 3
    @test acceleration_i(R, G, m, t, i) == [-0.25025, 0.0]

    @test acceleration(R, G, m, t) == [
        0.000625 0.0;
        -0.9985 0.0;
        -0.25025 0.0]
end

"""
Euler-Method: Compute the new position and
velocity for all the bodys in the system for the
next timestep.

# Arguments
- `Δt`: timestep length
- `t`: Index of time step into `V`, `R` and `A`
- `V`: Matrix of shape `[body, time, velocity vector]`
- `R`: Matrix of shape `[body, time, position vector]`
- `A`: Matrix of shape `[body, time, acceleration vector]`

# Returns
``v_new(t), r_new(t)`` Velocitiy and position of one body
at the next timestep
"""
function euler_method(Δt, t, V, R, A)
    v_new = V[:,t,:] + A[:,t,:] * Δt
    r_new = R[:,t,:] + V[:,t,:] * Δt

    return v_new, r_new
end

"""
Euler-Cormer- Method: Compute the new position and
velocity for all the bodys in the system for the
next timestep.

# Arguments
- `Δt`: timestep length
- `t`: Index of time step into `V`, `R` and `A`
- `V`: Matrix of shape `[body, time, velocity vector]`
- `R`: Matrix of shape `[body, time, position vector]`
- `A`: Matrix of shape `[body, time, acceleration vector]`

# Returns
``v_new(t), r_new(t)`` Velocitiy and position of one body
at the next timestep
"""
function euler_cormer(Δt, t, V, R, A)
    v_new = V[:,t,:] + A[:,t,:] * Δt
    r_new = R[:,t,:] + 0.5 * (V[:,t,:] + v_new) * Δt

    return v_new, r_new
end

"""
Verlet-Algotihm: Compute the new position and
velocity for all the bodys in the system for the
next timestep. Speical case for t = t_start.
(Here t_start = 1).

# Arguments
- `Δt`: timestep length
- `t`: Index of time step into `V`, `R` and `A`
- `V`: Matrix of shape `[body, time, velocity vector]`
- `R`: Matrix of shape `[body, time, position vector]`
- `A`: Matrix of shape `[body, time, acceleration vector]`

# Returns
``v_new(t), r_new(t)`` Velocitiy and position of one body
at the next timestep
"""
function verlet(Δt, t, V, R, A)

    if t == 1
        r_help = R[:,1,:] - V[:,1,:] * Δt + 0.5 * A[:,1,:] * Δt^2
        r_new = 2 * R[:,t,:] - r_help + A[:,t,:] * Δt^2
        v_new = (r_new - r_help) / (2 * Δt)
    else
        r_new = 2 * R[:,t,:] - R[:,t-1,:] + A[:,t,:] * Δt^2
        v_new = (r_new - R[:,t-1,:]) / (2 * Δt)
    end

    return v_new, r_new
end

#=
iterate over solvers; solve our example problem with each one
and animate each solution
=#

function animate_solver(solver, R, G, m, t, V, Δt)
    for time = 1:(t - 1)
        A[:, time, :] = acceleration(R, G, m, time)
        v_new, r_new = solver(Δt, time, V, R, A)

        V[:, time + 1, :] = v_new
        R[:, time + 1, :] = r_new
    end

    anim = @animate for t = 1:size(R, 2)
        Rx = [R[i, t, 1] for i in 1:size(R, 1)]
        Ry = [R[i, t, 2] for i in 1:size(R, 1)]
        scatter(Rx', Ry', xlims = (-5, 20), ylims = (-5, 5), label = ["Body $i" for i in 1:size(R, 1)])
    end every 2

    name = string(solver)
    gif(anim, "R_$name.gif", fps = 30)
end

for solver in [verlet, euler_method, euler_cormer]
    animate_solver(solver, R, G, m, t, V, Δt)
end

max_stars_repo_name: JoshuaSimon/N-Body-Problem

max_issues_repo_name: JoshuaSimon/N-Body-Problem

max_forks_repo_name: JoshuaSimon/N-Body-Problem


text:
using CSyntax
using LinearAlgebra

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))
include(joinpath(@__DIR__, "camera.jl"))

# init window
width, height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# camera
camera = PerspectiveCamera()
setposition!(camera, [0.0, 0.0, 2.0])

# vertex data
points = GLfloat[-0.5, -0.5, 0.0,
                  0.5, -0.5, 0.0,
                  0.5,  0.5, 0.0,
                  0.5,  0.5, 0.0,
                 -0.5,  0.5, 0.0,
                 -0.5, -0.5, 0.0]

texcoords = GLfloat[0.0, 0.0,
                    1.0, 0.0,
                    1.0, 1.0,
                    1.0, 1.0,
                    0.0, 1.0,
                    0.0, 0.0]

# create VBO
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

texcoords_vbo = GLuint(0)
@c glGenBuffers(1, &texcoords_vbo)
glBindBuffer(GL_ARRAY_BUFFER, texcoords_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(texcoords), texcoords, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoords_vbo)
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)

# create shader program
vert_shader = createshader(joinpath(@__DIR__, "blending.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "blending.frag"), GL_FRAGMENT_SHADER)
shader_prog = createprogram(vert_shader, frag_shader)

model_loc = glGetUniformLocation(shader_prog, "model")
view_loc = glGetUniformLocation(shader_prog, "view")
proj_loc = glGetUniformLocation(shader_prog, "proj")
glUseProgram(shader_prog)
glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))

# load texture
glActiveTexture(GL_TEXTURE0)
tex_a = load_texture(joinpath(@__DIR__, "blob.png"))
tex_b = load_texture(joinpath(@__DIR__, "blob2.png"))

glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
glEnable(GL_BLEND)
# glDisable(GL_DEPTH_TEST)

# enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CCW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

let
updatefps = FPSCounter()
# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)

    glDepthMask(GL_FALSE)

    glBindTexture(GL_TEXTURE_2D, tex_a)
    model_matrix = Matrix{GLfloat}(I, 4, 4)
    glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_matrix)
    glDrawArrays(GL_TRIANGLES, 0, 6)

    glBindTexture(GL_TEXTURE_2D, tex_b)
    model_matrix[1,4] = 0.5
    model_matrix[3,4] = 0.1
    glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_matrix)
    glDrawArrays(GL_TRIANGLES, 0, 6)

    glDepthMask(GL_TRUE)

    # check and call events
    GLFW.PollEvents()
    # move camera
    updatecamera!(window, camera)
    glUseProgram(shader_prog)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
"""
    stairs(x, y; kw...)

# Description

Draws a staircase plot on a new canvas.

The first (optional) vector `x` should contain the horizontal
positions for all the points. The second vector `y` should then
contain the corresponding vertical positions respectively.
This means that the two vectors must be of the same length and ordering.

# Usage

    stairs(x, y; $(keywords((style = :post,); add = (:canvas,))))

# Arguments

$(arguments(
    (; style = "specifies where the transition of the stair takes place (can be either `:pre` or `:post`)");
    add = (:x, :y, :canvas)
))

# Returns

A plot object of type `Plot{T<:Canvas}`.

# Author(s)

- Christof Stocker (github.com/Evizero)
- Dominique (github.com/dpo)

# Examples

```julia-repl
julia> stairs([1, 2, 4, 7, 8], [1, 3, 4, 2, 7], style = :post, title = "My Staircase Plot")
                 My Staircase Plot
     ┌────────────────────────────────────────┐
   7 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡄⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⢸│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠧⠤⠤⠤⠤⠼│
     │⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   1 │⣀⣀⣀⣀⣀⣸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
     └────────────────────────────────────────┘
     1                                        8
```

# See also

[`Plot`](@ref), [`scatterplot`](@ref), [`lineplot`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function stairs(X::AbstractVector, Y::AbstractVector; style::Symbol = :post, kw...)
    x_vex, y_vex = compute_stair_lines(X, Y, style)
    lineplot(x_vex, y_vex; kw...)
end

function stairs!(
    plot::Plot{<:Canvas},
    X::AbstractVector,
    Y::AbstractVector;
    style::Symbol = :post,
    kw...,
)
    x_vex, y_vex = compute_stair_lines(X, Y, style)
    lineplot!(plot, x_vex, y_vex; kw...)
end

function compute_stair_lines(X::AbstractVector, Y::AbstractVector, style::Symbol)
    if style == :post
        x_vex = similar(X, 2length(X) - 1)
        y_vex = similar(Y, 2length(X) - 1)
        x_vex[1] = X[1]
        y_vex[1] = Y[1]
        o = 0
        for i in 2:(length(X))
            x_vex[i + o] = X[i]
            x_vex[i + o + 1] = X[i]
            y_vex[i + o] = Y[i - 1]
            y_vex[i + o + 1] = Y[i]
            o += 1
        end
        return x_vex, y_vex
    elseif style == :pre
        x_vex = zeros(2length(X) - 1)
        y_vex = zeros(2length(X) - 1)
        x_vex[1] = X[1]
        y_vex[1] = Y[1]
        o = 0
        for i in 2:(length(X))
            x_vex[i + o] = X[i - 1]
            x_vex[i + o + 1] = X[i]
            y_vex[i + o] = Y[i]
            y_vex[i + o + 1] = Y[i]
            o += 1
        end
        return x_vex, y_vex
    end
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

half_pi(::Type{T}) where T = pi_interval(T) / 2
half_pi(x::T) where T<:AbstractFloat = half_pi(T)

two_pi(::Type{T})  where T = pi_interval(T) * 2

range_atan2(::Type{T}) where {T<:Real} = Interval(-(pi_interval(T).hi), pi_interval(T).hi)
half_range_atan2(::Type{T}) where {T} = (temp = half_pi(T); Interval(-(temp.hi), temp.hi) )
pos_range_atan2(::Type{T}) where {T<:Real} = Interval(zero(T), pi_interval(T).hi)

doc"""Finds the quadrant(s) corresponding to a given floating-point
number. The quadrants are labelled as 0 for x ∈ [0, π/2], etc.
For numbers very near a boundary of the quadrant, a tuple of two quadrants
is returned. The minimum or maximum must then be chosen appropriately.

This is a rather indirect way to determine if π/2 and 3π/2 are contained
in the interval; cf. the formula for sine of an interval in
Tucker, *Validated Numerics*."""

function find_quadrants(x::T) where {T<:AbstractFloat}
    temp = atomic(Interval{T}, x) / half_pi(x)
    (floor(temp.lo), floor(temp.hi))
end

function sin(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval{T}(-1, 1)

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return whole_range
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range  # in same quadrant but separated by almost 2pi
        lo = @round(sin(a.lo), sin(a.lo)) # Interval(sin(a.lo, RoundDown), sin(a.lo, RoundUp))
        hi = @round(sin(a.hi), sin(a.hi)) # Interval(sin(a.hi, RoundDown), sin(a.hi, RoundUp))
        return hull(lo, hi)

    elseif lo_quadrant==3 && hi_quadrant==0
        return @round(sin(a.lo), sin(a.hi)) # Interval(sin(a.lo, RoundDown), sin(a.hi, RoundUp))

    elseif lo_quadrant==1 && hi_quadrant==2
        return @round(sin(a.hi), sin(a.lo)) # Interval(sin(a.hi, RoundDown), sin(a.lo, RoundUp))

    elseif ( lo_quadrant == 0 || lo_quadrant==3 ) && ( hi_quadrant==1 || hi_quadrant==2 )
        return @round(min(sin(a.lo), sin(a.hi)), 1)
        # Interval(min(sin(a.lo, RoundDown), sin(a.hi, RoundDown)), one(T))

    elseif ( lo_quadrant == 1 || lo_quadrant==2 ) && ( hi_quadrant==3 || hi_quadrant==0 )
        return @round(-1, max(sin(a.lo), sin(a.hi)))
        # Interval(-one(T), max(sin(a.lo, RoundUp), sin(a.hi, RoundUp)))

    else #if( lo_quadrant == 0 && hi_quadrant==3 ) || ( lo_quadrant == 2 && hi_quadrant==1 )
        return whole_range
    end
end

function cos(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval(-one(T), one(T))

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return Interval(-one(T), one(T))
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant # Interval limits in the same quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range
        lo = @round(cos(a.lo), cos(a.lo))
        hi = @round(cos(a.hi), cos(a.hi))
        return hull(lo, hi)

    elseif lo_quadrant == 2 && hi_quadrant==3
        return @round(cos(a.lo), cos(a.hi))

    elseif lo_quadrant == 0 && hi_quadrant==1
        return @round(cos(a.hi), cos(a.lo))

    elseif ( lo_quadrant == 2 || lo_quadrant==3 ) && ( hi_quadrant==0 || hi_quadrant==1 )
        return @round(min(cos(a.lo), cos(a.hi)), 1)

    elseif ( lo_quadrant == 0 || lo_quadrant==1 ) && ( hi_quadrant==2 || hi_quadrant==3 )
        return @round(-1, max(cos(a.lo), cos(a.hi)))

    else#if ( lo_quadrant == 3 && hi_quadrant==2 ) || ( lo_quadrant == 1 && hi_quadrant==0 )
        return whole_range
    end
end

function tan(a::Interval{T}) where T
    isempty(a) && return a

    diam(a) > pi_interval(T).lo && return entireinterval(a)

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    lo_quadrant_mod = mod(lo_quadrant, 2)
    hi_quadrant_mod = mod(hi_quadrant, 2)

    if lo_quadrant_mod == 0 && hi_quadrant_mod == 1
        # check if really contains singularity:
        if hi_quadrant * half_pi(T) ⊆ a
            return entireinterval(a)  # crosses singularity
        end

    elseif lo_quadrant_mod == hi_quadrant_mod && hi_quadrant > lo_quadrant
        # must cross singularity
        return entireinterval(a)

    end

    # @show a.lo, a.hi
    # @show tan(a.lo), tan(a.hi)

    return @round(tan(a.lo), tan(a.hi))
end

function asin(a::Interval{T}) where T

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    return @round(asin(a.lo), asin(a.hi))
end

function acos(a::Interval{T}) where T

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    return @round(acos(a.hi), acos(a.lo))
end

function atan(a::Interval{T}) where T
    isempty(a) && return a

    return @round(atan(a.lo), atan(a.hi))
end

#atan2{T<:Real, S<:Real}(y::Interval{T}, x::Interval{S}) = atan2(promote(y, x)...)

function atan2(y::Interval{Float64}, x::Interval{Float64})
    (isempty(y) || isempty(x)) && return emptyinterval(Float64)

    atomic(Interval{Float64}, atan2(big53(y), big53(x)))
end

function atan2(y::Interval{BigFloat}, x::Interval{BigFloat})
    (isempty(y) || isempty(x)) && return emptyinterval(BigFloat)

    T = BigFloat

    # Prevent nonsense results when y has a signed zero:
    if y.lo == zero(T)
        y = Interval(zero(T), y.hi)
    end
    if y.hi == zero(T)
        y = Interval(y.lo, zero(T))
    end

    # Check cases based on x
    if x == zero(x)

        y == zero(y) && return emptyinterval(T)
        y.lo ≥ zero(T) && return half_pi(T)
        y.hi ≤ zero(T) && return -half_pi(T)
        return half_range_atan2(T)

    elseif x.lo > zero(T)

        y == zero(y) && return y
        y.lo ≥ zero(T) &&
            return @round(atan2(y.lo, x.hi), atan2(y.hi, x.lo)) # refinement lo bound
        y.hi ≤ zero(T) &&
            return @round(atan2(y.lo, x.lo), atan2(y.hi, x.hi))
        return @round(atan2(y.lo, x.lo), atan2(y.hi, x.lo))

    elseif x.hi < zero(T)

        y == zero(y) && return pi_interval(T)
        y.lo ≥ zero(T) &&
            return @round(atan2(y.hi, x.hi), atan2(y.lo, x.lo))
        y.hi < zero(T) &&
            return @round(atan2(y.hi, x.lo), atan2(y.lo, x.hi))
        return range_atan2(T)

    else # zero(T) ∈ x

        if x.lo == zero(T)
            y == zero(y) && return y

            y.lo ≥ zero(T) && return @round(atan2(y.lo, x.hi), half_range_atan2(BigFloat).hi)

            y.hi ≤ zero(T) && return @round(half_range_atan2(BigFloat).lo, atan2(y.hi, x.hi))
            return half_range_atan2(T)

        elseif x.hi == zero(T)
            y == zero(y) && return pi_interval(T)
            y.lo ≥ zero(T) && return @round(half_pi(BigFloat).lo, atan2(y.lo, x.lo))
            y.hi < zero(T) && return @round(atan2(y.hi, x.lo), -(half_pi(BigFloat).lo))
            return range_atan2(T)
        else
            y.lo ≥ zero(T) &&
                return @round(atan2(y.lo, x.hi), atan2(y.lo, x.lo))
            y.hi < zero(T) &&
                return @round(atan2(y.hi, x.lo), atan2(y.hi, x.hi))
            return range_atan2(T)
        end

    end
end


text:
# "Public" types that represent the file formats. These are used
# to communicate results externally, but are generally avoided for
# internal operations because they trigger excessive specialization
# and inference failures.

## DataFormat:
"""
`DataFormat{sym}()` indicates a known binary or text format of kind `sym`,
where `sym` is always a symbol. For example, a .csv file might have
`DataFormat{:CSV}()`.

An easy way to write `DataFormat{:CSV}` is `format"CSV"`.
"""
struct DataFormat{sym} end

macro format_str(s)
    :(DataFormat{$(Expr(:quote, Symbol(s)))})
end

formatname(::Type{DataFormat{sym}}) where sym = sym

abstract type Formatted{F<:DataFormat} end  # A specific file or stream

formatname(::Formatted{F}) where F<:DataFormat = formatname(F)

## File:

"""
`File{fmt}(filename)` indicates that `filename` is a file of known
[`DataFormat`](@ref) `fmt`.  For example, `File{format"PNG"}(filename)` would indicate a PNG
file.

!!! compat
    `File{fmt}(filename)` requires FileIO 1.6 or higher. The deprecated syntax `File(fmt, filename)` works
    on all FileIO 1.x releases.
"""
struct File{F<:DataFormat, Name} <: Formatted{F}
    filename::Name
end
File{F}(file::File{F}) where F<:DataFormat = file
File{DataFormat{sym}}(@nospecialize(file::Formatted)) where sym = throw(ArgumentError("cannot change the format of $file to $sym"))
File{F}(file::AbstractString) where F<:DataFormat = File{F,String}(String(file)) # canonicalize to limit type-diversity
File{F}(file) where F<:DataFormat = File{F,typeof(file)}(file)

# The docs are separated from the definition because of https://github.com/JuliaLang/julia/issues/34122
filename(@nospecialize(f::File)) = f.filename
"""
`filename(file)` returns the filename associated with [`File`](@ref) `file`.
"""
filename(::File)

file_extension(@nospecialize(f::File)) = splitext(filename(f))[2]
"""
`file_extension(file)` returns the file extension associated with [`File`](@ref) `file`.
"""
file_extension(::File)

## Stream:

"""
`Stream{fmt}(io, filename=nothing)` indicates that the stream `io` is
written in known format [`DataFormat`](@ref) `fmt`.
For example, `Stream{format"PNG"}(io)` would indicate PNG format.
If known, the optional `filename` argument can
be used to improve error messages, etc.

!!! compat
    `Stream{fmt}(io, ...)` requires FileIO 1.6 or higher. The deprecated syntax `Stream(fmt, io, ...)` works
    on all FileIO 1.x releases.
"""
struct Stream{F <: DataFormat, IOtype <: IO, Name} <: Formatted{F}
    io::IOtype
    filename::Name
end

Stream{F,IOtype}(io::IO, filename::AbstractString) where {F<:DataFormat,IOtype} = Stream{F, IOtype, String}(io, String(filename))
Stream{F,IOtype}(io::IO, filename)                 where {F<:DataFormat,IOtype} = Stream{F, IOtype, typeof(filename)}(io, filename)
Stream{F,IOtype}(io::IO)                           where {F<:DataFormat,IOtype} = Stream{F, IOtype}(io, nothing)

Stream{F,IOtype}(file::Formatted{F}, io::IO) where {F<:DataFormat,IOtype} = Stream{F,IOtype}(io, filename(file))
Stream{F,IOtype}(@nospecialize(file::Formatted), io::IO) where {F<:DataFormat,IOtype} =
    throw(ArgumentError("cannot change the format of $file to $(formatname(F)::Symbol)"))

Stream{F}(io::IO, args...) where {F<:DataFormat} = Stream{F, typeof(io)}(io, args...)
Stream{F}(file::File, io::IO) where {F<:DataFormat} = Stream{F, typeof(io)}(file, io)
Stream(file::File{F}, io::IO) where {F<:DataFormat} = Stream{F}(io, filename(file))

stream(@nospecialize(s::Stream)) = s.io
"`stream(s)` returns the stream associated with [`Stream`](@ref) `s`"
stream(::Stream)

filename(@nospecialize(s::Stream)) = s.filename
"""
`filename(stream)` returns a string of the filename
associated with [`Stream`](@ref) `stream`, or nothing if there is no file associated.
"""
filename(::Stream)

function file_extension(@nospecialize(f::Stream))
    fname = filename(f)
    (fname === nothing) && return nothing
    splitext(fname)[2]
end
"""
`file_extension(file)` returns a nullable-string for the file extension associated with [`Stream`](@ref) `stream`.
"""
file_extension(::Stream)

# Note this closes the stream. It's useful when you've opened
# the file to check the magic bytes, but don't want to leave
# a dangling stream.
function file!(strm::Stream{F}) where F
    f = filename(strm)
    f === nothing && error("filename unknown")
    close(strm.io)
    File{F}(f)
end

# Implement standard I/O operations for File and Stream
@inline function Base.open(@nospecialize(file::File{F}), @nospecialize(args...)) where F<:DataFormat
    fn = filename(file)
    Stream{F}(open(fn, args...), abspath(fn))
end
Base.close(@nospecialize(s::Stream)) = close(stream(s))

Base.position(@nospecialize(s::Stream)) = position(stream(s))
Base.seek(@nospecialize(s::Stream), offset::Integer) = (seek(stream(s), offset); s)
Base.seekstart(@nospecialize(s::Stream)) = (seekstart(stream(s)); s)
Base.seekend(@nospecialize(s::Stream)) = (seekend(stream(s)); s)
Base.skip(@nospecialize(s::Stream), offset::Integer) = (skip(stream(s), offset); s)
Base.eof(s::Stream) = eof(stream(s))

@inline Base.read(s::Stream, args...)  = read(stream(s), args...)
Base.read!(s::Stream, array::Array) = read!(stream(s), array)
@inline Base.write(s::Stream, args...) = write(stream(s), args...)
# Note: we can't sensibly support the all keyword. If you need that,
# call read(stream(s), ...; all=value) manually
Base.readbytes!(s::Stream, b) = readbytes!(stream(s), b)
Base.readbytes!(s::Stream, b, nb) = readbytes!(stream(s), b, nb)
Base.read(s::Stream) = read(stream(s))
Base.read(s::Stream, nb) = read(stream(s), nb)
Base.flush(s::Stream) = flush(stream(s))

Base.isreadonly(s::Stream) = isreadonly(stream(s))
Base.isopen(s::Stream) = isopen(stream(s))


text:
# functions to convert Julia expression to MATLAB statement strings

mstatement(s::Symbol) = string(s)

function write_mstatement(io::IO, s)
    print(io, string(s))
end

function write_mstatement(io::IO, s::ASCIIString)
    print(io, "'")
    print(io, s)
    print(io, "'")
end

function write_mstatement(io::IO, s::Char)
    print(io, "'")
    print(io, string(s))
    print(io, "'")
end

# Hack to support Julia 0.2 and 0.3
const matlab_infix_ops = union!(Set{Symbol}(),
    [:+, :-, :*, :/, :\, :^,
     :.+, :.-, :.*, :./, :.\, :.^,
     :|, :&])

function write_mstatement(io::IO, ex::Expr)
    h = ex.head
    a = ex.args
    na = length(a)

    if h == :(=)
        write_mstatement(io, a[1])
        print(io, " = ")
        write_mstatement(io, a[2])

    elseif h == :call || h == :ref
        f = a[1]

        if na == 1
            write_mstatement(io, f)
            print(io, "()")

        elseif na == 2
            write_mstatement(io, f)
            print(io, "(")
            write_mstatement(io, a[2])
            print(io, ")")

        elseif na == 3
            if (f in matlab_infix_ops)
                print(io, "(")
                write_mstatement(io, a[2])
                print(io, ") ")
                write_mstatement(io, f)
                print(io, " (")
                write_mstatement(io, a[3])
                print(io, ")")
            else
                write_mstatement(io, f)
                print(io, "(")
                write_mstatement(io, a[2])
                print(io, ", ")
                write_mstatement(io, a[3])
                print(io, ")")
            end

        else
            write_mstatement(io, f)
            print(io, "(")
            for i = 2 : na - 1
                write_mstatement(io, a[i])
                print(io, ", ")
            end
            write_mstatement(io, a[na])
            print(io, ")")
        end

    elseif h == :vcat || h == :hcat || h == :row || h == :tuple || h == :cell1d
        print(io, h == :cell1d ? "{" : "[")
        if na > 0
            sep = h == :vcat ? "; " : ", "
            for i = 1 : na-1
                write_mstatement(io, a[i])
                print(io, sep)
            end
            write_mstatement(io, a[na])
        end
        print(io, h == :cell1d ? "}" : "]")

    elseif h == :cell2d
        print(io, "{")
        if na > 0
            rows = a[1]
            cols = a[2]
            for i = 1:rows
                if i != 1
                    print(io, "; ")
                end
                for j = 0:cols-1
                    if j != 0
                        print(io, " ")
                    end
                    write_mstatement(io, a[2+j*rows+i])
                end
            end
        end
        print(io, "}")

    elseif h == :curly
        write_mstatement(io, a[1])
        print(io, "{")
        for i = 2 : na - 1
            write_mstatement(io, a[i])
            print(io, ", ")
        end
        write_mstatement(io, a[na])
        print(io, "}")

    elseif h == :.
        write_mstatement(io, a[1])
        print(io, ".")
        print(io, string(eval(a[2])))

    elseif h == :(:)
        write_mstatement(io, a[1])
        for i = 2:na
            print(io, ":")
            write_mstatement(io, a[i])
        end

    elseif h == symbol("'") || h == symbol(".'")
        print(io, "(")
        write_mstatement(io, a[1])
        print(io, ")")
        print(io, string(h))

    elseif h == :macrocall
        print(io, a[1])

    elseif h == :block
        for s in a
            write_mstatement(io, s)
            print(io, ";\n")
        end

    elseif h == :line
        nothing # ignore line numbers

    else
        error("Expression $(ex) is not recognized by mstatement.")
    end
end

function mstatement(ex::Expr)
    ss = IOBuffer()
    write_mstatement(ss, ex)
    bytestring(ss)
end


text:
@doc """ xgboost for train,val,test dataset without gridtune
""" ->
function model_xgboost(train, val, test, num_iter::Int64)
    train_X,train_Y = train[:,1:end-1],train[:,end]
    val_X,val_Y = val[:,1:end-1],val[:,end]
    test_X,test_Y = test[:,1:end-1],test[:,end]

    num_class = maximum(train_Y)+1
    m_tr,n_tr   = size(train_X)
    m_val,n_val = size(val_X)
    info("in the model, there are $num_class classes")
    info("train_X size is :$m_tr,$n_tr")
    info("val_X size is : $m_val, $n_tr")

    num_round = num_iter

    dtrain = DMatrix(train_X, label = train_Y)
    dval   = DMatrix(val_X,   label = val_Y)
    dtest  = DMatrix(test_X,  label = test_Y)

    watch_list = [(dtrain,"train"), (dval,"val")]
    param      = Dict{ASCIIString,Any}("max_depth"=>7,"eta"=>0.01,"nthread"=>50,
                                       "objective"=>"multi:softprob","silent"=>1,
                                       #"alpha"=>0.7,
                                       "sub_sample"=>0.9)#  ,"num_class"=>num_class)

    bst = xgboost(dtrain,num_round,watchlist=watch_list,param=param,num_class=num_class,
                  metrics=["mlogloss", "merror"],seed=2015)

    test_preds = XGBoost.predict(bst, test_X)

    yprob = reshape(test_preds, Int64(num_class), size(test_X,1))
    ind_prob = map(ind->indmax(yprob[:,ind]), 1:size(test_X,1))
    test_preds = ind_prob - 1

    return bst
end


text:
module DifferencesMatrix

using Test

import LinearAlgebra.diagm

function backwardDiff(dim::Int)
    return diagm(-1 => -ones(dim - 1), 0 => ones(dim))
end

@test backwardDiff(3) == [
  1 0 0
  -1 1 0
  0 -1 1
]

function forwardDiff(dim::Int)
    return diagm(0 => -ones(dim), 1 => ones(dim - 1))
end

@test forwardDiff(3) == [
  -1 1 0
  0 -1 1
  0 0 -1
]

@test backwardDiff(3) * forwardDiff(3) == [
  -1 1 0
  1 -2 1
  0 1 -2
]

function centeredDiff(dim::Int)
    return diagm(-1 => -ones(dim - 1), 0 => zeros(dim), 1 => ones(dim - 1))
end

@test centeredDiff(3) == backwardDiff(3) - (-forwardDiff(3))

export backwardDiff, forwardDiff, centeredDiff

end


text:
using HexapoleDataFit
using CoordinateTransformations

img = readImage("averaged_V+H_NO.dat")
rotimg = rotateImage(img, -102)
trnsimg = translateImage(rotimg, (50, 0))
trnssimg = translateImage(trnsimg, (0, 50))

showImage(translateImage(rotimg, CartesianIndex(0, 50)))
showImage(translateImage(rotimg, CartesianIndex(50, 0)))

# show (in Juno) rotated, original and translated image
showImage(img)
showImage(rotimg)
showImage(trnsimg)
showImage(trnssimg)

# add circles over an image
showImage(img)
plotRadius(90)
plotRadius(112)
plotRadius(130)

# add ellipses over an image
showImage(trnsimg)
ellipse = ParametricFormEllipse([240.0 220.0], [0.0 0.0], 131)
plotEllipse(ellipse)


text:
module FPGAUtils

function _is_kmod_loaded()
    # Throw an error if the module won't work on the system
    @static if Sys.iswindows() || Sys.isapple()
        error( "FPGAUtils only supports linux systems" )
    else
        isdir( "/sys/class/fpga_manager" ) || error( "fpga_manager kernel module is not loaded" )
    end
end

"""
    is_valid_device( device )

Test to see if `device` is a valid FPGA device loaded by the fpga_mapper kernel module.
"""
function is_valid_device( device )
    _is_kmod_loaded()
    isdir( "/sys/class/fpga_manager/$device" ) || error( "Device $device does not exist" )
end

const kmod_class_dir = "/sys/class/fpga_manager"

include( "constants.jl" )
include( "operations.jl" )

end


text:
## Example of menus with Qt
## add a top-level menu and submenu to which we add
## action, checkable action, radio-actions, separators
## and a popup menu

using PySide

w = Widget()
setWindowTitle(w, "Menu example")
lyt = VBoxLayout(w)
btn = Button(w, "Click for popup")
setLayout(w, lyt); push!(lyt, btn)

## main menubar
mb = MenuBar(w)

## add subment
submenu = Menu("File", mb)
addMenu(mb, submenu)

## add action
act = Action(w, "Action")
change_slot(act, () -> print("action triggered"))
addAction(submenu, act)

## add action with icon
act = Action(w, "Icon")
setIcon(act, Icon("/tmp/julia.gif"))   # nothing on the Mac
change_slot(act, () -> print("icon action"))
addAction(submenu, act)

## add separator
addSeparator(submenu)

## add checkbutton
cb = Action(w, "checkable")
cb[:setCheckable](true)
change_slot(cb, () -> println("checkable ", cb[:isChecked]() ? "is checked" : "is not checked"))
addAction(submenu, cb)

## add separator
addSeparator(submenu)

## add radio group
ag = ActionGroup(w)
ag[:setExclusive](true)
rbs = [Action(w, i) for i in ["Regular", "Medium", "Large"]]
map(rbs) do i
  addAction(ag, i)
  addAction(submenu, i)
  i[:setCheckable](true)
end
rbs[1][:setChecked](true)
qconnect(ag, :triggered, (action) -> println(text(action)))

## do popup on obtn
popup = Menu(btn)
act = Action(w, "fred")
change_slot(act, () -> println("fred"))
addAction(popup, act)

qconnect(btn, :customContextMenuRequested, (pt) ->  popup[:exec_](Qt[:QCursor]()[:pos]()))
btn[:setContextMenuPolicy](qt_enum("CustomContextMenu"))

raise(w)


text:
module Testbyrow

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@transform with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    @test @transform(df, @byrow :n = :i + :g) ≅ @transform(df, :n = :i + :g)
    @test @transform(df, @byrow :n = :t * string(:y)) ≅ @transform(df, :n = :t .* string.(:y))
    @test @transform(df, @byrow :n = :g == 1 ? 100 : 500) ≅ @transform(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @transform(df, @byrow :n = :g == 1 && :t == "a") ≅ @transform(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @transform(df, @byrow :n = first(:g)) ≅ @transform(df, :n = first.(:g))

    d = @transform df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @transform df @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)
    d = @transform df @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @transform df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @transform df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@transform! with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform!, missing])

    @test @transform!(copy(df), @byrow :n = :i + :g) ≅ @transform!(copy(df), :n = :i + :g)
    @test @transform!(copy(df), @byrow :n = :t * string(:y)) ≅ @transform!(copy(df), :n = :t .* string.(:y))
    @test @transform!(copy(df), @byrow :n = :g == 1 ? 100 : 500) ≅ @transform!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @transform!(copy(df), @byrow :n = :g == 1 && :t == "a") ≅ @transform!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @transform!(copy(df), @byrow :n = first(:g)) ≅ @transform!(copy(df), :n = first.(:g))

    d = @transform! df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @transform! df @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @transform! df @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @transform! df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @transform! df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@select with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select, missing])

    @test @select(df, @byrow :n = :i + :g) ≅ @select(df, :n = :i + :g)
    @test @select(df, @byrow :n = :t * string(:y)) ≅ @select(df, :n = :t .* string.(:y))
    @test @select(df, @byrow :n = :g == 1 ? 100 : 500) ≅ @select(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @select(df, @byrow :n = :g == 1 && :t == "a") ≅ @select(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @select(df, @byrow :n = first(:g)) ≅ @select(df, :n = first.(:g))

    d = @select df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @select df @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    d = @select df @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @select df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @select df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@select! with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select!, missing])

    @test @select!(copy(df), @byrow :n = :i + :g) ≅ @select!(copy(df), :n = :i + :g)
    @test @select!(copy(df), @byrow :n = :t * string(:y)) ≅ @select!(copy(df), :n = :t .* string.(:y))
    @test @select!(copy(df), @byrow :n = :g == 1 ? 100 : 500) ≅ @select!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @select!(copy(df), @byrow :n = :g == 1 && :t == "a") ≅ @select!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @select!(copy(df), @byrow :n = first(:g)) ≅ @select!(copy(df), :n = first.(:g))

    d = @select! copy(df) @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @select! copy(df) @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @select! copy(df) @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @select! copy(df) @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @select! copy(df) @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@with with @byrow" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])

    @test @with(df, @byrow :A * 1)   ==  df.A .* 1
    @test @with(df, @byrow :A * :B)  ==  df.A .* df.B

    t = @with df @byrow begin
        :A * 1
    end
    @test t == df.A .* 1

    t = @with df @byrow begin
        :A * :B
    end
    @test t == df.A .* df.B
end

@testset "where with @byrow" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    d = @where df begin
        @byrow :A > 1
        @byrow :B > 1
    end
    @test d ≅ @where(df, :A .> 1, :B .> 1)

    d = @where df @byrow begin
        :A > 1
        :B > 1
    end
    @test d ≅ @where(df, :A .> 1, :B .> 1)
end

@testset "orderby with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    d = @orderby df begin
        @byrow :c
        @byrow :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)

    d = @orderby df @byrow begin
        :c
        :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)
end

@testset "@combine with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :combine, missing])

    gd = groupby(df, :g)

    @test @combine(gd, @byrow :n = :i + :g) ≅ @combine(gd, :n = :i + :g)
    @test @combine(gd, @byrow :n = :t * string(:y)) ≅ @combine(gd, :n = :t .* string.(:y))
    @test @combine(gd, @byrow :n = :g == 1 ? 100 : 500) ≅ @combine(gd, :n = ifelse.(:g .== 1, 100, 500))
    @test @combine(gd, @byrow :n = :g == 1 && :t == "a") ≅ @combine(gd, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @combine(gd, @byrow :n = first(:g)) ≅ @combine(gd, :n = first.(:g))

    d = @combine gd @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @combine(gd, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @combine gd @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)
    d = @combine gd @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)

    d = @combine gd @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)

    d = @combine gd @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@by with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :combine, missing])

    @test @by(df, :g, @byrow :n = :i + :g) ≅ @by(df, :g, :n = :i + :g)
    @test @by(df, :g, @byrow :n = :t * string(:y)) ≅ @by(df, :g, :n = :t .* string.(:y))
    @test @by(df, :g, @byrow :n = :g == 1 ? 100 : 500) ≅ @by(df, :g, :n = ifelse.(:g .== 1, 100, 500))
    @test @by(df, :g, @byrow :n = :g == 1 && :t == "a") ≅ @by(df, :g, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @by(df, :g, @byrow :n = first(:g)) ≅ @by(df, :g, :n = first.(:g))

    d = @by df :g @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @by(df, :g, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @by df :g @byrow begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i .* :g, :n2 = :i .* :g)
end

end


text:
#!/usr/bin/env julia

# test polygon simplification

using Luxor, Colors

function circles(points)
    gsave()
    sethue(1.0,0,0)
    for p in points
        circle(p.x, p.y, .4, :fill) # show vertex locations
    end
    grestore()
end

function sinecurves()
    gsave()
    translate(100, 100)
    setline(0.25)
    sethue(0,0,0)
    x_vals = collect(0:pi/180: 20pi)
    # generate array of Points
    plist = [Point(d * 2pi , -sin(d) * cos(12 * d) * 8 * sin(d/10)) for d in x_vals]
    poly(plist, :stroke) # draw original
    circles(plist)
    text("original " * string(length(plist)) * " vertices" , 0, -15)
    for detail in [0.01, 0.05, 0.1, 0.2, 0.5, 0.75, 1.0, 2.0]
        translate(0, 50)
        simplified = simplify(plist, detail)
        poly(simplified, :stroke)
        circles(simplified)
        text("detail  $(detail), " * string(length(simplified)) * " vertices" , 0, -15)
    end
    grestore()
end

function test(pagewidth, pageheight)
    gsave()
    translate(100, pageheight/2)
    setline(0.25)
    sethue(0,0,0)
    x_vals = collect(0:pi/100: 4pi)
    polyline = [Point(d * 20pi, 5 * -sin(d) * cos(12 * d) * 8 * sin(d/10)) for d in x_vals]
    poly(polyline, :stroke)
    circles(polyline)
    text("original " * string(length(polyline)) * " vertices" , 0, -30)
    for detail in [0.01, 0.05, 0.075, 0.1, 0.2, 0.5, 0.75, 1.0, 2.0]
        translate(0, 60)
        polysimple = simplify(polyline, detail)
        poly(polysimple, :stroke)
        circles(polysimple)
        text(" detail ($detail), " * string(length(polysimple)) * " vertices" , 0, -30)
    end
    grestore()
end

function simplify_poly(fname)
    pagewidth  = 1190.0 # points
    pageheight = 1684.0 # points
    Drawing(pagewidth, pageheight, fname)
    sinecurves()
    test(pagewidth, pageheight)
    finish()
    println("finished test: output in $(fname)")
end

@time simplify_poly("/tmp/simplify-poly.pdf")


text:
# Additive Kernel References

composition_classes = (
    ExponentialClass,
    GammaExponentialClass,
    RationalClass,
    GammaRationalClass,
    MaternClass,
    ExponentiatedClass,
    PolynomialClass,
    PowerClass,
    LogClass,
    GammaLogClass,
    SigmoidClass
)

additive_kernels = (
    SquaredDistanceKernel,
    SineSquaredKernel,
    ChiSquaredKernel,
    ScalarProductKernel
)

composition_pairs = Dict(
    ExponentialClass      => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    GammaExponentialClass => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    RationalClass         => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    GammaRationalClass    => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    MaternClass           => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    ExponentiatedClass    => Dict(ScalarProductKernel=>true),
    PolynomialClass       => Dict(ScalarProductKernel=>true),
    PowerClass            => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    LogClass              => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    GammaLogClass         => Dict(SquaredDistanceKernel=>true, SineSquaredKernel=>true,
                                  ChiSquaredKernel=>true),
    SigmoidClass          => Dict(ScalarProductKernel=>true)
)

all_default_args = Dict(
    SquaredDistanceKernel => ([], Int[]),
    SineSquaredKernel     => ([π], Int[]),
    ChiSquaredKernel      => ([], Int[]),
    ScalarProductKernel   => ([], Int[]),
    ExponentialClass      => ([1,], Int[]),
    GammaExponentialClass => ([1,0.5],  Int[]),
    RationalClass         => ([1,1], Int[]),
    GammaRationalClass    => ([1,1,0.5], Int[]),
    MaternClass           => ([1,1], Int[]),
    ExponentiatedClass    => ([1,0], Int[]),
    PolynomialClass       => ([1,1], Int[3]),
    PowerClass            => ([1,0,1], Int[]),
    LogClass              => ([1,], Int[]),
    GammaLogClass         => ([1,0.5], Int[]),
    SigmoidClass          => ([1,0], Int[])
)

all_phifunctions = Dict(
    SquaredDistanceKernel => (x,y)     -> ((x-y)^2),
    SineSquaredKernel     => (p,x,y)   -> (sin(p*(x-y))^2),
    ChiSquaredKernel      => (x,y)     -> (x == y == 0) ? zero(typeof(x)) : ((x-y)^2/(x+y)),
    ScalarProductKernel   => (x,y)     -> x*y,
    ExponentialClass      => (α,z)     -> exp(-α*z^γ),
    GammaExponentialClass => (α,γ,z)   -> exp(-α*z^γ),
    RationalClass         => (α,β,z)   -> (1 + α*z)^(-β),
    GammaRationalClass    => (α,β,γ,z) -> (1 + α*z^γ)^(-β),
    MaternClass           => (ν,θ,z)   -> begin
                                              v1 = sqrt(2*ν) * z / θ
                                              v1 = v1 < eps(typeof(z)) ? eps(typeof(z)) : v1
                                              2*(v1/2)^ν * besselk(ν,v1)/gamma(ν)
                                          end,
    ExponentiatedClass    => (a,c,z)   -> exp(a*z+c),
    PolynomialClass       => (a,c,d,z) -> (a*z+c)^d,
    PowerClass            => (a,c,γ,z) -> (a*z+c)^γ,
    LogClass              => (α,z)     -> log(α*z+1),
    GammaLogClass         => (α,γ,z)   -> log(α*z^γ+1),
    SigmoidClass          => (a,c,z)   -> tanh(a*z+c)

)

all_kernelproperties = Dict( #|mercer |negdef |neg    |zero   |pos
    SineSquaredKernel =>      (false,  true,   false,  true,   true),
    ExponentiatedClass =>     (true,   false,  false,  false,  true),
    ScalarProductKernel =>    (true,   false,  true,   true,   true),
    SquaredDistanceKernel =>  (false,  true,   false,  true,   true),
    SigmoidClass =>           (false,  false,  true,   true,   true),
    ExponentialClass =>       (true,   false,  false,  false,  true),
    GammaLogClass =>          (false,  true,   false,  true,   true),
    GammaRationalClass =>     (true,   false,  false,  false,  true),
    LogClass =>               (false,  true,   false,  true,   true),
    MaternClass =>            (true,   false,  false,  false,  true),
    PolynomialClass =>        (true,   false,  true,   true,   true),
    PowerClass =>             (false,  true,   false,  true,   true),
    GammaExponentialClass =>  (true,   false,  false,  false,  true),
    ChiSquaredKernel =>       (false,  true,   false,  true,   true),
    RationalClass =>          (true,   false,  false,  false,  true)
 )

composition_kernels = (
    GaussianKernel,
    SquaredExponentialKernel,
    RadialBasisKernel,
    LaplacianKernel,
    PeriodicKernel,
    RationalQuadraticKernel,
    MaternKernel,
    MatérnKernel,
    PolynomialKernel,
    LinearKernel,
    SigmoidKernel
)

#=
all_test_args = Dict(
    SquaredDistanceKernel  => ([:t],[1],[0.5]),
    SineSquaredKernel      => ([:p, :t],[π,1],[2π,0.5]),
    ChiSquaredKernel       => ([:t],[1],[0.5]),
    ScalarProductKernel    => (Symbol[], Int[], Int[]),
    GammaExponentialClass  => ([:alpha, :gamma], [1,0.5], [2,1]),
    ExponentialClass       => ([:alpha], [1,], [2,]),
    GammaRationalClass     => ([:alpha, :beta, :gamma], [1,1,0.5], [2,2,1]),
    RationalClass          => ([:alpha, :beta], [1,1], [2,2]),
    MaternClass            => ([:nu, :theta], [1,1], [2,2]),
    ExponentiatedClass     => ([:a, :c], [1,0], [2,1]),
    PolynomialClass        => ([:a, :c, :d], [1,1,2], [2,2,3]),
    PowerClass             => ([:a, :c, :gamma], [1,0,1], [2,1,0.5]),
    GammaLogClass          => ([:alpha, :gamma], [1,0.5], [2,1]),
    LogClass               => ([:alpha], [1,], [2,]),
    SigmoidClass           => ([:a, :c], [1, 0], [2, 1])
)

all_testinputs = Dict(
    SquaredDistanceKernel  => ([1,1],[1,0],[0,1],[0,0],[-1,-1],[-1,0],[0,-1]),
    SineSquaredKernel      => ([1,1],[1,0],[0,1],[0,0],[-1,-1],[-1,0],[0,-1]),
    ChiSquaredKernel       => ([1,1],[1,0],[0,1],[0,0],[2,1],[1,2],[2,0],[0,2]),
    ScalarProductKernel    => ([1,1],[1,0],[0,1],[0,0],[-1,-1],[-1,0],[0,-1]),
    ExponentialClass       => ([0],[1e-10],[0.5],[1]),
    RationalClass => ([0],[1e-10],[0.5],[1]),
    MaternClass            => ([0],[1e-10],[0.5],[1]),
    PowerClass             => ([0],[1e-10],[0.5],[1]),
    LogClass               => ([0],[1e-10],[0.5],[1]),
    PolynomialClass        => ([0],[1e-10],[0.5],[1],[-1e-10],[-0.5],[-1]),
    ExponentiatedClass     => ([0],[1e-10],[0.5],[1],[-1e-10],[-0.5],[-1]),
    SigmoidClass           => ([0],[1e-10],[0.5],[1],[-1e-10],[-0.5],[-1])

)

all_phicases = Dict(
    SquaredDistanceKernel  => ([0.25],[0.5],[1]),
    SineSquaredKernel      => ([π,0.25],[π,0.5],[π,1]),
    ChiSquaredKernel       => ([0.25],[0.5],[1]),
    ScalarProductKernel    => (Int[],),
    ExponentialClass       => ([1,1], [1,0.5]),
    RationalClass => ([1,1,1], [1,1,0.5], [1,2,1], [2,2,0.5]),
    MaternClass            => ([1,1], [2,1]),
    PowerClass             => ([1,0,0.5], [1,0,1]),
    LogClass               => ([1,0.5], [1,1]),
    PolynomialClass        => ([1,1,2],[1,1,1]),
    ExponentiatedClass     => ([1,0],),
    SigmoidClass           => ([1,1],)
)

all_errorcases = Dict(
    SquaredDistanceKernel  => ([1.1], [0]),
    SineSquaredKernel      => ([π,1.1], [π,0], [0,1]),
    ChiSquaredKernel       => ([1.1], [0]),
    ExponentialClass       => ([0,1], [1,1.1], [1,0]),
    RationalClass => ([0,1,1], [1,0,1], [1,1,0], [1,1,1.1]),
    MaternClass            => ([1,0], [0,1]),
    PowerClass             => ([0,1,1], [0,-1,1], [1,0,0], [1,0,1.1]),
    LogClass               => ([0,1], [1,0], [1,1.1]),
    PolynomialClass        => ([0,1,1],[1,-0.1,1],[1,1,1.5],[1,1,0]),
    ExponentiatedClass     => ([0],),
    SigmoidClass           => ([0,1],[1,-0.1])
)

=#