text:
# Date internal algorithms
@test Dates.totaldays(0,2,28) == -307
@test Dates.totaldays(0,2,29) == -306
@test Dates.totaldays(0,3,1) == -305
@test Dates.totaldays(0,12,31) == 0
# Rata Die Days # start from 0001-01-01
@test Dates.totaldays(1,1,1) == 1
@test Dates.totaldays(1,1,2) == 2
@test Dates.totaldays(2013,1,1) == 734869

@test Dates.daysinmonth(2000,1) == 31
@test Dates.daysinmonth(2000,2) == 29
@test Dates.daysinmonth(2000,3) == 31
@test Dates.daysinmonth(2000,4) == 30
@test Dates.daysinmonth(2000,5) == 31
@test Dates.daysinmonth(2000,6) == 30
@test Dates.daysinmonth(2000,7) == 31
@test Dates.daysinmonth(2000,8) == 31
@test Dates.daysinmonth(2000,9) == 30
@test Dates.daysinmonth(2000,10) == 31
@test Dates.daysinmonth(2000,11) == 30
@test Dates.daysinmonth(2000,12) == 31
@test Dates.daysinmonth(2001,2) == 28

@test Dates.isleapyear(1900) == false
@test Dates.isleapyear(2000) == true
@test Dates.isleapyear(2004) == true
@test Dates.isleapyear(2008) == true
@test Dates.isleapyear(0) == true
@test Dates.isleapyear(1) == false
@test Dates.isleapyear(-1) == false
@test Dates.isleapyear(4) == true
@test Dates.isleapyear(-4) == true

# Create "test" check manually
test = Dates.DateTime(Dates.UTM(63492681600000))
# Test DateTime construction by parts
@test Dates.DateTime(2013) == test
@test Dates.DateTime(2013,1) == test
@test Dates.DateTime(2013,1,1) == test
@test Dates.DateTime(2013,1,1,0) == test
@test Dates.DateTime(2013,1,1,0,0) == test
@test Dates.DateTime(2013,1,1,0,0,0) == test
@test Dates.DateTime(2013,1,1,0,0,0,0) == test
test = Dates.Date(Dates.UTD(734869))
# Test Date construction by parts
@test Dates.Date(2013) == test
@test Dates.Date(2013,1) == test
@test Dates.Date(2013,1,1) == test

# Test various input types for Date/DateTime
test = Dates.Date(1,1,1)
@test Dates.Date(int8(1),int8(1),int8(1)) == test
@test Dates.Date(uint8(1),uint8(1),uint8(1)) == test
@test Dates.Date(int16(1),int16(1),int16(1)) == test
@test Dates.Date(uint8(1),uint8(1),uint8(1)) == test
@test Dates.Date(int32(1),int32(1),int32(1)) == test
@test Dates.Date(uint32(1),uint32(1),uint32(1)) == test
@test Dates.Date(int64(1),int64(1),int64(1)) == test
@test Dates.Date('\x01','\x01','\x01') == test
@test Dates.Date(true,true,true) == test
@test_throws ArgumentError Dates.Date(false,true,false)
@test Dates.Date(false,true,true) == test - Dates.Year(1)
@test_throws ArgumentError Dates.Date(true,true,false)
@test Dates.Date(uint64(1),uint64(1),uint64(1)) == test
@test Dates.Date(0xffffffffffffffff,uint64(1),uint64(1)) == test - Dates.Year(2)
@test Dates.Date(int128(1),int128(1),int128(1)) == test
@test Dates.Date(170141183460469231731687303715884105727,int128(1),int128(1)) == test - Dates.Year(2)
@test Dates.Date(uint128(1),uint128(1),uint128(1)) == test
@test Dates.Date(big(1),big(1),big(1)) == test
@test Dates.Date(big(1),big(1),big(1)) == test
# Potentially won't work if can't losslessly convert to Int64
@test Dates.Date(BigFloat(1),BigFloat(1),BigFloat(1)) == test
@test Dates.Date(complex(1),complex(1),complex(1)) == test
@test Dates.Date(float64(1),float64(1),float64(1)) == test
@test Dates.Date(float32(1),float32(1),float32(1)) == test
@test Dates.Date(float16(1),float16(1),float16(1)) == test
@test Dates.Date(Rational(1),Rational(1),Rational(1)) == test
@test_throws InexactError Dates.Date(BigFloat(1.2),BigFloat(1),BigFloat(1))
@test_throws InexactError Dates.Date(1 + im,complex(1),complex(1))
@test_throws InexactError Dates.Date(1.2,1.0,1.0)
@test_throws InexactError Dates.Date(1.2f0,1.f0,1.f0)
@test_throws InexactError Dates.Date(3//4,Rational(1),Rational(1)) == test

# Months, days, hours, minutes, seconds, and milliseconds must be in range
@test_throws ArgumentError Dates.Date(2013,0,1)
@test_throws ArgumentError Dates.Date(2013,13,1)
@test_throws ArgumentError Dates.Date(2013,1,0)
@test_throws ArgumentError Dates.Date(2013,1,32)
@test_throws ArgumentError Dates.DateTime(2013,0,1)
@test_throws ArgumentError Dates.DateTime(2013,13,1)
@test_throws ArgumentError Dates.DateTime(2013,1,0)
@test_throws ArgumentError Dates.DateTime(2013,1,32)
@test_throws ArgumentError Dates.DateTime(2013,1,1,24)
@test_throws ArgumentError Dates.DateTime(2013,1,1,-1)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,-1)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,60)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,0,-1)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,0,60)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,0,0,-1)
@test_throws ArgumentError Dates.DateTime(2013,1,1,0,0,0,1000)

# Test DateTime traits
a = Dates.DateTime(2000)
b = Dates.Date(2000)
@test Dates.calendar(a) == Dates.ISOCalendar
@test Dates.calendar(b) == Dates.ISOCalendar
@test Dates.precision(a) == Dates.Millisecond
@test Dates.precision(b) == Dates.Day
@test string(typemax(Dates.DateTime)) == "146138512-12-31T23:59:59"
@test string(typemin(Dates.DateTime)) == "-146138511-01-01T00:00:00"
@test typemax(Dates.DateTime) - typemin(Dates.DateTime) == Dates.Millisecond(9223372017043199000)
@test string(typemax(Dates.Date)) == "252522163911149-12-31"
@test string(typemin(Dates.Date)) == "-252522163911150-01-01"
@test convert(Real,b) == 730120
@test convert(Float64,b) == 730120.0
@test convert(Int32,b) == 730120
@test convert(Real,a) == 63082368000000
@test convert(Float64,a) == 63082368000000.0
@test convert(Int64,a) == 63082368000000
# Date-DateTime conversion/promotion
@test Dates.DateTime(a) == a
@test Dates.Date(a) == b
@test Dates.DateTime(b) == a
@test Dates.Date(b) == b
@test a == b
@test a == a
@test b == a
@test b == b
@test !(a < b)
@test !(b < a)
c = Dates.DateTime(2000)
d = Dates.Date(2000)
@test ==(a,c)
@test ==(c,a)
@test ==(d,b)
@test ==(b,d)
@test ==(a,d)
@test ==(d,a)
@test ==(b,c)
@test ==(c,b)
b = Dates.Date(2001)
@test b > a
@test a < b
@test a != b
@test Dates.Date(Dates.DateTime(Dates.Date(2012,7,1))) == Dates.Date(2012,7,1)

y = Dates.Year(1)
m = Dates.Month(1)
w = Dates.Week(1)
d = Dates.Day(1)
h = Dates.Hour(1)
mi = Dates.Minute(1)
s = Dates.Second(1)
ms = Dates.Millisecond(1)
@test Dates.DateTime(y) == Dates.DateTime(1)
@test Dates.DateTime(y,m) == Dates.DateTime(1,1)
@test Dates.DateTime(y,m,d) == Dates.DateTime(1,1,1)
@test Dates.DateTime(y,m,d,h) == Dates.DateTime(1,1,1,1)
@test Dates.DateTime(y,m,d,h,mi) == Dates.DateTime(1,1,1,1,1)
@test Dates.DateTime(y,m,d,h,mi,s) == Dates.DateTime(1,1,1,1,1,1)
@test Dates.DateTime(y,m,d,h,mi,s,ms) == Dates.DateTime(1,1,1,1,1,1,1)
@test Dates.DateTime(Dates.Day(10),Dates.Month(2),y) == Dates.DateTime(1,2,10)
@test Dates.DateTime(Dates.Second(10),Dates.Month(2),y,Dates.Hour(4)) == Dates.DateTime(1,2,1,4,0,10)
@test Dates.DateTime(Dates.Year(1),Dates.Month(2),Dates.Day(1),
                     Dates.Hour(4),Dates.Second(10)) == Dates.DateTime(1,2,1,4,0,10)
@test Dates.Date(y) == Dates.Date(1)
@test Dates.Date(y,m) == Dates.Date(1,1)
@test Dates.Date(y,m,d) == Dates.Date(1,1,1)
@test Dates.Date(m) == Dates.Date(1,1,1)
@test Dates.Date(d,y) == Dates.Date(1,1,1)
@test Dates.Date(d,m) == Dates.Date(1,1,1)
@test Dates.Date(m,y) == Dates.Date(1,1,1)


text:
function calculate_bisector_span(λrest::T, wav::AA{T,1}) where T<:AF
    minw = minimum(filter(!isnan, wav))
    return abs(minw - wav[1])/wav[1] * (c /100.0)
end

function calculate_bisector_span(λrest::T, wav::AA{T,2}) where T<:AF
    out = zeros(size(wav,2))
    for i in eachindex(out)
        out[i] = calculate_bisector_span(λrest, wav[:,i])
    end
    return out
end

function calculate_bisector_extreme(λrest::T, wav::AA{T,2}, bis::AA{T,2}) where T<:AF
    out = zeros(size(wav,2))
    for i in eachindex(out)
        out[i] = calculate_bisector_extreme(λrest, wav[:,i], bis[:,i])
    end
    return out
end

function calculate_bisector_extreme(λrest::T, wav::AA{T,1}, bis::AA{T,1}) where T<:AF
    ind1 = searchsortednearest(wav, 0.4)
    ind2 = searchsortednearest(wav, 0.8)
    wpos = mean(wav[ind1:ind2])
    return (λrest - wpos)/λrest * c/100.0
end

function calculate_bisector_lslope(λrest::T, wav::AA{T,1}, bis::AA{T,1}) where T<:AF
    dλ = (minimum(wav[.!isnan.(wav)]) - wav[1])/λrest * c/100.0
    ind1 = findfirst(wav .== minimum(wav[.!isnan.(wav)]))

    dF = bis[ind1] - bis[1]
    return dF/dλ
end

function calculate_bisector_lslope(λrest::T, wav::AA{T,2}, bis::AA{T,2}) where T<:AF
    out = zeros(size(wav,2))
    for i in eachindex(out)
        out[i] = calculate_bisector_lslope(λrest, wav[:,i], bis[:,i])
    end
    return out
end

function measure_bisector(xs::AA{T,1}, ys::AA{T,1}; interpolate::Bool=true,
                          top::T=0.99, len::Integer=100) where T<:AF
    if interpolate
        return measure_bisector_interpolate(xs, ys, top=top, len=len)
    else
        return measure_bisector_loop(xs, ys, top=top, len=len)
    end
end

function measure_bisector_interpolate(xs::AA{T,1}, ys::AA{T,1}; top::T=0.99,
                                      len::Integer=100, max_loop::Int=20) where T<:AF
    # check lengths and normalization
    @assert length(xs) == length(ys)

    # normalize the spec, find bottom of line
    ys ./= maximum(ys)
    botind = argmin(ys)
    depths = range(ys[botind], top, length=len)

    # find left and right halves
    lind = findfirst(ys .< top)
    rind = findlast(ys .< top)
    lspec = reverse(ys[lind:botind])
    rspec = ys[botind:rind]
    lwav = reverse(xs[lind:botind])
    rwav = xs[botind:rind]

    # make sure lspec is sorted
    num_loop = 0
    unsorted = !issorted(lspec)
    while unsorted
        num_loop += 1
        ldiff = diff(lspec)
        inds = BitArray(vcat(0, ldiff .< 0))
        lspec[inds] .= NaN
        itp = LinearInterpolation(reverse(lwav[.!inds]),
                                  reverse(lspec[.!inds]),
                                  extrapolation_bc=Flat())
        lspec .= itp.(lwav)
        unsorted = !issorted(lspec)
        if num_loop > max_loop
            break
        end
    end

    # make sure rspec is sorted
    num_loop = 0
    unsorted = !issorted(rspec)
    while unsorted
        num_loop += 1
        rdiff = diff(rspec)
        inds = BitArray(vcat(0, rdiff .< 0))
        rspec[inds] .= NaN
        itp = LinearInterpolation(rwav[.!inds], rspec[.!inds],
                                  extrapolation_bc=Flat())
        rspec .= itp.(rwav)
        unsorted = !issorted(rspec)
        if num_loop > max_loop
            break
        end
    end

    # interpolate wavelengths onto intensity grid
    lspline = LinearInterpolation(lspec, lwav, extrapolation_bc=Flat())
    rspline = LinearInterpolation(rspec, rwav, extrapolation_bc=Flat())
    wavs = (lspline(depths) .+ rspline(depths)) ./ 2.0
    return wavs, depths
end

function measure_bisector_loop(xs::AA{T,1}, ys::AA{T,1}; top::T=0.99,
                               len::Integer=100) where T<:AF
    # normalize the spec, find bottom of line
    ys ./= maximum(ys)

    # assign depths to measure bisector at
    dep = range(one(T)-minimum(ys)-0.01, one(T) - top, length=len)

    # set iterators
    nccf = Int(length(xs) ÷ 2)
    L = nccf
    R = nccf

    # allocate memory
    xL = zeros(len)
    xR = zeros(len)
    wav = zeros(len)

    # loop over depths
    for d in eachindex(dep)
        y = one(T) - dep[d]
        while((ys[L] < y) & (L > 0))
            L -= 1
        end

        while ((ys[R] < y) & (R < length(xs)))
            R += 1
        end

        if ((y > maximum(ys[1:nccf])) | (y > maximum(ys[nccf+1:end])))
            L = 0
            R = length(xs)
        end

        if L == 0
            xL[d] = xL[d-1]
        else
            mL = (xs[L+1] - xs[L]) / (ys[L+1] - ys[L])
            xL[d] = xs[L] + mL * (y - ys[L])
        end

        if R == length(xs)
            xR[d] = xR[d-1]
        else
            mR = (xs[R-1] - xs[R]) / (ys[R-1] - ys[R])
            xR[d] = xs[R] + mR * (y - ys[R])
        end
        wav[d] = (xL[d] + xR[d]) / 2.0
    end
    return wav, one(T) .- dep
end

function bisector_uncertainty(wav::AA{T,1}, bis::AA{T,1}) where T<:AF
    dF = diff(bis)
    dv = diff(wav)
    return one(T)/sqrt(2.0) .* dF ./ abs.(dF ./ dv)
end

ext: jl

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using StochasticDiffEq, Plots, DiffEqDevTools, DiffEqProblemLibrary
using ParameterizedFunctions, DiffEqMonteCarlo, Base.Test, OrdinaryDiffEq
using Plots; pyplot()
using BenchmarkTools, DiffEqNoiseProcess

################################################################################

### Van Der Pol

van = @ode_def VanDerPol2 begin
  dy = μ*((1-x^2)*y - x)
  dx = 1*y
end μ

σ = (du,u,p,t) -> begin
  for i in 1:2
    du[i] = 3.0 #Additive
  end
end

ode_prob = ODEProblem(van,[0;2.],(0.0,6.3),(1e5,))
prob = SDEProblem(van,σ,[0;2.],(0.0,6.3),(1e5,))
monte_prob = MonteCarloProblem(prob)
N = 100; fails = zeros(5,4); times = zeros(5,4)
offset = 0

#####

@time sol0 =solve(ode_prob,Tsit5();abstol=1e-6,reltol=1e-3)
p0 = plot(sol0,plotdensity=20000,denseplot=true,ylims=[-10,10],title="ODE Solution",
          xtickfont = font(16, "Latin Modern Roman"),
          titlefont = font(20, "Latin Modern Roman"),
          ytickfont = font(16, "Latin Modern Roman"),
          guidefont = font(18, "Latin Modern Roman"),
          legendfont = font(14, "Latin Modern Roman"))

srand(101)
test_dt = 0.0001
t = prob.tspan[1]:test_dt:prob.tspan[2]
brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
np = NoiseGrid(t,brownian_values,brownian_values2)
fixed_noise_prob = SDEProblem(van,σ,[0;2.],(0.0,6.3),(1e5,),noise=np)

@time sol1 =solve(fixed_noise_prob,SOSRA();abstol=1,reltol=1/2^1)
p1 = plot(sol1,plotdensity=20000,denseplot=true,ylims=[-10,10],
          title="High tolerance SOSRA",
          legend = false,
          xtickfont = font(16, "Latin Modern Roman"),
          titlefont = font(20, "Latin Modern Roman"),
          ytickfont = font(16, "Latin Modern Roman"),
          guidefont = font(18, "Latin Modern Roman"),
          legendfont = font(14, "Latin Modern Roman"))

#@time sol1 =solve(prob,ImplicitEM();dt=1/1000000) # Errors

@time sol2 =solve(fixed_noise_prob,SOSRA2();abstol=1/2,reltol=1/2^2)
plot(sol2,plotdensity=20000,denseplot=true,ylims=[-10,10])
#@time sol3 =solve(prob,SRA1();abstol=1/2^7,reltol=1/2^5)
@time sol4 =solve(fixed_noise_prob,SRA3();abstol=1/2^3,reltol=1/2^3)
p2 = plot(sol4,plotdensity=20000,denseplot=true,ylims=[-10,10],
          title="High Tolerance SRA3",
          legend = false,
          xtickfont = font(16, "Latin Modern Roman"),
          titlefont = font(20, "Latin Modern Roman"),
          ytickfont = font(16, "Latin Modern Roman"),
          guidefont = font(18, "Latin Modern Roman"),
          legendfont = font(14, "Latin Modern Roman"))
@time sol42 =solve(fixed_noise_prob,SRA3();abstol=1/2^6,reltol=1/2^4)
p3 = plot(sol42,plotdensity=20000,denseplot=true,ylims=[-10,10],
          title="Low Tolerance SRA3",
          legend = false,
          xtickfont = font(16, "Latin Modern Roman"),
          titlefont = font(20, "Latin Modern Roman"),
          ytickfont = font(16, "Latin Modern Roman"),
          guidefont = font(18, "Latin Modern Roman"),
          legendfont = font(14, "Latin Modern Roman"))

@time sol52 =solve(fixed_noise_prob,SKenCarp();abstol=16,reltol=16)
p5 = plot(sol52,plotdensity=20000,denseplot=true,ylims=[-10,10])

diff(sol52.t)

plot(p0,p1,p2,p3,layout=grid(2,2),size=(1200,800))

savefig("additive_van_der_pol.png")
savefig("additive_van_der_pol.pdf")

######

#### Setup
@time sol1 =solve(prob,SOSRA();abstol=10,reltol=1/2^1)
@time sol2 =solve(prob,SOSRA2();abstol=1/2,reltol=1/2^2)
@time sol3 =solve(prob,SRA1();abstol=1/2^7,reltol=1/2^5)
@time sol4 =solve(prob,SRA3();abstol=1/2^6,reltol=1/2^3)
@time sol5 =solve(prob,SKenCarp();abstol=1,reltol=1)

#### Estimates

prob_quick = SDEProblem(van,σ,[0;2.],(0.0,1.0),(1e5,))
@time sol6 =solve(prob_quick,EM();dt=5e-7) # Unstable
@time sol6 =solve(prob_quick,EM();dt=1e-7) # Stable?

42.087122 *6.3 * 100 # 7 hours 21 minutes....

@time sol6 =solve(prob,EM();dt=1e-7) # Stable?

346.117152 * 100 # 9 hours 36 minutes 52 seconds

@time sol6 =solve(prob,EM();dt=5e-8) # Stable?

@time sol6 =solve(prob_quick,ImplicitEM();dt=1e-7) # Unstable
@time sol6 =solve(prob_quick,ImplicitEM();dt=5e-8) # Stable?

363.668494 * 6.3 * 100 # 2 days 15 hours 38 minutes and 31 seconds...

#### Monte Carlo
println("SOSRA time")
@time sol1 =solve(monte_prob,SOSRA();num_monte=N,abstol=10,reltol=1/2^1,verbose=false,save_everystep=false)
println("SOSRA2 time")
@time sol2 =solve(monte_prob,SOSRA2();num_monte=N,abstol=1,reltol=1/2^2,verbose=false,save_everystep=false)
println("SRA1 time")
@time sol3 =solve(monte_prob,SRA1();num_monte=N,abstol=1/2^6,reltol=1/2^(5),verbose=false,save_everystep=false)
println("SRA3 time")
@time sol4 =solve(monte_prob,SRA3();num_monte=N,abstol=1/2^6,reltol=1/2^3,verbose=false,save_everystep=false)
println("SKenCarp time")
@time sol5 =solve(monte_prob,SKenCarp();num_monte=N,abstol=16,reltol=16,verbose=false,save_everystep=false)

sum(sol.retcode == :Success for sol in sol5)

@time sol5 =solve(monte_prob,SKenCarp();num_monte=N,abstol=4,reltol=4,verbose=false,save_everystep=false)

sum(sol.retcode == :Success for sol in sol5)

@time sol5 =solve(monte_prob,SKenCarp(extrapolant=:trivial);num_monte=N,
                  abstol=4,reltol=4,verbose=false,save_everystep=false)

sum(sol.retcode == :Success for sol in sol5)

println("EM time")
# dt = 1e-7 has 15 failures
@time sol5 =solve(monte_prob,EM();num_monte=N,dt=1e-8,save_everystep=false)

sum(sol.retcode == :Success for sol in sol5)

#######

# Eigenvalue estimates (put inside of the perform_step)

integrator.f(@muladd(t + c₀[2]*dt),H0[2],ftmp)
integrator.f(@muladd(t + c₀[3]*dt),H0[3],gtmp)
eig_est = dt*(norm(ftmp-gtmp)/norm(H0[2]-H0[3]))/5
eig_est > 5 && println("$t")


text:
module OpenVR
  const openvr_version = "1.2.10"
  const openvr_dirname = "openvr-$openvr_version"
  const module_dir = dirname(dirname(@__FILE__))
  const src_dir = joinpath(module_dir,"src")
  const openvr_dir = joinpath(module_dir,"deps",openvr_dirname)
  const os_dir_and_so_ext = begin
    local os_dirname = "unknown"
    local os_so_extension = ""
    Sys.islinux()   &&                        (os_dirname = "linux64"; os_so_extension=".so")
    Sys.isapple()   &&                        (os_dirname = "osx32"  ; os_so_extension=".dylib")
    Sys.iswindows() && Sys.ARCH == :x86_32 && (os_dirname = "win32"  ; os_so_extension=".dll")
    Sys.iswindows() && Sys.ARCH == :x86_64 && (os_dirname = "win64"  ; os_so_extension=".dll")
    (os_dirname, os_so_extension)
  end
  const openvr_bin_dir = joinpath(openvr_dir,"lib",os_dir_and_so_ext[1])
  const stage1_path = joinpath(openvr_bin_dir,"libvrstage1$(os_dir_and_so_ext[2])")
  const stage2_path = joinpath(openvr_bin_dir,"libvrstage2$(os_dir_and_so_ext[2])")

  using ModernGL
  #  add Cxx#master

  module VR_enums
    using CxxWrap
    # @wrapmodule(joinpath("/home/christianl/src/openvr/samples/bin/linux64","libhellovr_preamble"))

    # https://stackoverflow.com/questions/32518131/reference-function-of-outer-module
    import ...OpenVR: stage1_path
    @wrapmodule(stage1_path)

    function __init__()
      @initcxx
    end
  end

  for n in names(VR_enums; all=true, imported=false) # skip `#eval` and `#include`
    startswith(string(n),"__cxxwrap") && continue
    # skip `#eval` and `#include` etc.
    n ∈ [:eval, :include, :__init__, :VR_enums, Symbol("#include"), Symbol("#eval"), Symbol("#__init__")] && continue
    Core.eval(@__MODULE__,:(const $n = VR_enums.$n))
  end

  const ETrackedPropertyErrorRef = Ref{ETrackedPropertyError}

  const UInt8Ptr = Ptr{UInt8}

  # const VR0 = VR_enums
  # include("openvr/hellovr_opengl_julia/vrstructs.jl")
  include(joinpath(src_dir,"vrstructs.jl"))
  # include("openvr/hellovr_opengl_julia/HelloVRstructs.jl")

  # create Ref{T} and Ptr{T}, since we cannot use these names directly from the C++ shared library
  for n in names(@__MODULE__; all=true, imported=false)
    T = getfield(@__MODULE__,n)
    # only DataTypes
    typeof(T) != DataType && continue
    # only those which were defined in this module (by checking their prefix)
    Tname = string(T.name)
    Mname = occursin(".",Tname) ? split(Tname,".")[end-1] : ""
    (Mname == String(nameof(@__MODULE__)) || Mname == String(nameof(VR_enums))) || continue
    # skip `#eval` and `#include`
    n ∈ [Symbol("#include"), Symbol("#eval")] && continue
    # build expression
    refname = Symbol(String(n)*"Ref")
    # ptrname = Symbol(String(n)*"Ptr")
    refexpr = :(const $refname = Ref{$T})
    # ptrexpr = :(const $ptrname = Ptr{$T})
    Core.eval(@__MODULE__,refexpr)
    # display(n)
    # Core.eval(@__MODULE__,ptrexpr)
    # TODO: check all field-types on their bitstypeness
    #   if any, then the layout is guaranteed to not match
    ~isbitstype(T) && println("WARNING $(nameof(T)) is NOT a bitstype anymore. The memory layout might not match what C/C++ expects.")
  end

  const RenderModel_tPtrRef = Ref{Ptr{RenderModel_t}} # TODO: figure out how T** template should work
  const RenderModel_TextureMap_tPtrRef = Ref{Ptr{RenderModel_TextureMap_t}}

  struct FramebufferDesc
    m_nDepthBufferId :: GLuint
    m_nRenderTextureId :: GLuint
    m_nRenderFramebufferId :: GLuint
    m_nResolveTextureId :: GLuint
    m_nResolveFramebufferId :: GLuint
  end

  using CxxWrap
  # @wrapmodule(joinpath("/home/christianl/src/openvr/samples/bin/linux64","libhellovr_julia"))
  # CxxWrap.wrapmodule(joinpath("/home/christianl/src/openvr/samples/bin/linux64","libhellovr_julia"),:define_julia_module,@__MODULE__)
  CxxWrap.wrapmodule(stage2_path,:define_julia_module,@__MODULE__)
  # const module_functions = CxxWrap.get_module_functions(@__MODULE__)

  # ptr(x :: T) where T = reinterpret(Ptr{T},pointer_from_objref(x))
  # nulltpr(T) = reinterpret(Ptr{T},C_NULL)

  mkRefOrPtr(x::Ptr{Nothing},::Val{Ref{U}}) where U = reinterpret(Ptr{U},x)
  mkRefOrPtr(x::U,::Val{U}) where U = x
  mkRefOrPtr(x,::Val) = Ref(x) # ptr(x)

  # overload for each function that takes a bitstype, allocates it with Ref and passes it as such
  for f ∈ filter(f -> any(T -> T <: Ref && ~(T <: Ptr), f.argument_types), CxxWrap.get_module_functions(@__MODULE__))
    argtypes = f.argument_types
    argsymbols = map((i) -> Symbol(:arg,i[1]), enumerate(argtypes))

    # this makes the smart pointer type union, taken from CxxWrap.build_function_expression
    function map_julia_arg_type_named(fname, t)
      if fname ∈ CxxWrap.__excluded_names
        return t
      end
      return CxxWrap.map_julia_arg_type(t)
    end

    # BUG: I got this looping with nested Ref{Ref{Ref{...}}} when passing a Ptr{Nothing}, where a Ref{Something} was needed
    # fname = mod === nothing ? f.name : (f.name,mod) # TODO: multi-module support?
    fname = CxxWrap.process_fname(f.name) # matches on constructors and other functions (and different modules, but that is to be done)
    args = [map_julia_arg_type_named(f.name, t) for t in f.argument_types]
    newargs = [T <: Ref ? :($arg)           : :($arg :: $T) for (arg,T) in  zip(argsymbols,args)]
    wraps   = [T <: Ref ? :(mkRefOrPtr($arg,Val($T))) : :($arg)       for (arg,T) in  zip(argsymbols,args)]
    expr = :($fname($(newargs...))::$(f.return_type)=$fname($(wraps...)))
    Core.eval(@__MODULE__,expr)
  end

  function getUnion(y::VROverlayIntersectionMaskPrimitive_t)
    off = fieldoffset(VROverlayIntersectionMaskPrimitive_t,2)
    x = Ref(y)
    @GC.preserve x begin
    	x.m_nPrimitiveType == OverlayIntersectionPrimitiveType_Rectangle && return @GC.preserve x unsafe_load(reinterpret(Ptr{IntersectionMaskRectangle_t},pointer_from_objref(x)+off))
    	x.m_nPrimitiveType == OverlayIntersectionPrimitiveType_Circle    && return @GC.preserve x unsafe_load(reinterpret(Ptr{IntersectionMaskCircle_t}   ,pointer_from_objref(x)+off))
    end
  	nothing
  end

  function getUnion(y::VREvent_t)
    x = Ref(y)
    off = fieldoffset(VREvent_t,4)

    # this is reconstructed from the openvr.h comments
    (( t == VREvent_ButtonPress
    || t == VREvent_ButtonUnpress
    || t == VREvent_ButtonTouch
    || t == VREvent_ButtonUntouch
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Controller_t}, pointer_from_objref(x) + off))
    (( t == VREvent_DualAnalog_Press
    || t == VREvent_DualAnalog_Unpress
    || t == VREvent_DualAnalog_Touch
    || t == VREvent_DualAnalog_Untouch
    || t == VREvent_DualAnalog_Move
    || t == VREvent_DualAnalog_ModeSwitch1
    || t == VREvent_DualAnalog_ModeSwitch2
    || t == VREvent_DualAnalog_Cancel
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_DualAnalog_t}, pointer_from_objref(x) + off))
    (( t == VREvent_MouseMove
    || t == VREvent_MouseButtonDown
    || t == VREvent_MouseButtonUp
    || t == VREvent_TouchPadMove
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Mouse_t}, pointer_from_objref(x) + off))
    (( t == VREvent_FocusEnter
    || t == VREvent_FocusLeave
    || t == VREvent_OverlayFocusChanged
    || t == VREvent_DashboardThumbSelected
    || t == VREvent_DashboardRequested
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Overlay_t}, pointer_from_objref(x) + off))
    (( t == VREvent_ScrollDiscrete
    || t == VREvent_ScrollSmooth
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Scroll_t}, pointer_from_objref(x) + off))
    (( t == VREvent_InputFocusCaptured
    || t == VREvent_InputFocusReleased
    || t == VREvent_SceneFocusLost
    || t == VREvent_SceneFocusGained
    || t == VREvent_SceneApplicationChanged
    || t == VREvent_SceneFocusChanged
    || t == VREvent_InputFocusChanged
    || t == VREvent_SceneApplicationSecondaryRenderingStarted
    || t == VREvent_SceneApplicationUsingWrongGraphicsAdapter
    || t == VREvent_ActionBindingReloaded
    || t == VREvent_Quit
    || t == VREvent_ProcessQuit
    || t == VREvent_QuitAborted_UserPrompt
    || t == VREvent_QuitAcknowledged
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Process_t}, pointer_from_objref(x) + off))
    (( t == VREvent_RenderToast
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_Notification_t}, pointer_from_objref(x) + off))
    (( t == VREvent_ShowUI
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_ShowUI_t}, pointer_from_objref(x) + off))
    (( t == VREvent_Input_HapticVibration
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_HapticVibration_t}, pointer_from_objref(x) + off))
    (( t == VREvent_Input_BindingLoadFailed
    || t == VREvent_Input_BindingLoadSuccessful
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_InputBindingLoad_t}, pointer_from_objref(x) + off))
    (( t == VREvent_Input_ActionManifestLoadFailed
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_InputActionManifestLoad_t}, pointer_from_objref(x) + off))
    (( t == VREvent_Input_ProgressUpdate
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_ProgressUpdate_t}, pointer_from_objref(x) + off))
    (( t == VREvent_SpatialAnchors_PoseUpdated
    || t == VREvent_SpatialAnchors_DescriptorUpdated
    || t == VREvent_SpatialAnchors_RequestPoseUpdate
    || t == VREvent_SpatialAnchors_RequestDescriptorUpdate
    )) && return @GC.preserve x unsafe_load(reinterpret(Ptr{VREvent_SpatialAnchor_t}, pointer_from_objref(x) + off))

    # ?? VREvent_Reserved_t reserved;
    # ?? VREvent_Status_t status;
    # ?? VREvent_Keyboard_t keyboard;
    # ?? VREvent_Ipd_t ipd;
    # ?? VREvent_Chaperone_t chaperone;
    # ?? VREvent_PerformanceTest_t performanceTest;
    # ?? VREvent_TouchPadMove_t touchPadMove;
    # ?? VREvent_SeatedZeroPoseReset_t seatedZeroPoseReset;
    # ?? VREvent_Screenshot_t screenshot;
    # ?? VREvent_ScreenshotProgress_t screenshotProgress;
    # ?? VREvent_ApplicationLaunch_t applicationLaunch;
    # ?? VREvent_EditingCameraSurface_t cameraSurface;
    # ?? VREvent_MessageOverlay_t messageOverlay;
    # ?? VREvent_Property_t property;
    # ?? VREvent_WebConsole_t webConsole;
  end

  function GetStringTrackedDeviceProperty(vrsystem,deviceindex,prop)::String
    err = TrackedProp_Success
    len = GetStringTrackedDeviceProperty(vrsystem,deviceindex,prop,Cstring(C_NULL),0,err)
    buf = Array{UInt8,1}(undef,len)
    GetStringTrackedDeviceProperty(vrsystem,deviceindex,prop,Cstring(pointer(buf)),len,err)
    err != TrackedProp_Success && println("Warning: GetStringTrackedDeviceProperty failed with error $(err)")
    return unsafe_string(pointer(buf))
  end

  function __init__()
    @initcxx
  end
end


text:
using NumericalMethodsforEngineers, Gadfly

old = pwd()
ProjDir = Pkg.dir("NumericalMethodsforEngineers", "examples", "NumericalMethodsforEngineers", "Ch07", "IVP")
cd(ProjDir)

f(x::Float64, y::Vector{Float64}) = [y[2], 2.0*y[1]-3.0*y[2]+3.0*x^2]
steps = 4
h = 0.05

x = 0.0
y = [1.0, 0.0]

r = Array{Float64,2}[]
push!(r, euler(f, x, y, steps, h))
push!(r, modified_euler(f, x, y, steps, h))
push!(r, mid_point_euler(f, x, y, steps, h))
push!(r, runga_kutta_4(f, x, y, steps, h))
println(r)

xs = Float64[(i-1)*h for i in 1:steps+1]
pa = Plot[]
titles = ["Euler", "Modified_Euler", "Mid_Point_Euler", "Runga_Kutta_4"]
for i in 1:4
  y1 = r[i][:, 2]
  y2 = r[i][:, 3]
  p = plot(
    layer(x=xs, y=y1, Geom.line,
      color=repeat([symbol("y1")], inner=[3])),
    layer(x=xs, y=y2, Geom.line,
      color=repeat([symbol("y2")], inner=[1])),
    Scale.color_discrete_manual("darkred", "red", "darkblue","darkgreen"),
    Guide.colorkey("Legend"),
    Guide.xlabel("x", orientation=:horizontal),
    Guide.ylabel("w(x)", orientation=:vertical),
    Guide.title(titles[i]))
  push!(pa, p)
end

draw(SVG("Ex7.3.svg", 8inch, 9.5inch), vstack(pa...))
# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@osx ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
  isfile("Ex7.3.svg") &&
    run(`open -a $(JULIA_SVG_BROWSER) "Ex7.3.svg"`)
  end : println()

cd(old)


text:
using Distributions, SearchDiscoveryReplication, NLopt, GalacticOptim
using BSON
using PDMats # required for BSON to load data again

##############################################
## Run simulation for main tables

# Set parameter values
β0 = [1.0 , -1.0,3.5]
cs0 = 0.03
cd0 = 0.06
dChars = MvNormal([2., 3.5], [3.0, 1.0])
outDummy = true
dV = Normal(0, 0)
nCons = 2000
# By having many products in total, it is unlikely that the same product is shown to two
# different consumers
nProdTot = 1000000
nProdCons = 30
nA0 = 1

# Run simulation
d, wSD, outFI, outRS, outWf, outW = runSimulation(β0, cs0, cd0, dChars, outDummy, dV, nCons, nProdTot,
										nProdCons,nA0;
										seed =24,
										options = estimOptions(algo = Opt(:LN_BOBYQA,1)))

# Save results as bson

BSON.@save "gen/results_spec0.bson" wSD outFI outRS outWf outW β0 cs0 cd0 dChars outDummy dV nA0
# BSON.@load "gen/results_spec0.bson" wSD outFI outRS outWf outW β0 cs0 cd0 dChars outDummy dV nA0

##############################################
## Gather results  and print to tables
estimOut, welfareChangesOut, demandOut, dStats = gatherResults("gen/results_spec0.bson")

# TABLE 2: Estimation Results
colnames = 	["\$ \\beta_2 \$", "\$ \\beta_1 / | \\beta_2 |\$",
				"\$ \\beta_3 / | \\beta_2 | \$ ",
				"\$ c_s / | \\beta_2 | \$", "\$ c_d / | \\beta_2 |\$"]
rownames = ["\\emph{SD}","\\emph{DS1}","\\emph{DS2}","\\emph{RS}","\\emph{FI}"]
td = deepcopy(estimOut')[:,1:end-1]

# Normalize with price coefficient
for i in [1,3,4,5]
	global td[:,i] = td[:,i] ./ abs.(td[:,2])
end
td = td[:,[2,1,3,4,5]]

# nSearches and share outside Option
colnames = prepend!(colnames,["\\#Searches","Purchases (\\%)"])
dStats[:,2] .= (1 .- dStats[:,2]) .* 100  # Make into % share purchases

td = hcat(dStats,td)

title =  "Estimated Coefficients and Search Set Size"
notes = "Estimation from a simulated dataset with 2,000 consumers and 30 products per consumer.
Characteristics are independent draws (across consumers and products) from
\$x_{1j} \\sim N(2,3.0)\$, \$x_{2j} \\sim N(3.5,1.0)\$ and \$y_j \\sim N(0,1)\$, , with parameters in the
estimated models denoted by \$c^{RS}=c_s\$, \$c_j^{DS1} = c_s\$ and \$c_J^{DS2} = c_s + c_d h_j\$. The third characteristic is an outside dummy.
The data is generated based on the \\emph{SD} model with \$ n_d=|A_0|=1\$. The first two columns are based either
on the generated data (SD) or estimated by generating 5,000 search paths for each consumer. "

writeLatexTable(td,"./gen/table2_coefs.tex",title,
				colnames , rownames, notes;
				label = "tab:coef_estim")

# TABLE 3: Counterfactuals
colnames = repeat(["\$\\Delta CS \$", "\$\\Delta D_1\$","\$\\Delta D_5 \$"],2)

td = zeros(size(welfareChangesOut[1],1),1)
for i in 1:2
	global td = hcat(td,welfareChangesOut[1][:,1+i],demandOut[1][:,i],demandOut[2][:,i])
end
td = td[:,2:end]

title = "Counterfactuals"
notes = "Results from simulated counterfactuals based on Table \\ref{tab:coef_estim}, where (i) all costs are set to zero and (ii) the price for the 5th
is reduced by 1 \\% for each consumer. All changes are expressed in \\% relative to the baseline. Demand and consumer
surplus are calculated by averaging across 5,000 simulated search paths for each consumer."

writeLatexTable(td,"./gen/table3_counterfactuals.tex",title,
				colnames , rownames, notes;
				label = "tab:counterfactuals" )

# Insert table headers
flines = open(readlines,"gen/table3_counterfactuals.tex")
addedLines = ["&& \\multicolumn{3}{c}{Remove costs} & \\multicolumn{3}{c}{\$\\Delta p_5 = -1\\%\$} \\\\ ",
				"\\cline{3-5} \\cline{6-8}"]

flines = vcat(flines[1:5],addedLines,flines[6:end])
f = open("gen/table3_counterfactuals.tex","w")
for l in flines
	println(f,l)
end
close(f)


text:
dv = @data([NA, 3, 2, 5, 4])

dv

mean(dv)

dropna(dv)

mean(dropna(dv))

convert(Array, dv) #fails...

convert(Array, dropna(dv)) # works fine...

# or overwrite the NA value like:

dv[1] = 3

#then

convert(Array, dv) # also works

# replaciung and converting in one go:

dv = @data([NA, 3, 2, 5, 4])

mean(convert(Array, dv, 11))

# similar for n-dinemsianal arrays...

dm = @data([NA 0.0; 0.0 1.0])
dm * dm

## The Data frame type

df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"]) ## how to build a data frame
df

# can also be done in stages

df = DataFrame()
df[:A] = 1:8
df[:B] = ["M", "F", "F", "F", "M", "M", "F", "M"]

df

# check size

nrows = size(df, 1)
ncols = size(df,2)

head(df)

df[1:3, :B]

# getting a summary

describe(df)

# access columns via this
mean(df[1])
median(df[1])

# or via names

mean(df[:A])
median(df[:A])

# columnwise functions:

df = DataFrame(A = 1:4, B = randn(4))
colwise(cumsum, df) # for all columns

cumsum(df[:A]) # here just the one column returns just one vector

## more complex things

using RDatasets
using Requests
url = "http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"

# read without error checking

df = readtable(Requests.get_streaming(url))

# and with error checking

r = Requests.get_streaming(url; timeout = 30.0)

if r.response.status / 100 != 2
        error("Error Downloading Data")
end
df = readtable(r)

df = DataFrame(A = 1:10)

writetable("output.csv", df)

a = DataFrame(ID = [1, 2], Name = ["A", "B"])
b = DataFrame(ID = [1, 3], Job = ["Doctor", "Lawywer"])

iris = df

by(iris, :Species, size)

#time stuff

now()
Date(2018, 6, 18)
Date("20180618", DateFormat("yyyymmdd"))
Date("2018-06-18", DateFormat("y-m-d"))

t = Date(now())
t2 = Date(now() - Dates.Day(7))

Dates.firstdayofweek(Date(now() - Dates.Day(1)))

ceil(DateTime(now()), Dates.Minute(30))

Dates.dayname(t2)

dr = Dates.Date(2018):Dates.Date(now())
Dates.dayofweek(dr) .== Dates.Tue

sum(1:10) do x
        2x
end

a = 1:10
filter(a) do x
        isodd(x)
end

dr = Dates.Date(now() - Dates.Day(12)):Dates.Date(now());

datums = filter(dr) do x
        ((Dates.dayofweek(x) == Dates.Mon) || (Dates.dayofweek(x) == Dates.Tue)) && (Dates.Week(x) < Dates.Week(now()))
end

Dates.firstdayofweek(Date(now()-Dates.Day(5))) + Dates.Day(1)
Dates.firstdayofweek(Date(now()-Dates.Day(5)))

round(convert(Int,(ceil(DateTime(now()), Dates.Minute(30)) - DateTime(now()))) / 1000)

maxVolume = 50000
timeBlock = 60 * 30

currentVolume = 20000
timeCurrent = 1800 - round(convert(Int,(ceil(DateTime(now()), Dates.Minute(30)) - DateTime(now()))) / 1000)

(currentVolume/timeCurrent) / (maxVolume/timeBlock)


text:
"""
    showgeometry!(scene, vs, ns; arrow = 0.5)

Show pointcloud with normals.
"""
function showgeometry!(scene, vs::Array{SVector{3,F},1}, ns::Array{SVector{3,F},1}; arrow = 0.5, kwargs...) where F
    plns = normalsforplot(vs, ns, arrow)
    scatter!(scene, vs; kwargs...)
    linesegments!(scene, plns, color = :blue)
    #cam3d!(scene)
    scene
end

function showgeometry!(scene, vs, ns; arrow = 0.5, kwargs...)
    vsn = [SVector{3, Float64}(i) for i in vs]
    nsn = [SVector{3, Float64}(i) for i in ns]
    showgeometry!(scene, vsn, nsn; arrow=arrow, kwargs...)
end

function showgeometry!(scene, m; arrow = 0.5, kwargs...)
    vsn = [SVector{3, Float64}(i) for i in m.vertices]
    nsn = [SVector{3, Float64}(i) for i in m.normals]
    showgeometry!(scene, vsn, nsn; arrow=arrow, kwargs...)
end

function showgeometry(m; arrow=0.5, kwargs...)
    s = Scene()
    showgeometry!(s, m; arrow=arrow, kwargs...)
end

function showgeometry(vs, ns; arrow=0.5, kwargs...)
    s = Scene()
    showgeometry!(s, vs, ns; arrow=arrow, kwargs...)
end

function showcandlength(ck)
    for c in ck
        println("candidate length: $(length(c.inpoints))")
    end
end

function showshapes!(s, pointcloud, candidateA; plotleg=true, texts=nothing, kwargs...)
    colscheme = ColorSchemes.gnuplot
    colA = get.(Ref(colscheme), range(0, stop=1, length=(size(candidateA,1)+1)))
    colA = deleteat!(colA, 1)
    texts_ = []
    for i in 1:length(candidateA)
        ind = candidateA[i].inpoints
        push!(texts_, RANSAC.strt(candidateA[i].candidate.shape)*"$i")
        scatter!(s, pointcloud.vertices[ind], color = colA[i]; kwargs...)
    end

    usetexts = texts === nothing ? texts_ : texts
    if plotleg
        if s.attributes[:show_axis][]
            sl = legend(s.plots[2:end], usetexts)
        else
            sl = legend(s.plots, usetexts)
        end
        sn = Scene(clear=false, show_axis = false, resolution=(1920,1080))
        sn.center=false
        return vbox(s, sl, parent=sn)
    else
        return s
    end
end

function showshapes(pointcloud, candidateA; plotleg=true, kwargs...)
    sc = Scene()
    showshapes!(sc, pointcloud, candidateA; plotleg=plotleg, kwargs...)
end

function getrest(pc)
    return findall(pc.isenabled)
end

function showtype(l)
    for t in l
        println(t.candidate.shape)
    end
end

function showbytype!(s, pointcloud, candidateA, plotleg=true; kwargs...)
    colors_ = []
    texts_ = []
    for i in eachindex(candidateA)
        c = candidateA[i]
        ind = c.inpoints
        if c.candidate.shape isa FittedCylinder
            colour = :red
            push!(colors_, :red)
            push!(texts_, "Cylinder$i")
        elseif c.candidate.shape isa FittedSphere
            colour = :green
            push!(colors_, :green)
            push!(texts_, "Sphere$i")
        elseif c.candidate.shape isa FittedPlane
            colour = :orange
            push!(colors_, :orange)
            push!(texts_, "Plane$i")
        elseif c.candidate.shape isa FittedCone
            colour = :blue
            push!(colors_, :blue)
            push!(texts_, "Cone$i")
        elseif c.candidate.shape isa AbstractTranslationalSurface
            colour = :purple
            push!(colors_, :purple)
            push!(texts_, "Translational$i")
        end
        scatter!(s, pointcloud.vertices[ind], color = colour; kwargs...)
    end
    if plotleg
        if s.attributes[:show_axis][]
            sl = legend(s.plots[2:end], texts_)
        else
            sl = legend(s.plots, texts_)
        end
        sn = Scene(clear=false, show_axis = false, resolution=(1920,1080))
        sn.center=false
        return vbox(s, sl, parent=sn)
    else
        return s
    end
end

function showbytype(pointcloud, candidateA, plotleg=true; kwargs...)
    sc = Scene()
    showbytype!(sc, pointcloud, candidateA, plotleg; kwargs...)
end

function plotshape(shape::FittedShape; kwargs...)
    plotshape!(Scene(), shape; kwargs...)
end

function plotshape!(sc, shape::FittedPlane; scale=(1.,1.), color=(:blue, 0.1))
    # see project2plane
    o_z = normalize(shape.normal)
    o_x = normalize(arbitrary_orthogonal2(o_z))
    o_y = normalize(cross(o_z, o_x))

    p1 = shape.point
    p2 = p1 + scale[1]*o_x
    p3 = p1 + scale[1]*o_x + scale[2]*o_y
    p4 = p1 + scale[2]*o_y

    mesh!(sc, [p1,p2,p3], color=color, transparency=true)
    mesh!(sc, [p1,p3,p4], color=color, transparency=true)
    sc
end

function plotshape!(sc, shape::FittedSphere; scale=(1.,), color=(:blue, 0.1))
    mesh!(sc, Sphere(Point(shape.center), scale[1]*shape.radius), color=color, transparency=true)
end

function plotshape!(sc, shape::FittedCylinder; scale=(1.,), color=(:blue, 0.1))
    o = Point(shape.center)
    extr = o+scale[1]*Point(normalize(shape.axis))
    mesh!(sc, Cylinder(o, extr, shape.radius), color=color, transparency=true)
end

function shiftplane!(sc, p::FittedPlane, dist; kwargs...)
    newo = p.point+dist*normalize(p.normal)
    newp = FittedPlane(true, newo, p.normal)
    plotshape!(sc, newp; kwargs...)
end

function plotimplshape!(sc, shape::ImplicitPlane; scale=(1.,1.), color=(:blue, 0.1))
    fp = FittedPlane(true, shape.point, shape.normal)
    plotshape!(sc, fp; scale=scale, color=color)
end

function plotimplshape!(sc, shape::ImplicitSphere; scale=(1.,1.), color=(:blue, 0.1))
    fs = FittedSphere(true, shape.center, shape.radius, true)
    plotshape!(sc, fs; scale=scale, color=color)
end

function plotimplshape!(sc, shape::ImplicitCylinder; scale=(1.,1.), color=(:blue, 0.1))
    fc = FittedCylinder(true, shape.axis, shape.center, shape.radius, true)
    plotshape!(sc, fc; scale=scale, color=color)
end

function plotimplshape!(sc, shape::ImplicitTranslational; scale=(1.,1.), color=(:blue, 0.1))
    @warn "Plotting of ImplicitTranslational is not implemented yet!"
    return sc
end

function plotimplshape(shape::CSGBuilding.AbstractImplicitSurface; scale=(1.,1.), color=(:blue, 0.1))
    s = Scene()
    plotimplshape!(s, shape; scale=scale, color=color)
end

function givelargest(scoredshapes)
    sizes = [size(s.inpoints, 1) for s in scoredshapes]
    mind = argmax(sizes)
    println("Best: $mind. - $(scoredshapes[mind])")
    return scoredshapes[mind]
end

"""
    drawcircles!(sc, points, r; kwargs...)

Use `strokecolor=:blue` or what you want.
"""
function drawcircles!(sc, points, r; kwargs...)
    ps = [Makie.Point2f0(p) for p in points]
    poly!(sc, [Circle(p, Float32(r)) for p in ps], color = (:blue, 0.0), transparency = true, strokewidth=1; kwargs...)
    sc
end

function wframe(ps, trs; kwargs...)
	vrts = [Point3(i[1], i[2], 0.) for i in ps]
	fces = [Face{3,Int}(i) for i in trs]
	tm = HomogenousMesh(vrts, fces, [], [], [], [], [])
	Makie.wireframe(tm; kwargs...)
end

function plotspantree!(s, points, tree; kwargs...)
	for (i, obj) in enumerate(tree)
		p2 = points[obj]
		lines!(s, (x->x[1]).(p2), (x->x[2]).(p2); kwargs...)
	end
	s
end


text:
# All functionality corresponding to Predictionerror as Loss function Γ.

## Constructors:

abstract type AbstractMCDTSpredictionType end

abstract type AbstractPredictionLoss{P} end

abstract type AbstractPredictionMethod end

abstract type AbstractLocalPredictionMethod{T} <: AbstractPredictionMethod end

"""
    Prediction_error <: AbstractLoss

    Constructor for the Prediction_error loss function.

    ## Fieldnames
    * `PredictionType::MCDTSpredictionType`: Determines the prediction type by
      setting a prediction-method and the way the prediction error is measured,
      see [`MCDTSpredictionType`](@ref).
    * `threshold::Float`: A threshold for the sufficient minimum prediction error
      for the current embedding. When the prediction error, specified in
      `PredictionType`, falls below this threshold in an embedding cycle the
      embedding stops.
    * `samplesize::Real = 1.`: the fraction of all phase space points
      to be considered in the computation of the prediction error under the given
      `PredictionType`.
    * `error_wheights::Vector{Real} = [0 1]`: The wheights for determining the prediction
      error. The first element of this Vector is the wheight for the insample prediction
      error and the second elements corresponds to the wheight for the out-of-sample
      prediction error. By default only the out-of-sample error will be used (i.e. wheights [0 1]).
      For specifying the prediction horizon see [`local_model`](@ref)

    ## Defaults
    * When calling `Prediction_error()`, a Prediction_error-object is created,
      which uses the threshold 0, i.e. no threshold, a zeroth-order predictor
      (see [`MCDTSpredictionType`](@ref), [`PredictionLoss`](@ref) &
      [`local_model`](@ref)) and a full phase space sample (samplesize=1).
"""
struct Prediction_error <: AbstractLoss
    PredictionType::AbstractMCDTSpredictionType
    threshold::AbstractFloat
    samplesize::Real
    error_wheights::Vector{Real}
    # Constraints and Defaults
    Prediction_error(x ,y=0, z=1, zz=[0;1.]) = begin
        @assert y ≥ 0 "A positive threshold for the prediciton loss must be chosen."
        @assert 0 < z ≤ 1. "The samplesize must be in the interval (0 1]."
        @assert length(zz) == 2 && 0 <= zz[1] ≤ 1 && 0 <= zz[2] ≤ 1 "The errorweights
        must be passed in a vector of length 2. The elements of this vector need to be
        in the interval [0 1]."
        new(x, y, z, zz)
    end

end
Prediction_error() = Prediction_error(MCDTSpredictionType())

"""
    MCDTSpredictionType <: AbstractMCDTSpredictionType

    Constructor, which determines the way how predictions are made technically.

    ## Fieldnames
    * `loss::AbstractPredictionLoss`: Indicates the way of computing the prediction error.
       See [`PredictionLoss`](@ref) for information on how to construct this object.
    * `method::AbstractPredictionMethod`: The method based on the state space reconstruction,
       which makes the actual prediction. See [`local_model`](@ref)

    ## Default settings
    * When calling `MCDTSpredictionType()` a MCDTSpredictionType-object is constructed
      with a `local_zeroth`-predictor [`local_model`](@ref), using 2 nearest neighbors
      and a 1-step-ahead-prediction. The loss-function is the root mean squared prediction
      error over all components [`PredictionLoss`](@ref).
"""
struct MCDTSpredictionType <: AbstractMCDTSpredictionType
    loss::AbstractPredictionLoss
    method::AbstractPredictionMethod
    # Constraints and Defaults
    MCDTSpredictionType(x,y) = new(x,y)
    MCDTSpredictionType(x) = new(x,local_model())
    MCDTSpredictionType() = new(PredictionLoss(), local_model())
end

"""
    PredictionLoss <: AbstractPredictionLoss

    Constructor, which indicates the way of computing the prediction error. This
    object is used for the constructor, which determines the way how predictions are
    made methodologically [`MCDTSpredictionType`](@ref).

    ## Fieldnames
    * `type::Int` is an integer, which encodes the type of prediction error:

    - For `type = 1` the root mean squared prediction error over the first component,
      i.e. the timeseries, which needs to be predicted, is used. (default)
    - For `type = 2` the root mean squared prediction error over all components
      (dimensionality of the state space) is used.
    - For `type = 3` the mean Kullback-Leibler Distance of the predicted and the true
      values of the first component, i.e. the timeseries, which needs to be predicted, is used.
    - For `type = 4` the mean Kullback-Leibler Distance of the predicted and the true
      values over all components (dimensionality of the state space) is used.

    ## Defaults
    * When calling `PredictionLoss()` a PredictionLoss-object is constructed with
      fieldname `type = 1` (≡root mean squared prediction error over all components)
"""
struct PredictionLoss{t} <: AbstractPredictionLoss{t}
    type::Int
    # Constraints and Defaults
    PredictionLoss(x=1) = begin
        @assert x == 1 || x == 2 || x == 3 || x == 4
        new{x}(x)
    end
end

"""
    local_model <: AbstractLocalPredictionMethod

    Constructor, which indicates the local state space prediction model.

    ## Fieldnames
    * `method::String`: Could be `"zeroth"` (averaged `Tw`-step-ahead image of the
     `KNN`-nearest neighbors) or `"linear"` (local linear regression on the
     `KNN`-nearest neighbors).
    * `KNN::Int`: The number of considered nearest neighbors.
    * `Tw_out::Int` : The prediction horizon in sampling units for the out-of-sample prediction.
    * `Tw_in::Int` : The prediction horizon in sampling units for the in-sample prediction.
    * `trials::Int` : The number of different out-of-sample prediction trials.

    ## Defaults
    * When calling `local_model()` a local_model-object is constructed with a zeroth
      order prediction scheme, 2 nearest neighbors, a 1-step-ahead in-sample prediction,
      a 10-step-ahead out-of-sample prediction and 20 out-of-sample trials.
    * When calling `local_model(method)` a local_model-object is constructed with a
      `method`-prediction scheme, 2 nearest neighbors, a 1-step-ahead in-sample prediction,
      a 10-step-ahead out-of-sample prediction and 20 out-of-sample trials.
    * When calling `local_model(method,KNN)` a local_model-object is constructed with a
     `method`-prediction scheme, `KNN` nearest neighbors, a 1-step-ahead in-sample prediction,
     a 10-step-ahead out-of-sample prediction and 20 out-of-sample trials.
"""
struct local_model{m} <: AbstractLocalPredictionMethod{m}
    method::String
    KNN::Int
    Tw_out::Int
    Tw_in::Int
    trials::Int

    # Constraints and Defaults
    local_model(x="zeroth", y=2, z=10, zz=1, zzz=20) = begin
        @assert x in ["zeroth", "linear"]
        @assert y > 0
        @assert z > 0
        @assert zz > 0
        @assert zzz > 0
        m = Symbol(x)
        new{m}(x,y,z,zz,zzz)
    end
end

## Functions:

"""
    Return the loss based on a `Tw`-step-ahead local-prediction.
"""
function compute_loss(Γ::Prediction_error, Λ::AbstractDelayPreselection, dps::Vector{P}, Y_act::Dataset{D, T}, Ys, τs, w::Int, ts::Int, τ_vals, ts_vals; metric=Euclidean(), kwargs...) where {P, D, T}

    PredictionLoss = Γ.PredictionType.loss
    PredictionMethod = Γ.PredictionType.method
    samplesize = Γ.samplesize
    error_wheights = Γ.error_wheights

    max_idx = get_max_idx(Λ, dps, τ_vals, ts_vals, ts) # get the candidate delays
    isempty(max_idx) && return Float64[], Int64[], []

    costs_insample = zeros(Float64, length(max_idx))
    costs_out_of_sample = zeros(Float64, length(max_idx))
    for (i,τ_idx) in enumerate(max_idx)
        # create candidate phase space vector for this peak/τ-value
        tau_trials = (τ_vals...,τs[τ_idx-1],)
        ts_trials = (ts_vals...,ts,)
        Y_trial = genembed(Ys, tau_trials.*(-1), ts_trials)
        Y_trial1 = deepcopy(Y_trial)

        # make an in-sample prediction for Y_trial (if needed)
        if error_wheights[1]>0
            prediction_insample, ns, temp = insample_prediction(PredictionMethod, Y_trial; samplesize, w, metric, i_cycle=length(τ_vals), kwargs...)
            # compute loss/costs
            costs_insample[i] = compute_costs_from_prediction(PredictionLoss, prediction_insample, Y_trial, PredictionMethod.Tw_in, ns)
        end
        # make an out-of-sample prediction for Y_trial (if needed)
        if error_wheights[2]>0
            costs_out_of_sample[i] = out_of_sample_prediction(PredictionMethod, PredictionLoss, Y_trial; w, metric, i_cycle=length(τ_vals), kwargs...)
        end
    end
    costs = error_wheights[1]*costs_insample .+ error_wheights[2]*costs_out_of_sample

    return costs, max_idx, [nothing for i in max_idx]
end

"""
    insample_prediction(pred_meth::AbstractPredictionMethod, Y::AbstractDataset{D, ET};
            samplesize::Real = 1, K::Int = 3, w::Int = 1, Tw::Int = 1, metric = Euclidean()) → prediction

    Compute an in-sample `Tw`-time-steps-ahead prediction of the data `Y`, using
    the prediction method `pred_meth`. `w` is the Theiler window and `K` the nearest
    neighbors used.

    * `Y`: Dataset (Nt x N_embedd)
    * `K`: Nearest Neighbours
    * `w`: Theiler window
    * `Tw`: Prediction horizon
    * `metric`: Metric for NN search
    * `samplesize`: fraction of considered points in the trajectory
    * `i_cycle`: Which embedding cycling we are predicting for (for future use)

    Note: In case of a local linear prediction method `pred_meth` the number of
    nearest neighbours used gets adapted to 2(D+1) - with D the embedding dimension,
    if the provided `K` is lower than that number.")
"""
function insample_prediction(pred_meth::AbstractLocalPredictionMethod, Y::AbstractDataset{D, ET}; samplesize::Real= 1, w::Int = 1, metric = Euclidean(), i_cycle::Int=1, kwargs...) where {D, ET}

    Tw = pred_meth.Tw_in # total time horizon
    NN = length(Y)-Tw
    if samplesize==1
        ns = 1:NN
        Nfp = length(ns)
    else
        Nfp = Int(floor(samplesize*NN)) # number of considered fiducial points
        ns = sample(vec(1:NN), Nfp, replace = false)  # indices of fiducial points
    end
    prediction_new = deepcopy(Y[ns,:]) # intitial trajectory up to the prediction time horizon
    prediction_old = deepcopy(Y) # intitial trajectory prediction is based on
    for i = 1:Tw
        insample_prediction!(pred_meth, prediction_old, prediction_new, ns; w, K = pred_meth.KNN, Tw_step = i)
    end
    return prediction_new, ns, nothing
end

function insample_prediction!(pred_meth::AbstractLocalPredictionMethod{:zeroth}, prediction_old::AbstractDataset{D, ET},
                prediction_new::AbstractDataset{D, ET}, ns::Union{AbstractRange, AbstractVector};
                w::Int = 1, metric = Euclidean(), K::Int=1, Tw_step::Int=1) where {D, ET}

    NN = length(prediction_new)
    vtree = KDTree(prediction_old[1:NN+Tw_step-1], metric)
    ns_act = ns .+ (Tw_step -1)
    allNNidxs, _ = DelayEmbeddings.all_neighbors(vtree, prediction_old[ns_act], ns_act, K, w)
    ϵ_ball = zeros(ET, K, D) # preallocation
    # loop over each fiducial point
    for (i,v) in enumerate(ns_act)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v

        # determine neighborhood one time step ahead
        @inbounds for (k, j) in enumerate(NNidxs)
            ϵ_ball[k, :] .= prediction_old[j+1] # consider 1-step ahead prediction
        end
        # take the average as a prediction
        prediction_new[i] = mean(ϵ_ball; dims=1)
    end
    for (i,v) in enumerate(ns_act)
        prediction_old[v+1] = prediction_new[i] # update trajectory with the predicted 1-step ahead values
    end
end
function insample_prediction!(pred_meth::AbstractLocalPredictionMethod{:linear}, prediction_old::AbstractDataset{D, ET},
                prediction_new::AbstractDataset{D, ET}, ns::Union{AbstractRange, AbstractVector};
                w::Int = 1, metric = Euclidean(), K::Int=1, Tw_step::Int=1) where {D, ET}

    if K < 2*(D+1)
        K = 2*(D+1)
    end
    NN = length(prediction_new)
    vtree = KDTree(prediction_old[1:NN+Tw_step-1], metric)
    ns_act = ns .+ (Tw_step -1)
    allNNidxs, _ = DelayEmbeddings.all_neighbors(vtree, prediction_old[ns_act], ns_act, K, w)
    prediction = zeros(ET, D) # preallocation
    ϵ_ball = zeros(ET, K, D) # preallocation
    b  = zeros(D) # preallocation
    ar_coeffs = zeros(D, D) # preallocation

    # loop over each fiducial point
    for (i,v) in enumerate(ns_act)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v
        A = ones(K,D) # datamatrix for later linear equation to solve for AR-process
        # determine neighborhood one time step ahead
        @inbounds for (k, j) in enumerate(NNidxs)
            ϵ_ball[k, :] .= prediction_old[j + 1] # consider 1-step ahead prediction
            A[k,:] = prediction_old[j]
        end

        namess = ["X"*string(z) for z = 1:D]
        ee = Meta.parse.(namess)
        formula_expression = Term(:Y) ~ sum(term.(ee))

        for j = 1:D
            data = DataFrame()
            for (cnt,var) in enumerate(namess)
                data[!, var] = A[:,cnt]
            end
            data.Y = ϵ_ball[:,j]

            ols = lm(formula_expression, data)
            b[j] = coef(ols)[1]
            for k = 1:D
                ar_coeffs[j,k] = coef(ols)[k+1]
            end
            prediction[j] = prediction_old[v,:]'*ar_coeffs[j,:] + b[j]
        end
        prediction_new[i] = prediction
    end
    for (i,v) in enumerate(ns_act)
        prediction_old[v+1] = prediction_new[i] # update trajectory with the predicted values
    end
end

"""
    out_of_sample_prediction(pred_meth::AbstractLocalPredictionMethod, pred_loss::AbstractPredictionLoss, Y::AbstractDataset{D, ET};
            K::Int = 3, w::Int = 1, Tw::Int = 1, metric = Euclidean()) → average_cost

    Compute an out-of-sample `Tw`-time-steps-ahead prediction of the data `Y`, using
    the prediction method `pred_meth`. `w` is the Theiler window and `K` the nearest
    neighbors used.

    * `Y`: Dataset (Nt x N_embedd)
    * `K`: Nearest Neighbours
    * `w`: Theiler window
    * `Tw`: Prediction horizon
    * `metric`: Metric for NN search
    * `i_cycle`: Which embedding cycling we are predicting for (for future use)

    Note: In case of a local linear prediction method `pred_meth` the number of
    nearest neighbours used gets adapted to 2(D+1) - with D the embedding dimension,
    if the provided `K` is lower than that number.")
"""
function out_of_sample_prediction(pred_meth::AbstractLocalPredictionMethod, pred_loss::AbstractPredictionLoss,
        Y::AbstractDataset{D, ET}; w::Int = 1, metric = Euclidean(), i_cycle::Int=1, kwargs...) where {D, ET}

    Tw = pred_meth.Tw_out
    num_of_trials = pred_meth.trials
    N = length(Y)
    NN = N-Tw

    # split data into training (90%) and test set (10%)
    N_train = Int(ceil(NN*0.9))
    N_test = NN-N_train

    if num_of_trials ≥ N_test
        println("The test set has not sufficient size to satisfy the given number of out-of-sample trials. This has been adjusted now.")
        num_of_trials = N_test
    end

    ns = sample(vec(N_train+1:N_train+N_test), num_of_trials, replace = false)  # starting indices of trial
    costs = zeros(num_of_trials)
    Threads.@threads for i = 1:num_of_trials
        costs[i] = cost_from_out_of_sample_prediction(pred_meth, pred_loss, Y, ns[i]; w, Tw, metric, kwargs...)
    end
    # return the average of the costs for each trial
    return mean(costs)
end

"""
    Compute costs for the out-of-sample prediction
"""
function cost_from_out_of_sample_prediction(pred_meth::AbstractLocalPredictionMethod{:zeroth}, pred_loss::AbstractPredictionLoss,
        Y::AbstractDataset{D, ET}, ns::Int; w::Int = 1, Tw::Int = 1, metric = Euclidean(), verbose::Bool=false, kwargs...) where {D, ET}

    predicted = iterated_local_zeroth_prediction(Y[1:ns], pred_meth.KNN, Tw; metric, theiler = w, verbose)
    # compute loss/costs
    return compute_costs_from_prediction(pred_loss, predicted, Y[ns+1:ns+Tw], 0, 1:length(ns+1:ns+Tw))
end
function cost_from_out_of_sample_prediction(pred_meth::AbstractLocalPredictionMethod{:linear}, pred_loss::AbstractPredictionLoss,
        Y::AbstractDataset{D, ET}, ns::Int; w::Int = 1, Tw::Int = 1, metric = Euclidean(), verbose::Bool=false, kwargs...) where {D, ET}

    predicted = iterated_local_linear_prediction(Y[1:ns], pred_meth.KNN, Tw; metric, theiler = w, verbose)
    # compute loss/costs
    return compute_costs_from_prediction(pred_loss, predicted, Y[ns+1:ns+Tw], 0, 1:length(ns+1:ns+Tw))
end

"""
    Compute the in-sample prediction costs based on the loss-metric determined
    by PredictionLoss
"""
function compute_costs_from_prediction(PredictionLoss::AbstractPredictionLoss{1}, prediction::AbstractDataset{D, T},
                            Y::AbstractDataset{D, T}, Tw::Int, ns::Union{AbstractRange,AbstractVector}) where {D, T}

    NN = length(ns)
    @assert length(prediction) == length(ns)
    costs = zeros(T, NN, D)
    @inbounds for (i,v) in enumerate(ns)
        costs[i,:] = (Vector(prediction[i]) .- Vector(Y[v+Tw])).^2
    end
    c = sqrt.(mean(costs; dims=1))
    return c[1]
end
function compute_costs_from_prediction(PredictionLoss::AbstractPredictionLoss{2}, prediction::AbstractDataset{D, T},
                            Y::AbstractDataset{D, T}, Tw::Int, ns::Union{AbstractRange,AbstractVector}) where {D, T}

    NN = length(ns)
    @assert length(prediction) == length(ns)
    costs = zeros(T, NN, D)
    @inbounds for (i,v) in enumerate(ns)
        costs[i,:] = (Vector(prediction[i]) .- Vector(Y[v+Tw])).^2
    end
    c = sqrt.(mean(costs; dims=1))
    return mean(c)
end
function compute_costs_from_prediction(PredictionLoss::AbstractPredictionLoss{3}, prediction::AbstractDataset{D, T},
                            Y::AbstractDataset{D, T}, Tw::Int, ns::Union{AbstractRange,AbstractVector}) where {D, T}

    NN = length(ns)
    @assert length(prediction) == length(ns)
    costs = zeros(D)
    for j = 1:D
        costs[j] = compute_KL_divergence(Vector(prediction[:,j]),Y[ns .+ Tw,j])
    end
    return costs[1]
end
function compute_costs_from_prediction(PredictionLoss::AbstractPredictionLoss{4}, prediction::AbstractDataset{D, T},
                            Y::AbstractDataset{D, T}, Tw::Int, ns::Union{AbstractRange,AbstractVector}) where {D, T}

    NN = length(ns)
    @assert length(prediction) == length(ns)
    costs = zeros(D)
    for j = 1:D
        costs[j] = compute_KL_divergence(Vector(prediction[:,j]),Y[ns .+ Tw,j])
    end
    return mean(costs)
end

"""
    local_zeroth_prediction(Y::Dataset, K::Int = 5; kwargs...) → x_pred, e_expect

    Perform a "zeroth" order prediction for the time horizon `Tw` (default = 1). Based
    on `K` nearest neighbours of the last point of the given trajectory `Y`, the
    `Tw`-step ahead prediction is simply the mean of the images of these `K`-nearest
    neighbours. The output `x_pred` is, thus, the `Tw`-step ahead prediction vector.
    The function also returns `e_expect`, the expected error on the prediction `x_pred`,
    computed as the mean of the RMS-errors of all `K`-neighbours-errors.

    Keywords:
    * `metric = Euclidean()`: Metric used for distance computation
    * `theiler::Int = 1`: Theiler window for excluding serially correlated points from
       the nearest neighbour search.
    * `Tw::Int = 1`: The prediction time in sampling units. If `Tw > 1`, a multi-step
      prediction is performed.

"""
function local_zeroth_prediction(Y::Dataset{D,T}, K::Int = 5;
    metric = Euclidean(), theiler::Int = 1, Tw::Int = 1) where {D,T}

    NN = length(Y)
    ns = 1:NN
    vs = Y[ns]
    vtree = KDTree(Y[1:length(Y)-Tw,:], metric)
    allNNidxs, _ = DelayEmbeddings.all_neighbors(vtree, vs, ns, K, theiler)

    ϵ_ball = zeros(T, K, D) # preallocation
    # consider NNs of the very last point of the trajectory
    NNidxs = allNNidxs[end] # indices of k nearest neighbors to v
    # determine neighborhood `Tw` time steps ahead
    @inbounds for (i, j) in enumerate(NNidxs)
        ϵ_ball[i, :] .= Y[j + Tw]
    end
    # take the average as a prediction
    prediction = mean(ϵ_ball; dims=1)
    # predicted prediction error
    error_predict = sum((ϵ_ball .- prediction).^2; dims=1)

    return vec(prediction), vec(error_predict)
end

"""
    iterated_local_zeroth_prediction(Y::Dataset, K::Int = 5, Tw::Int = 2; kwargs...) → Y_predict
    Perform an iterated one step forecast over `Tw` time steps using the local zeroth
    prediction algorithm. `Y_predict` is a Vector of dimension like `Y`.

    Keywords:
    * `metric = Euclidean()`: Metric used for distance computation
    * `theiler::Int = 1`: Theiler window for excluding serially correlated points from
       the nearest neighbour search.
    * `verbose::Bool = false`: When set to `true`, the function prints the actual time
      step, which it is computing.
"""
function iterated_local_zeroth_prediction(Y::Dataset{D,T}, K::Int = 5, Tw::Int = 2;
    metric = Euclidean(), theiler::Int = 1, verbose::Bool = false) where {D,T}

    @assert Tw > 0 "Time horizon must be a positive integer"
    N = length(Y)
    predicted_trajectory = deepcopy(Y)
    for Th = 1:Tw
        if verbose
            println("Compute prediction for time step $Th")
        end
        # iterated one step
        predicted, _ = local_zeroth_prediction(predicted_trajectory, K;
                                            theiler, metric)
        Base.push!(predicted_trajectory, predicted)
    end
    return predicted_trajectory[N+1:end]
end

"""
    local_linear_prediction(Y::Dataset, K::Int = 5; kwargs...) → x_pred, e_expect

    Perform a prediction for the time horizon `Tw` (default = 1) by a locally linear
    fit. Based on `K` nearest neighbours of the last point of the given trajectory
    `Y`, we fit a linear model to these points and their `Tw`-step ahead images. The
    output `x_pred` is, thus, the `Tw`-step ahead prediction vector.
    The function also returns `e_expect`, the expected error on the prediction `x_pred`,
    computed as the mean of the RMS-errors of all `K`-neighbours-errors.

    Keywords:
    * `metric = Euclidean()`: Metric used for distance computation
    * `theiler::Int = 1`: Theiler window for excluding serially correlated points from
       the nearest neighbour search.
    * `Tw::Int = 1`: The prediction time in sampling units. If `Tw > 1`, a multi-step
      prediction is performed.
"""
function local_linear_prediction(Y::Dataset{D,T}, K::Int = 5;
    metric = Euclidean(), theiler::Int = 1, Tw::Int = 1) where {D,T}

    if K < 2*(D+1)
        K = 2*(D+1)
    end

    NN = length(Y)
    ns = 1:NN
    vs = Y[ns]
    vtree = KDTree(Y[1:length(Y)-Tw,:], metric)
    allNNidxs, _ = DelayEmbeddings.all_neighbors(vtree, vs, ns, K, theiler)

    ϵ_ball = zeros(T, K, D) # preallocation
    A = ones(K,D) # datamatrix for later linear equation to solve for AR-process
    # consider NNs of the very last point of the trajectory
    NNidxs = allNNidxs[end] # indices of k nearest neighbors to v

    # determine neighborhood `Tw` time steps ahead
    @inbounds for (i, j) in enumerate(NNidxs)
        ϵ_ball[i, :] .= Y[j + Tw]
        A[i,:] = Y[j]
    end

    # make local linear model of the last point of the trajectory
    prediction = zeros(D)
    b  = zeros(D)
    ar_coeffs = zeros(D, D)
    namess = ["X"*string(i) for i = 1:D]
    ee = Meta.parse.(namess)
    formula_expression = Term(:Y) ~ sum(term.(ee))

    for i = 1:D
        data = DataFrame()
        for (cnt,var) in enumerate(namess)
            data[!, var] = A[:,cnt]
        end
        data.Y = ϵ_ball[:,i]

        ols = lm(formula_expression, data)
        b[i] = coef(ols)[1]
        for j = 1:D
            ar_coeffs[i,j] = coef(ols)[j+1]
        end
        prediction[i] = Y[NN,:]'*ar_coeffs[i,:] + b[i]
    end

    # predicted prediction error
    error_predict = sum((ϵ_ball .- prediction').^2; dims=1)

    return vec(prediction), vec(error_predict)
end

"""
    iterated_local_linear_prediction(Y::Dataset, K::Int = 5, Tw::Int = 2; kwargs...) → Y_predict

    Perform an iterated one step forecast over `Tw` time steps using the local linear
    prediction algorithm. `Y_predict` is a Vector of dimension like `Y`.

    Keywords:
    * `metric = Euclidean()`: Metric used for distance computation
    * `theiler::Int = 1`: Theiler window for excluding serially correlated points from
       the nearest neighbour search.
    * `verbose::Bool = false`: When set to `true`, the function prints the actual time
      step, which it is computing.
"""
function iterated_local_linear_prediction(Y::Dataset{D,T}, K::Int = 5, Tw::Int = 2;
    metric = Euclidean(), theiler::Int = 1, verbose::Bool = false) where {D,T}

    @assert Tw > 0 "Time horizon must be a positive integer"
    N = length(Y)
    predicted_trajectory = deepcopy(Y)
    for Th = 1:Tw
        if verbose
            println("Compute prediction for time step $Th")
        end
        # iterated one step
        predicted, _ = local_linear_prediction(predicted_trajectory, K;
                                            theiler, metric)
        Base.push!(predicted_trajectory, predicted)
    end
    return predicted_trajectory[N+1:end]
end

max_stars_repo_path: src/Cost functions/Prediction_error.jl

max_issues_repo_path: src/Cost functions/Prediction_error.jl

max_forks_repo_path: src/Cost functions/Prediction_error.jl


text:
using Test, SerializationCaches, OrderedCollections, Serialization

@testset "`SerializationCache`" begin
    mktempdir() do tmp
        in_memory_limit = 5
        file_limit, file_gc_ratio = 10, 0.1
        total_limit = in_memory_limit + file_limit
        offset = 3
        cache = SerializationCache(tmp; in_memory_limit=in_memory_limit,
                                   file_limit=file_limit, file_gc_ratio=file_gc_ratio)
        range = 1:(total_limit + offset)
        for i in range
            @test i == fetch!(() -> i, cache, string(i))
        end
        file_keys = string.((offset + 1):(offset + file_limit))
        @test cache.file_keys == OrderedSet(file_keys)
        in_memory_keys = string.((offset + file_limit + 1):(offset + file_limit + in_memory_limit))
        @test all(keys(cache.in_memory) .== in_memory_keys)
        let hits = 0
            for r in (range, reverse(range))
                hits += length(r)
                for i in r
                    @test i == fetch!(cache, string(i)) do
                        hits -= 1
                        return i
                    end
                end
            end
            @test hits == 10
            @test cache.file_keys == OrderedSet(["15", "14", "13", "12", "11", "10", "9", "8", "7", "6"])
            @test all(keys(cache.in_memory) .== ["5", "4", "3", "2", "1"])
        end
        put_range = (total_limit + offset):(total_limit + offset + total_limit)
        for i in put_range
            put!(cache, string(i), ()-> i)
        end
        @test all(keys(cache.in_memory) .== ["29", "30", "31", "32", "33"])
        @test cache.file_keys == OrderedSet(["19", "20", "21", "22", "23", "24", "25", "26", "27", "28"])

        @test isdir(set_up_cache_path(joinpath(tmp, "big/made/up/path")))
    end
end

struct NoFile
    i::Int
end

# Test that when `in_memory_limit` is 0, deserializing an existing file does not
# (redundantly) call `put`
function Base.put!(cache::SerializationCache, key::AbstractString, item::NoFile;
                   directly_to_file::Bool=false)
   throw(error("Should not hit this!"))
end

@testset "No in_memory_limit" begin
    mktempdir() do tmp
        in_memory_limit = 0
        file_limit, file_gc_ratio = 10, 0.1

        # Set up existing cache
        range = 1:file_limit
        for i in range
            serialize(joinpath(tmp, "$i.jls"), NoFile(i))
        end

        # Set up cache
        cache = SerializationCache(tmp; in_memory_limit=0,
                                   file_limit=file_limit, file_gc_ratio=file_gc_ratio)

        # Now test that `fetch!()` never calls `put`, since it should be able to
        # grab all files
        for i in range
            @test NoFile(i) == fetch!(() -> (error("shouldn't need to call this...")), cache, string(i))
        end
    end
end

@testset "No `file_limit`" begin
    mktempdir() do tmp
        in_memory_limit = 0
        cache_file_limit_5 = SerializationCache(tmp; in_memory_limit=in_memory_limit,
                                                file_limit=5, file_gc_ratio=0.1)
        range = 1:10
        for i in range
            @test i == fetch!(() -> i, cache_file_limit_5, string(i))
        end
        @test length(readdir(cache_file_limit_5.path)) == 5
    end

    mktempdir() do tmp
        in_memory_limit = 0
        cache_no_file_limit = SerializationCache(tmp; in_memory_limit=in_memory_limit,
                                                 file_limit=nothing, file_gc_ratio=0.1)
        range = 1:10
        for i in range
            @test i == fetch!(() -> i, cache_no_file_limit, string(i))
        end
        @test length(readdir(cache_no_file_limit.path)) == length(range)
    end
end


text:
@testset "CHELSEA BioClim" begin
    using RasterDataSources: rasterurl

    @test rastername(CHELSA{BioClim}, 5) == "CHELSA_bio10_05.tif"

    bioclim_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "CHELSA", "BioClim")
    @test rasterpath(CHELSA{BioClim}) == bioclim_path
    @test rasterpath(CHELSA{BioClim}, 5) == joinpath(bioclim_path, "CHELSA_bio10_05.tif")

    @test rasterurl(CHELSA) ==
        URI(scheme="https", host="os.zhdk.cloud.switch.ch", path="/envicloud/chelsa/chelsa_V1/")
    @test rasterurl(CHELSA{BioClim}, 5) ==
        URI(scheme="https", host="os.zhdk.cloud.switch.ch", path="/envicloud/chelsa/chelsa_V1/climatologies/bio/CHELSA_bio10_05.tif")

    raster_path = joinpath(bioclim_path, "CHELSA_bio10_05.tif")
    @test getraster(CHELSA{BioClim}, 5) == raster_path
    @test getraster(CHELSA{BioClim}, (5,)) == (raster_path,)
    @test isfile(raster_path)
end


text:
##############################
### Get Edges of the Graph ###
##############################

function get_edges(nodes::Dict{Int,T},roadways::Vector{OpenStreetMapX.Way}) where T<:Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF}
    oneway_roads = map(OpenStreetMapX.oneway,roadways)
    reverse_roads = map(OpenStreetMapX.reverseway,roadways)
	classes = OpenStreetMapX.classify_roadways(roadways)
    edges = Dict{Tuple{Int,Int},Int}()
    for i = 1:length(roadways)
        for j = 2:length(roadways[i].nodes)
            n0 = roadways[i].nodes[j-1]
            n1 = roadways[i].nodes[j]
            start = n0 * !reverse_roads[i] + n1 * reverse_roads[i]
            fin = n0 * reverse_roads[i] + n1 * !reverse_roads[i]
			edges[(start,fin)] = classes[roadways[i].id]
            oneway_roads[i] || (edges[(fin,start)] = classes[roadways[i].id])
        end
    end
	return collect(keys(edges)), collect(values(edges))
end

#################################
### Get Vertices of the Graph ###
#################################

function get_vertices(edges::Array{Tuple{Int64,Int64},1})
    graph_nodes = unique(reinterpret(Int, edges))
    vertices = Dict{Int,Int}(zip(graph_nodes, 1:length(graph_nodes)))
end

###################################
### Get Distances Between Edges ###
###################################

function distance(nodes::Dict{Int,T},edges::Array{Tuple{Int64,Int64},1}) where T<:Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF}
    distances = Float64[]
    for edge in edges
        dist = OpenStreetMapX.distance(nodes[edge[2]],nodes[edge[1]])
        push!(distances,dist)
    end
    return distances
end

####################################################
###	For Each Feature Find the Nearest Graph Node ###
####################################################

function features_to_graph(m::OpenStreetMapX.MapData, features::Dict{Int,Tuple{String,String}}) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})
    features_to_nodes = Dict{Int,Int}()
    sizehint!(features_to_nodes,length(features))
    for (key,value) in features
        if !haskey(m.v,key)
            features_to_nodes[key] = OpenStreetMapX.nearest_node(m,features[key])
        else
            features_to_nodes[key] = key
        end
    end
    return features_to_nodes
end

#########################################
### Find Routes - Auxiliary Functions ###
#########################################

### Dijkstra's Algorithm ###
function dijkstra(m::OpenStreetMapX.MapData, w::SparseArrays.SparseMatrixCSC{Float64,Int64}, start_vertex::Int)
    return LightGraphs.dijkstra_shortest_paths(m.g, start_vertex, w)
end

### Transpose distances to times ###

function network_travel_times(m::OpenStreetMapX.MapData, class_speeds::Dict{Int,Float64} = OpenStreetMapX.SPEED_ROADS_URBAN)
    @assert length(m.e) == length(m.w.nzval)
    indices = [(m.v[i],m.v[j]) for (i,j) in m.e]
    w = Array{Float64}(undef,length(m.e))
    for i = 1:length(w)
        w[i] = 3.6 * (m.w[indices[i]]/class_speeds[m.class[i]])
    end
    return w
end

### Create a Sparse Matrix for a given vector of weights ###

function create_weights_matrix(m::OpenStreetMapX.MapData,weights::Vector{Float64})
    w = SparseArrays.spzeros(length(m.v), length(m.v))
    for (i,edge) in enumerate(m.e)
        w[m.v[edge[1]],m.v[edge[2]]] = weights[i]
    end
    w
end

### Get velocities matrix ###

function get_velocities(m::OpenStreetMapX.MapData,
            class_speeds::Dict{Int,Float64} = OpenStreetMapX.SPEED_ROADS_URBAN)
    @assert length(m.e) == length(m.w.nzval)
    indices = [(m.v[i],m.v[j]) for (i,j) in m.e]
    V = SparseArrays.spzeros(length(m.v), length(m.v))
    for (i,ind) in enumerate(indices)
        V[ind[1],ind[2]] = class_speeds[m.class[i]]/3.6
    end
    V
end

### Extract route from Dijkstra results object ###

function extract_route(dijkstra::LightGraphs.DijkstraState{Float64,Int64}, startIndex::Int, finishIndex::Int)
    route = Int[]
    distance = dijkstra.dists[finishIndex]
    if distance != Inf
        index = finishIndex
        push!(route, index)
        while index != startIndex
            index = dijkstra.parents[index]
            push!(route, index)
        end
    end
    reverse!(route)
    return route, distance
end

### Extract nodes ID's from route object ###

function get_route_nodes(m::OpenStreetMapX.MapData, route_indices::Array{Int64,1})
    route_nodes = Array{Int}(undef,length(route_indices))
    v = Dict{Int,Int}(reverse.(collect(m.v)))
    for n = 1:length(route_nodes)
        route_nodes[n] = v[route_indices[n]]
    end
    return route_nodes
end

### Generate an ordered list of edges traversed in route ###

function route_edges(m::OpenStreetMapX.MapData, route_nodes::Vector{Int})
	e = Array{Int}(undef,length(route_nodes)-1)
	for i = 2:length(route_nodes)
		e[i-1] = findfirst(isequal((route_nodes[i-1], route_nodes[i])), m.e)
	end
	return e
end

### Calculate distance with a given weights ###

calculate_distance(m::OpenStreetMapX.MapData, weights::SparseArrays.SparseMatrixCSC{Float64,Int64}, route_indices::Array{Int64,1}) = sum(weights[(route_indices[i-1], route_indices[i])] for i = 2:length(route_indices))

#####################################
### Find Route with Given Weights ###
#####################################

function find_route(m::OpenStreetMapX.MapData, node0::Int, node1::Int,
                    weights::SparseArrays.SparseMatrixCSC{Float64,Int64};
                    routing::Symbol = :astar, heuristic::Function = (u,v) -> zero(Float64),
                    get_distance::Bool = false, get_time::Bool = false)
    result = Any[]
	start_vertex = m.v[node0]
	finish_vertex= m.v[node1]
    if routing == :dijkstra
        dijkstra_result = OpenStreetMapX.dijkstra(m, weights, start_vertex)
        route_indices, route_values = OpenStreetMapX.extract_route(dijkstra_result, start_vertex, finish_vertex)
        route_nodes = OpenStreetMapX.get_route_nodes(m, route_indices)
        push!(result, route_nodes, route_values)
    elseif routing == :astar
        route_indices, route_values = OpenStreetMapX.a_star_algorithm(m.g, start_vertex, finish_vertex, weights, heuristic)
        route_nodes = OpenStreetMapX.get_route_nodes(m, route_indices)
        push!(result, route_nodes, route_values)
    else
        @warn "routing module declared wrongly - a star algorithm will be used instead!"
        route_indices, route_values = OpenStreetMapX.a_star_algorithm(m.g, start_vertex, finish_vertex, weights, heuristic)
        route_nodes = OpenStreetMapX.get_route_nodes(m, route_indices)
        push!(result, route_nodes, route_values)
    end
    if get_distance
        if isempty(route_indices)
            distance = Inf
        elseif length(route_indices) == 1
            distance = 0
        else
            distance = OpenStreetMapX.calculate_distance(m, m.w, route_indices)
        end
        push!(result, distance)
    end
    if get_time
        w = OpenStreetMapX.create_weights_matrix(m,network_travel_times(m))
        if isempty(route_indices)
            route_time = Inf
        elseif length(route_indices) == 1
            route_time = 0
        else
            route_time = OpenStreetMapX.calculate_distance(m, w, route_indices)
        end
        push!(result, route_time)
    end
    return result
end

#########################################################
### Find Route Connecting 3 Points with Given Weights ###
#########################################################

function find_route(m::OpenStreetMapX.MapData, node0::Int, node1::Int, node2::Int,
                    weights::SparseArrays.SparseMatrixCSC{Float64,Int64};
                    routing::Symbol = :astar, heuristic::Function = (u,v) -> zero(Float64),
                    get_distance::Bool = false, get_time::Bool = false)
    result = Any[]
    route1 = OpenStreetMapX.find_route(m, node0, node1, weights,
                                        routing = routing, heuristic = heuristic,
                                        get_distance = get_distance, get_time = get_time)
    route2 = OpenStreetMapX.find_route(m, node1, node2, weights,
                                        routing = routing, heuristic = heuristic,
                                        get_distance = get_distance, get_time = get_time)
    push!(result,vcat(route1[1],route2[1]))
    for i = 2:length(route1)
        push!(result,route1[i] + route2[i])
    end
    return result
end

###########################
###  ###
###########################
"""
    shortest_route(m::MapData, node1::Int, node2::Int; routing::Symbol = :astar)

Find Shortest route between `node1` and `node2` on map `m`.

"""
function shortest_route(m::MapData, node1::Int, node2::Int; routing::Symbol = :astar)
    route_nodes, distance, route_time = OpenStreetMapX.find_route(m,node1,node2,m.w,
                                                                routing = routing, heuristic = (u,v) -> OpenStreetMapX.get_distance(u, v, m.nodes, m.n),
                                                                get_distance =false, get_time = true)
    return route_nodes, distance, route_time
end
"""
    shortest_route(m::MapData, node1::Int, node2::Int, node3::Int; routing::Symbol = :astar)

Find Shortest route between `node1` and `node2` and `node3` on map `m`.

"""
function shortest_route(m::MapData, node1::Int, node2::Int, node3::Int; routing::Symbol = :astar)
    route_nodes, distance, route_time = OpenStreetMapX.find_route(m,node1,node2, node3, m.w,
                                                                routing = routing, heuristic = (u,v) -> OpenStreetMapX.get_distance(u, v, m.nodes, m.n),
                                                                get_distance =false, get_time = true)
    return route_nodes, distance, route_time
end

"""
    ffastest_route(m::MapData, node1::Int, node2::Int;
                        routing::Symbol = :astar,
                        speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN)

Find fastest route between `node1` and `node2`  on map `m` with assuming `speeds` for road classes.

"""
function fastest_route(m::MapData, node1::Int, node2::Int;
                        routing::Symbol = :astar,
                        speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN)
    w = OpenStreetMapX.create_weights_matrix(m,network_travel_times(m, speeds))
    route_nodes, route_time, distance = OpenStreetMapX.find_route(m, node1, node2, w,
                                                                routing = routing,
																heuristic = (u,v) -> OpenStreetMapX.get_distance(u, v, m.nodes, m.n) / maximum(values(speeds)),
                                                                get_distance = true, get_time = false)
    return route_nodes, distance, route_time
end

"""
    fastest_route(m::MapData, node1::Int, node2::Int, node3::Int;
                        routing::Symbol = :astar,
                        speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN)

Find fastest route between `node1` and `node2` and `node3`  on map `m` with assuming `speeds` for road classes.

"""
function fastest_route(m::MapData, node1::Int, node2::Int, node3::Int;
                        routing::Symbol = :astar,
                        speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN)
    w = OpenStreetMapX.create_weights_matrix(m,network_travel_times(m, speeds))
    route_nodes, route_time, distance = OpenStreetMapX.find_route(m, node1, node2, node3, w,
                                                                routing = routing, heuristic = (u,v) -> OpenStreetMapX.get_distance(u, v, m.nodes, m.n) / maximum(values(speeds)),
                                                                get_distance = true, get_time = false)
    return route_nodes, distance, route_time
end

###########################################
### Find  waypoint minimizing the route ###
###########################################

### Approximate solution ###

function find_optimal_waypoint_approx(m::OpenStreetMapX.MapData, weights::SparseArrays.SparseMatrixCSC{Float64,Int64}, node0::Int, node1::Int, waypoints::Dict{Int,Int})
    dists_start_waypoint = LightGraphs.dijkstra_shortest_paths(m.g, m.v[node0], weights).dists
    dists_waypoint_fin = LightGraphs.dijkstra_shortest_paths(m.g, m.v[node1], weights).dists
    node_id = NaN
    min_dist = Inf
    for (key,value) in waypoints
        dist  = dists_start_waypoint[m.v[value]] + dists_waypoint_fin[m.v[value]]
        if dist < min_dist
            min_dist = dist
            node_id = value
        end
    end
    return node_id
end

### Exact solution ###

function find_optimal_waypoint_exact(m::OpenStreetMapX.MapData, weights::SparseArrays.SparseMatrixCSC{Float64,Int64}, node0::Int, node1::Int, waypoints::Dict{Int,Int})
    dists_start_waypoint = LightGraphs.dijkstra_shortest_paths(m.g, m.v[node0], weights).dists
    node_id = NaN
    min_dist = Inf
    for (key,value) in waypoints
        dist_to_fin = LightGraphs.dijkstra_shortest_paths(m.g, m.v[value], weights).dists[m.v[node1]]
        dist  = dists_start_waypoint[m.v[value]] + dist_to_fin
        if dist < min_dist
            min_dist = dist
            node_id = value
        end
    end
    return node_id
end

########################################################################
### Find Nodes Within Driving Time or Distance - Auxiliary Functions ###
########################################################################

### Bellman Ford's Algorithm ###
function bellman_ford(m::OpenStreetMapX.MapData, w::SparseArrays.SparseMatrixCSC{Float64,Int64}, start_vertices::Vector{Int})
    return LightGraphs.bellman_ford_shortest_paths(m.g, start_vertices, w)
end

### Filter vertices from bellman_fordStates object ###

function filter_vertices(vertices::Dict{Int,Int}, weights::Vector{Float64}, limit::Float64)
    if limit == Inf
        @assert length(vertices) == length(weights)
        return keys(vertices), weights
    end
    indices = Int[]
    distances = Float64[]
    for vertex in keys(vertices)
        distance = weights[vertices[vertex]]
        if distance < limit
            push!(indices, vertex)
            push!(distances, distance)
        end
    end
    return indices, distances
end

##############################################################################
### Extract Nodes from bellman_fordStates Object Within an (Optional) Limit ###
### Based on Weights													   ###
##############################################################################

function nodes_within_weights(m::OpenStreetMapX.MapData, weights::SparseArrays.SparseMatrixCSC{Float64,Int64}, start_indices::Vector{Int}, limit::Float64=Inf)
	start_vertices = [m.v[i] for i in start_indices]
    bellman_ford = OpenStreetMapX.bellman_ford(m, weights, start_vertices)
    return OpenStreetMapX.filter_vertices(m.v, bellman_ford.dists, limit)
end

nodes_within_weights(nodes::Dict{Int,T}, m::OpenStreetMapX.MapData, weights::SparseArrays.SparseMatrixCSC{Float64,Int64}, loc::T, limit::Float64=Inf,locrange::Float64=500.0) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF}) = OpenStreetMapX.nodes_within_weights(m, weights, nodes_within_range(nodes, loc, network, locrange), limit)

##############################################################################
### Extract Nodes from bellman_fordStates Object Within an (Optional) Limit ###
### Based on Driving Distance											   ###
##############################################################################

function nodes_within_driving_distance(m::OpenStreetMapX.MapData, start_indices::Vector{Int}, limit::Float64=Inf)
    start_vertices = [m.v[i] for i in start_indices]
    bellman_ford = OpenStreetMapX.bellman_ford(m, m.w, start_vertices)
    return OpenStreetMapX.filter_vertices(m.v, bellman_ford.dists, limit)
end

nodes_within_driving_distance(nodes::Dict{Int,T}, m::OpenStreetMapX.MapData, loc::T, limit::Float64=Inf,locrange::Float64=500.0) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})= OpenStreetMapX.nodes_within_driving_distance(m, nodes_within_range(nodes, loc ,network, locrange), limit)

##############################################################################
### Extract Nodes from bellman_fordStates Object Within an (Optional) Limit ###
### Based on Driving Time												   ###
##############################################################################

function nodes_within_driving_time(m::OpenStreetMapX.MapData, start_indices::Vector{Int}, limit::Float64=Inf, speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN)
	w = OpenStreetMapX.create_weights_matrix(m,network_travel_times(m, speeds))
	start_vertices = [m.v[i] for i in start_indices]
    bellman_ford = OpenStreetMapX.bellman_ford(m, w, start_vertices)
    return OpenStreetMapX.filter_vertices(m.v, bellman_ford.dists, limit)
end

function nodes_within_driving_time(nodes::Dict{Int,T}, m::OpenStreetMapX.MapData, loc::T, limit::Float64=Inf, locrange::Float64=500.0, speeds::Dict{Int,Float64}=SPEED_ROADS_URBAN) where T<:(Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF})
	w = OpenStreetMapX.create_weights_matrix(m,network_travel_times(m, speeds))
	return OpenStreetMapX.nodes_within_driving_time(m,nodes_within_range(nodes, loc, m,locrange),limit,speeds)
end

"""
    generate_point_in_bounds(m::MapData)

Generates a random pair of Latitude-Longitude coordinates within
boundaries of map `m`
"""
function generate_point_in_bounds(m::MapData)
    boundaries = m.bounds
    (rand() * (boundaries.max_y -  boundaries.min_y) + boundaries.min_y,
    rand() * (boundaries.max_x -  boundaries.min_x) + boundaries.min_x)
end

"""
    point_to_nodes(point::Tuple{Float64,Float64}, m::MapData)

Converts a pair Latitude-Longitude of coordinates
`point` to a node on a map `m`
The result is a node indentifier.

"""
function point_to_nodes(point::Tuple{Float64,Float64}, m::MapData)
    pointLLA = LLA(point[1],point[2])
    point_to_nodes(pointLLA, m)
end
"""
    point_to_nodes(point::LLA, m::MapData)

Converts a pair of coordinates LLA (Latitude-Longitude-Altitude) `point` to a node on a map `m`
The result is a node indentifier.

"""
function point_to_nodes(point::LLA, m::MapData)
    nearest_node(m,OpenStreetMapX.ENU(point, m.bounds))
end


text:
using CUDA.CUDNN: CUDNN_BN_MIN_EPSILON, cudnnBatchNormalizationBackward,
                  cudnnBatchNormalizationForwardInference, CUDNN_BATCHNORM_SPATIAL,
                  cudnnBatchNormalizationForwardTraining

# TODO: replace with new cudnn normalization interface
# https://github.com/JuliaGPU/CUDA.jl/blob/master/lib/cudnn/normalization.jl

mutable struct BNCache
  mean
  ivar
end

BNCache() = BNCache(nothing, nothing)

@inline _wsize(x::AbstractArray{<:Any,N}) where N = ntuple(i -> i == N-1 ? size(x, N-1) : 1, N)

function batchnorm(g::Nothing, b::Nothing, x::DenseCuArray,
                   running_mean, running_var, momentum; kws...)
  affine_sz = _wsize(x)
  g = fill!(similar(x, affine_sz), 1)
  b = fill!(similar(x, affine_sz), 0)
  return batchnorm(g, b, x, running_mean, running_var, momentum; kws...)
end

# NOTE: CuDNN supports only 4D and 5D Tensors for BatchNorm Operations
# so reshape a 2D Tensor into 4D
function batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T,2},
                   running_mean, running_var, momentum; kws...) where T<:CUDNNFloat
  x = reshape(x, 1, 1, size(x, 1), size(x, 2))
  y = batchnorm(g, b, x, running_mean, running_var, momentum; kws...)
  return dropdims(y, dims = (1, 2))
end

function batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::Union{DenseCuArray{T,4},DenseCuArray{T,5}},
                   running_mean, running_var, momentum; kws...) where T<:CUDNNFloat
  cudnnBNForward!(similar(x), g, b, x, running_mean, running_var, momentum; kws...)
end

function cudnnBNForward!(y::DenseCuArray{T}, g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T},
                        running_mean, running_var, momentum;
                        cache = nothing,
                        alpha = T(1), beta = T(0),
                        eps = T(1e-5),
                        training = true,
                        affine = true,
                        track_stats = true) where T<:CUDNNFloat
  dims = _wsize(x)
  if eps < CUDNN_BN_MIN_EPSILON
    @warn "eps $eps is too small for CuDNN, setting to CUDNN_BN_MIN_EPSILON=$CUDNN_BN_MIN_EPSILON"
    eps = CUDNN_BN_MIN_EPSILON
  end

  if running_mean === nothing || running_var === nothing
    running_mean !== running_var && throw(ArgumentError("both or neither of running_mean and running_var must be nothing"))
    if track_stats || !training
      running_mean = fill!(similar(x, dims), 0)
      running_var = fill!(similar(x, dims), 1)
    end
  end

  xd = cudnnTensorDescriptor(x)
  yd = cudnnTensorDescriptor(y)
  gd = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, cudnnDataType(T), Cint(length(dims)), dim4(dims,Val(CUDNN_TENSOR_NCHW)))

  if training
    if !track_stats
      running_mean = CU_NULL
      running_var = CU_NULL
    end

    if cache !== nothing
      mean = fill!(similar(x, dims), 0)
      ivar = fill!(similar(x, dims), 1)
    else
      mean = CU_NULL
      ivar = CU_NULL
    end

    cudnnBatchNormalizationForwardTraining(handle(), CUDNN_BATCHNORM_SPATIAL, scalingParameter(T, alpha), scalingParameter(T, beta), xd, x, yd, y, gd, g, b, momentum, running_mean, running_var, eps, mean, ivar)

    if cache !== nothing
      cache.mean = mean
      cache.ivar = ivar
    end
  else
    cudnnBatchNormalizationForwardInference(handle(), CUDNN_BATCHNORM_SPATIAL, scalingParameter(T, alpha), scalingParameter(T, beta), xd, x, yd, y, gd, g, b, running_mean, running_var, eps)
  end
  return y
end

function ∇batchnorm(g::Nothing, b::Nothing, x::DenseCuArray, dy::DenseCuArray,
                    running_mean, running_var, momentum; kws...)
  affine_sz = _wsize(x)
  g = fill!(similar(x, affine_sz), 1)
  b = fill!(similar(x, affine_sz), 0)
  return ∇batchnorm(g, b, x, dy, running_mean, running_var, momentum; kws...)
end

function ∇batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T, 2}, dy::DenseCuArray{T, 2},
            running_mean, running_var, momentum;
            kws...) where T<:CUDNNFloat
  dg, db, dx = ∇batchnorm(g, b, reshape(x, 1, 1, size(x, 1), size(x, 2)), reshape(dy, 1, 1, size(dy, 1),
                          size(dy, 2)), running_mean, running_var, momentum; kws...)
  (dg, db, dropdims(dx, dims = (1, 2)))
end

function ∇batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                    running_mean, running_var, momentum;
                    affine=true, kws...) where T<:CUDNNFloat
  dg = similar(g)
  db = similar(b)
  dx = similar(x)
  cudnnBNBackward!(dg, g, db, dx, x, dy, running_mean, running_var, T(momentum); kws...)
  if affine
    (dg, db, dx)
  else
    # CUDNN always calculates dg and db, therefore we just have to drop them
    (nothing, nothing, dx)
  end
end

function cudnnBNBackward!(dg::DenseCuArray{T}, g::DenseCuArray{T}, db::DenseCuArray{T},
                          dx::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                          running_mean, running_var,
                          momentum; cache = nothing, eps = T(1e-5),
                          alpha = T(1), beta = T(0),
                          dalpha = T(1), dbeta = T(0), training = true,
                          track_stats = true) where T<:CUDNNFloat
  if !track_stats
    running_mean = CU_NULL
    running_var = CU_NULL
  end

  xd = cudnnTensorDescriptor(x)
  dyd = cudnnTensorDescriptor(dy)
  dxd = cudnnTensorDescriptor(dx)
  gd = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, cudnnDataType(T), Cint(length(_wsize(x))), dim4(_wsize(x),Val(CUDNN_TENSOR_NCHW)))
  if cache !== nothing
    @debug "fetching mean and ivar from the cache"
    mean, ivar = cache.mean, cache.ivar
  else
    mean, ivar = CU_NULL, CU_NULL
  end

  if eps < CUDNN_BN_MIN_EPSILON
    @warn "eps $eps is too small for CuDNN, setting to CUDNN_BN_MIN_EPSILON=$CUDNN_BN_MIN_EPSILON"
    eps = CUDNN_BN_MIN_EPSILON
  end

  cudnnBatchNormalizationBackward(handle(), CUDNN_BATCHNORM_SPATIAL,
        scalingParameter(T, alpha), scalingParameter(T, beta), scalingParameter(T, dalpha), scalingParameter(T, dbeta),
        xd, x, dyd, dy, dxd, dx, gd, g, dg, db, eps, mean, ivar)
end