text:
using ReinforcementLearning, ReinforcementLearningEnvironments, Flux, BenchmarkTools

b = CircularArrayBuffer{Float64}(84, 84, 10)

function show_msg(desc, x)
    println(repeat("=", 50))
    println(desc)
    println(repeat("=", 50))
    display(x)
    println(repeat("=", 50))
    println()
end

show_msg("push a frame into buffer", @benchmark(push!($b, $(rand(84, 84)))))
show_msg("view a frame", @benchmark(select_frame($b, $(rand(1:10)))))
show_msg("consecutive view a batch", @benchmark(consecutive_view($b, $(rand(1:5, 32)), 4)))
show_msg("get an element", @benchmark(getindex($b, $(rand(1:84)), $(rand(1:84)), $(rand(1:10)))))

m = Dense(84*84, 512)

show_msg("view a frame and do matrix mul", @benchmark($m(vec(select_frame($b, $(rand(1:10)))))))
show_msg("view a frame and do matrix mul (use CartesianIndex)", @benchmark($m(vec(view($b, :, :, $(rand(1:10)))))))
show_msg("get a frame and do matrix mul", @benchmark($m(vec(getindex($b, :, :, $(rand(1:10)))))))
show_msg("do matrix mul", @benchmark($m(reshape($b, :, 10))))
show_msg("do matrix mul(for comparison)", @benchmark($m(reshape($(rand(84, 84, 10)), :, 10))))

#=
==================================================
push a frame into buffer
==================================================
BenchmarkTools.Trial:
  memory estimate:  128 bytes
  allocs estimate:  4
  --------------
  minimum time:     1.449 μs (0.00% GC)
  median time:      1.472 μs (0.00% GC)
  mean time:        1.489 μs (0.00% GC)
  maximum time:     7.425 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10
==================================================

==================================================
view a frame
==================================================
BenchmarkTools.Trial:
  memory estimate:  64 bytes
  allocs estimate:  1
  --------------
  minimum time:     8.857 ns (0.00% GC)
  median time:      11.057 ns (0.00% GC)
  mean time:        19.323 ns (22.62% GC)
  maximum time:     4.618 μs (99.41% GC)
  --------------
  samples:          10000
  evals/sample:     999
==================================================

==================================================
consecutive view a batch
==================================================
BenchmarkTools.Trial:
  memory estimate:  4.59 KiB
  allocs estimate:  80
  --------------
  minimum time:     1.331 μs (0.00% GC)
  median time:      1.411 μs (0.00% GC)
  mean time:        1.923 μs (20.34% GC)
  maximum time:     585.943 μs (99.59% GC)
  --------------
  samples:          10000
  evals/sample:     10
==================================================

==================================================
get an element
==================================================
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     3.617 ns (0.00% GC)
  median time:      3.637 ns (0.00% GC)
  mean time:        3.701 ns (0.00% GC)
  maximum time:     21.961 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
==================================================

==================================================
view a frame and do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  31.98 KiB
  allocs estimate:  6
  --------------
  minimum time:     294.195 μs (0.00% GC)
  median time:      318.465 μs (0.00% GC)
  mean time:        330.413 μs (0.57% GC)
  maximum time:     8.884 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
view a frame and do matrix mul (use CartesianIndex)
==================================================
BenchmarkTools.Trial:
  memory estimate:  33.33 KiB
  allocs estimate:  52
  --------------
  minimum time:     354.701 μs (0.00% GC)
  median time:      380.825 μs (0.00% GC)
  mean time:        393.057 μs (0.53% GC)
  maximum time:     8.943 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
get a frame and do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  87.95 KiB
  allocs estimate:  39
  --------------
  minimum time:     342.392 μs (0.00% GC)
  median time:      367.611 μs (0.00% GC)
  mean time:        381.534 μs (1.16% GC)
  maximum time:     8.892 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  316.27 KiB
  allocs estimate:  17
  --------------
  minimum time:     1.307 ms (0.00% GC)
  median time:      1.354 ms (0.00% GC)
  mean time:        1.398 ms (0.90% GC)
  maximum time:     10.018 ms (0.00% GC)
  --------------
  samples:          3547
  evals/sample:     1
==================================================

==================================================
do matrix mul(for comparison)
==================================================
BenchmarkTools.Trial:
  memory estimate:  315.95 KiB
  allocs estimate:  8
  --------------
  minimum time:     1.121 ms (0.00% GC)
  median time:      1.163 ms (0.00% GC)
  mean time:        1.202 ms (1.02% GC)
  maximum time:     10.235 ms (0.00% GC)
  --------------
  samples:          4121
  evals/sample:     1
==================================================
=#


text:
abstract type InterleaveMarker end

struct Iyes <: InterleaveMarker end
struct Ino <: InterleaveMarker end

struct InterleavedImage{T,N,AA1<:AbstractArray{T,N}, AA2<:AbstractArray{T,N}, IMS<:NTuple{N,InterleaveMarker}} <: AbstractArray{T,N}
    oddA::AA1
    evenA::AA2
    imarkers::IMS
end

function InterleavedImage(oddA::AbstractArray{T,N}, evenA::AbstractArray{T,N}, idim::Int=N) where {T,N}
    if size(oddA) != size(evenA)
        error("The two image arrays must be of the same size")
    end
    if idim > N || idim < 1
        error("Interleaved dim out of range")
    end
    ims = ((ifelse(x==idim, Iyes(), Ino()) for x = 1:N)...,)
    return InterleavedImage(oddA, evenA, ims)
end

imarkers(img::InterleavedImage) = img.imarkers
interleaved_dim(A::InterleavedImage) = findfirst(isa.(imarkers(A), Iyes))
oddchild(A::InterleavedImage) = A.oddA
evenchild(A::InterleavedImage) = A.evenA

Base.IndexStyle(::Type{<:InterleavedImage}) = IndexCartesian()

_size(szs::Tuple, curm::Iyes, ims...) = (first(szs)*2, _size(Base.tail(szs), ims...)...)
_size(szs::Tuple, curm::Ino, ims...) = (first(szs), _size(Base.tail(szs), ims...)...)
_size(szs::Tuple{}, ims...) = ()

size(B::InterleavedImage) = _size(size(B.oddA), B.imarkers...)

#returns the appropriate child array and translates the query I
#into an index for that array
function arr_idx(img::InterleavedImage{T,N}, I::Tuple) where {T,N}
    markers = imarkers(img)
    child = _chooseimage(img, I, markers...)
    return child, _arr_idx(markers, I)
end

@inline _chooseimage(img, I, ::Ino, rest...) = _chooseimage(img, Base.tail(I), rest...)
@inline _chooseimage(img, I, ::Iyes, rest...) = return isodd(first(I)) ? oddchild(img) : evenchild(img)
@inline _chooseimage(img, I) = error("no yes markers found")

_arr_idx(::Tuple{}, ::Tuple{}) = ()
_arr_idx(markers, I) = (_arr_idx1(first(markers), first(I)), _arr_idx(Base.tail(markers), Base.tail(I))...)
_arr_idx1(::Iyes, i) = i>>1 + isodd(i)
_arr_idx1(::Ino, i)  = i

function getindex(img::InterleavedImage{T,N}, I::Vararg{Int, N}) where {T,N}
    chosenA, idx = arr_idx(img, (I...,))
    return chosenA[idx...]
end

function setindex!(img::InterleavedImage{T,N}, v, I::Vararg{Int, N}) where {T,N}
    chosenA, idx = arr_idx(img, (I...,))
    chosenA[idx...] = v
end

#utilities for ImageMeta and AxisArray types
#Note: this ditches img2's properites in favor of img1's. Will fix if it turns out to be a problem.
InterleavedImage(img1::ImageMeta{T,N}, img2::ImageMeta{T,N}, idim::Int=N) where {T,N} =
    ImageMeta(InterleavedImage(arraydata(img1), arraydata(img2), idim), properties(img1))
InterleavedImage(img1::AxisArray{T,N}, img2::AxisArray{T,N}, idim::Int=N) where {T,N} =
    match_axisspacing(InterleavedImage(arraydata(img1), arraydata(img2), idim), img1)

ext: jl


text:
# Copyright (c) 2015 Rex Computing and Isaac Yonemoto"
# see LICENSE.txt
# this work was supported in part by DARPA Contract D15PC00135

#an i64array made for fsizesize.  Helps multiple dispatch and generated functions
#optimize for certain complex array operations.

__cell_length(FSS) = 1 << (FSS - 6)

type I64Array{FSS}
  a::Array{UInt64,1}
  function I64Array(a::Array{UInt64,1})
    FSS < 7 && throw(ArgumentError("invalid FSS == $FSS < 7"))
    _al = __cell_length(FSS)
    length(a) != _al && throw(ArgumentError("invalid array length, should be $_al != $(length(a))"))
    new(a)
  end
end

Base.rand{FSS}(::Type{I64Array{FSS}}) = I64Array{FSS}(rand(UInt64, __cell_length(FSS)))

#performs a carried add on an unsigned integer array, the old-fashioned way.
function __carried_add_dumb(carry::UInt64, v1::I64Array, v2::I64Array)
  #first perform a direct sum on the integer arrays
  res = v1.a + v2.a
  #iterate downward from the least significant word
  for idx = length(v1.a):-1:2
    #if it looks like it's lower than it should be, then make it okay.
    @inbounds (res[idx] <= v1.a[idx]) && (res[idx - 1] += 1)
  end
  #add the last word with a carry.
  @inbounds carry += res[1] < v1.a[1] ? 0 : 1
  (carry, res)
end

#performs a carried add on an unsigned integer array, using type dispatch
function __carried_add_td{FSS}(carry::UInt64, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)
  res = v1.a + v2.a
  for idx = l:-1:2
    @inbounds (res[idx] <= v1.a[idx]) && (res[idx - 1] += 1)
  end
  @inbounds carry += res[1] < v1.a[1] ? 0 : 1
  (carry, res)
end

#performs a carried add on an unsigned integer array, using a generated function
@generated function __carried_add_gn{FSS}(carry::UInt64, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)

  loop = :()
  for idx = l:-1:2
    loop = :($loop; @inbounds (res[$idx] <= v1.a[$idx]) && (res[$idx - 1] += 1))
  end

  quote
    res = v1.a + v2.a

    $loop

    @inbounds carry += (res[1] < v1.a[1]) ? 0 : 1
    (carry, res)
  end
end

#doubly unroll this guy
@generated function __carried_add_du{FSS}(carry::UInt64, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)
  ex = :(res = zeros(UInt64, $l))
  for idx = l:-1:1
    ex = :($ex; @inbounds res[$idx] = v1.a[$idx] + v2.a[$idx])
  end
  for idx = l:-1:2
    ex = :($ex; @inbounds (res[$idx] <= v1.a[$idx]) && (res[$idx - 1] += 1))
  end
  ex = :($ex; @inbounds carry += res[1] < v1.a[1] ? 0 : 1)
  return :($ex; (carry, res))
end

const PREALLOC_ARRAY = zeros(UInt64, 32)

@generated function __carried_add_dupa{FSS}(carry::UInt64, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)
  ex = :(@simd for idx = 1:$l
           @inbounds PREALLOC_ARRAY[idx] = v1.a[idx] + v2.a[idx]
         end)
  for idx = l:-1:2
    ex = :($ex; @inbounds (PREALLOC_ARRAY[$idx] <= v1.a[$idx]) && (PREALLOC_ARRAY[$idx - 1] += 1))
  end
  ex = :($ex; @inbounds carry += PREALLOC_ARRAY[1] < v1.a[1] ? 0 : 1)
  return :($ex; (carry, slice(PREALLOC_ARRAY, 1:$l)))
end

@generated function __carried_add_res{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)
  loop = :()
  for idx = l:-1:2
    loop = :($loop; @inbounds (res.a[$idx] <= v1.a[$idx]) && (res.a[$idx - 1] += 1))
  end

  quote
    res.a = v1.a + v2.a

    $loop

    @inbounds carry += (res.a[1] < v1.a[1]) ? 0 : 1
    carry
  end
end

@generated function __carried_add_res2{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)

  loop1 = :()
  for idx = l:-1:1
    loop1 = :($loop1; @inbounds res.a[$idx] = v1.a[$idx] + v2.a[$idx] )
  end

  loop2 = :()
  for idx = l:-1:2
    loop2 = :($loop2; @inbounds (res.a[$idx] <= v1.a[$idx]) && (res.a[$idx - 1] += 1))
  end

  quote
    $loop1

    $loop2

    @inbounds (res.a[1] < v1.a[1]) && (carry += 1)
    carry
  end
end

@generated function __carried_add_lahf{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)

  loop1 = :()
  for idx = l:-1:1
    loop1 = :($loop1; @inbounds carry = ccall( (:carried_add_lahf, "./bin/libcarry.so"), UInt64, (UInt64, Ptr{UInt64}, Ptr{UInt64}, Ptr{UInt64}),
      carry,
      pointer(res.a, $idx),
      pointer(v1.a, $idx),
      pointer(v2.a, $idx)))
  end

  loop1
end

@generated function __carried_add_jnc{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)

  loop1 = :()
  for idx = l:-1:1
    loop1 = :($loop1; @inbounds carry = ccall( (:carried_add_jnc, "./bin/libcarry.so"), UInt64, (UInt64, Ptr{UInt64}, Ptr{UInt64}, Ptr{UInt64}),
      carry,
      pointer(res.a, $idx),
      pointer(v1.a, $idx),
      pointer(v2.a, $idx)))
  end

  loop1
end

@generated function __carried_add_loop{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)

  quote
    ccall((:carried_add_loop2, "./bin/libcarry.so"), UInt64, (UInt64, Ptr{UInt64}, Ptr{UInt64}, Ptr{UInt64}, Int64),
    carry,
    pointer(res.a, $l),
    pointer(v1.a, $l),
    pointer(v2.a, $l),
    $l)
  end
end

@generated function __carried_add_asm_gen{FSS}(carry::UInt64, res::I64Array{FSS}, v1::I64Array{FSS}, v2::I64Array{FSS})
  l = __cell_length(FSS)
  fnname = string("carried_add_", FSS)

  quote
    ccall(($(fnname), "./bin/libcarry.so"), UInt64, (UInt64, Ptr{UInt64}, Ptr{UInt64}, Ptr{UInt64}),
    carry,
    pointer(res.a, $l),
    pointer(v1.a, $l),
    pointer(v2.a, $l))
  end
end

const mFSS = 11

a = rand(I64Array{mFSS})
b = rand(I64Array{mFSS})

answer = __carried_add_dumb(zero(UInt64), a, b)[2]
answer != __carried_add_td(zero(UInt64), a, b)[2] && println("error, td")
answer != __carried_add_gn(zero(UInt64), a, b)[2] && println("error, gn")
answer != __carried_add_du(zero(UInt64), a, b)[2] && println("error, du")
answer != __carried_add_dupa(zero(UInt64), a, b)[2] && println("error, dupa")

res = I64Array{mFSS}(zeros(UInt64, __cell_length(mFSS)))
__carried_add_lahf(zero(UInt64), res, a, b)
answer != res.a && println("error, lahf")

res = I64Array{mFSS}(zeros(UInt64, __cell_length(mFSS)))
__carried_add_jnc(zero(UInt64), res, a, b)
answer != res.a && println("error, jnc")

res = I64Array{mFSS}(zeros(UInt64, __cell_length(mFSS)))
__carried_add_loop(zero(UInt64), res, a, b)
answer != res.a && println("error, loop")

res = I64Array{mFSS}(zeros(UInt64, __cell_length(mFSS)))
__carried_add_asm_gen(zero(UInt64), res, a, b)
answer != res.a && println("error, asmgen")

const count = 1000000
const alist = I64Array{mFSS}[rand(I64Array{mFSS}) for idx = 1:count]
const blist = I64Array{mFSS}[rand(I64Array{mFSS}) for idx = 1:count]

function testfn(f)
  z = zero(UInt64)
  for idx = 1:count
    (carry, r) = f(z, alist[idx], blist[idx])
  end
end

function testfncache(f)
  z = zero(UInt64)
  res = I64Array{mFSS}(zeros(UInt64,__cell_length(mFSS)))
  for idx = 1:count
    carry = f(z, res, alist[idx], blist[idx])
  end
end

const GLOBAL_RES = I64Array{mFSS}(zeros(UInt64,__cell_length(mFSS)))
function testfnglobal(f)
  z = zero(UInt64)
  for idx = 1:count
    carry = f(z, GLOBAL_RES, alist[idx], blist[idx])
  end
end

println("----")
println("warning:  results may suffer compiler penalty")
@time testfn(__carried_add_dumb)
@time testfn(__carried_add_td)
@time testfn(__carried_add_gn)
@time testfn(__carried_add_du)
@time testfn(__carried_add_dupa)
@time testfncache(__carried_add_res)
@time testfnglobal(__carried_add_res)
@time testfncache(__carried_add_res2)
@time testfnglobal(__carried_add_res2)
@time testfncache(__carried_add_lahf)
@time testfnglobal(__carried_add_lahf)
@time testfncache(__carried_add_jnc)
@time testfnglobal(__carried_add_jnc)
@time testfncache(__carried_add_loop)
@time testfnglobal(__carried_add_loop)
@time testfncache(__carried_add_asm_gen)
@time testfnglobal(__carried_add_asm_gen)

println("\n----")
println("post-JITted")
print("dumb carried add                    ")
@time testfn(__carried_add_dumb)
print("with semi-smart type dispatch       ")
@time testfn(__carried_add_td)
print("as a generated fn with unrolling    ")
@time testfn(__carried_add_gn)
print("doubly-unrolled generated funtion   ")
@time testfn(__carried_add_du)
print("doubly-unrolled, dumb preallocation ")
@time testfn(__carried_add_dupa)
println("\nsingle unroll function call preallocation")
print("local                               ")
@time testfncache(__carried_add_res)
print("global                              ")
@time testfnglobal(__carried_add_res)
println("\ndouble unroll function call preallocation")
print("local                               ")
@time testfncache(__carried_add_res2)
print("global                              ")
@time testfnglobal(__carried_add_res2)
println("\nusing unrolled assembly (lahf)")
print("local                               ")
@time testfnglobal(__carried_add_lahf)
print("global                              ")
@time testfnglobal(__carried_add_lahf)
println("\nusing unrolled assembly (jnc)")
print("local                               ")
@time testfnglobal(__carried_add_jnc)
print("global                              ")
@time testfnglobal(__carried_add_jnc)
println("\nusing c-looped assembly (lahf)")
print("local                               ")
@time testfnglobal(__carried_add_loop)
print("global                              ")
@time testfnglobal(__carried_add_loop)
println("\nusing python-unrolled assembly")
print("local                               ")
@time testfnglobal(__carried_add_asm_gen)
print("global                              ")
@time testfnglobal(__carried_add_asm_gen)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Test, GpABC, OrdinaryDiffEq, Distances, Distributions
import Random

@testset "Model selection test" begin

	Random.seed!(4)
	threshold_schedule = [20.0, 15.0, 10.0]
	max_iter = 1e4
	n_particles = 200

	#
	# Experimental data - from ABCSysBio example at
	# https://github.com/jamesscottbrown/abc-sysbio/blob/master/examples/SBML/Example1/SIR.ipynb
	#
	times = [0.0, 0.6, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0]
	data = [[ 20.     ,  10.     ,   0.     ],
	       [  0.12313,  13.16813,   9.42344],
	       [  0.12102,   7.17251,  11.18957],
	       [  0.09898,   2.36466,  10.0365 ],
	       [  0.37887,   0.92019,   6.87117],
	       [  1.00661,   0.61958,   4.44955],
	       [  1.20135,   0.17449,   3.01271],
	       [  1.46433,   0.28039,   1.76431],
	       [  1.37789,   0.0985 ,   1.28868],
	       [  1.57073,   0.03343,   0.81813],
	       [  1.4647 ,   0.28544,   0.52111],
	       [  1.24719,   0.10138,   0.22746],
	       [  1.56065,   0.21671,   0.19627]]

	data = hcat(data...)

	#
	# Define priors and initial conditions for each model
	#

	#
	# Priors and initial conditions - these are model-specfic as each model can
	# have different numbers of parameters/species
	#
	priors1 = [Uniform(0.2, 5.0) for i in 1:4]
	ic1 = [20.0, 10.0, 0.0]

	priors2 = vcat([Uniform(0.2, 5.0) for i in 1:4], Uniform(0.2, 10.0))
	ic2 = [20.0, 0.0, 10.0, 0.0]

	priors3 = [Uniform(0.2, 5.0) for i in 1:6]
	ic3 = [20.0, 10.0, 0.0]

	priors = [priors1, priors2, priors3]
	modelprior = DiscreteUniform(1, 3)

	# p = (alpha, gamma, d, v)
	# x = (S, I, R)
	function model1(dx, x, p, t)
	    dx[1] = p[1] - p[2]*x[1]*x[2] - p[3]*x[1] # dS/dt = alpha - gamma*S*I - d*S
	    dx[2] = p[3]*x[1]*x[2] - p[4]*x[2] - p[3]*x[2] # dI/dt = gamma*S*I - v*I - d*I
	    dx[3] = p[4]*x[2] - p[3]*x[3] # dR/dt = v*I - d*R
	end

	# p = (alpha, gamma, d, v, delta)
	# x = (S, L, I, R)
	function model2(dx, x, p, t)
	    dx[1] = p[1] - p[2]*x[1]*x[3] - p[3]*x[1] # dS/dt = alpha - gamma*S*I - d*S
	    dx[2] = p[2]*x[1]*x[3] - p[5]*x[2] - p[3]*x[2] # dL/dt = gamma*S*I - delta*L - d*L
	    dx[3] = p[5]*x[2] - p[4]*x[3] - p[3]*x[3] # dI/dt = delta*L - v*I - d*I
	    dx[4] = p[4]*x[3] - p[3]*x[4] # dR/dt = v*I - d*R
	end

	ics = [ic1, ic2, ic3]

	# Define simulator functions for each model

	simulator1(params) = Array{Float64,2}(
	    solve(ODEProblem(model1, ics[1], (times[1], times[end]), params), RK4(); saveat=times, force_dtmin=true))

	# Model2 contains the species L, which is not measured - we remove it from the returned ODE solution
	# so that it can be compared to the reference data "data", which only contains S, I and R
	simulator2(params) = Array{Float64,2}(
	    solve(ODEProblem(model2, ics[2], (times[1], times[end]), params), RK4(); saveat=times, force_dtmin=true))[[1,3,4],:]

	# Model to test dead behaviour
	simulator3(params) = rand(Float64, size(data))

	simulators = [simulator1, simulator2, simulator3]

	#
	# For tests on output shapes - the number of parameters in each model
	#
	n_params = [4,5,6]

	function test_ms_output(ms_res::ModelSelectionOutput, is_simulation::Bool)
		# 2nd dimension of population should always be the number of parameters
		@test all(vcat([[size(ms_res.smc_outputs[m].population[i],2) == n_params[m]
			for i in 1:length(ms_res.threshold_schedule)]
					for m in 1:ms_res.M]...))

		# First dimension of population should always be the number of accepted particles
		@test all(vcat([[size(ms_res.smc_outputs[m].population[i],1) == ms_res.n_accepted[i][m]
			for i in 1:length(ms_res.threshold_schedule)]
				for m in 1:ms_res.M]...))

		# First dimension of distances should always be the number of accepted particles
		@test all(vcat([[size(ms_res.smc_outputs[m].distances[i],1) == ms_res.n_accepted[i][m]
			for i in 1:length(ms_res.threshold_schedule)]
					for m in 1:ms_res.M]...))

		# First dimension of weights should always be the number of accepted particles
		@test all(vcat([[size(ms_res.smc_outputs[m].weights[i],1) == ms_res.n_accepted[i][m]
			for i in 1:length(ms_res.threshold_schedule)]
					for m in 1:ms_res.M]...))

		# Total number of accepted particles at each population can't be more than n_particles
		@test all([sum(arr) <= n_particles for arr in ms_res.n_accepted])

		# There should be as many populations, distances as weights as populations
		@test all([size(ms_res.smc_outputs[m].population,1) == length(ms_res.threshold_schedule) for m in 1:ms_res.M])
		@test all([size(ms_res.smc_outputs[m].distances,1) == length(ms_res.threshold_schedule) for m in 1:ms_res.M])
		@test all([size(ms_res.smc_outputs[m].weights,1) == length(ms_res.threshold_schedule) for m in 1:ms_res.M])

		# Can't have more than max_iter tries in total at each population - this is only true for simulation
		if is_simulation
			@test all([sum([ms_res.smc_outputs[m].n_tries[i] for m = 1:ms_res.M]) <= max_iter for i=1:length(ms_res.threshold_schedule)])
		end

		# Weights must sum to 1 for all models in all cases where the model accepted at least one particle
		weight_sums = vcat([[sum(ms_res.smc_outputs[m].weights[i]) for i = 1:length(ms_res.threshold_schedule)] for m in 1:ms_res.M]...)
		nonzero_weights = vcat([[size(ms_res.smc_outputs[m].weights[i],1) > 0 for i = 1:length(ms_res.threshold_schedule)] for m in 1:ms_res.M]...)
		@test all(weight_sums[nonzero_weights] .== 1.0)

		# Check that dead models are properly dead
		zero_acceptance_populations = [findall(ms_res.smc_outputs[m].n_accepted .== 0) for m in 1:ms_res.M]
		for m in 1:ms_res.M
			if size(zero_acceptance_populations[m],1) > 1
				first_dead_pop = zero_acceptance_populations[m][2]
				@test all(ms_res.smc_outputs[m].n_accepted[first_dead_pop:end] .== 0)
				@test all(ms_res.smc_outputs[m].n_tries[first_dead_pop:end] .== 0)
				@test all([pop == zeros(0,ms_res.smc_outputs[m].n_params) for pop in ms_res.smc_outputs[m].population[first_dead_pop:end]])
				@test all([dist == zeros(0) for dist in ms_res.smc_outputs[m].distances[first_dead_pop:end]])
				@test all([we == zeros(0) for we in ms_res.smc_outputs[m].weights[first_dead_pop:end]])
			end
		end

		# All SMC outputs should have maching threshold schedules
		@test all([ms_res.threshold_schedule == ms_res.smc_outputs[m].threshold_schedule for m in 1:ms_res.M])
	end

	#
	# Do model selection using simulation
	ms_res = SimulatedModelSelection(data,
		simulators,
		priors,
		threshold_schedule,
		n_particles)
	test_ms_output(ms_res, true)

	# Same with emulation; no re-training
	n_design_points = 200
	ms_res = EmulatedModelSelection(data,
		simulators,
		priors,
		threshold_schedule,
		n_particles,
		n_design_points)
	test_ms_output(ms_res, false)

	# With emulator re-training on
	ms_res = EmulatedModelSelection(data,
		simulators,
		priors,
		threshold_schedule,
		n_particles,
		n_design_points;
		emulator_retraining=PreviousPopulationThresholdRetraining(200, 50, 10),
		emulated_particle_selection=MeanVarEmulatedParticleSelection())
	test_ms_output(ms_res, false)

	# Repeat above with thresholds that are too small for any particles to be accepted. Algorithm will terminate early
	ms_res = SimulatedModelSelection(data,
		simulators,
		priors,
		[0.4, 0.2], # threshold_schedule,
		n_particles)
	@test !ms_res.completed_all_populations

	ms_res = EmulatedModelSelection(data,
		simulators,
		priors,
		[0.4, 0.2], # threshold_schedule,
		n_particles,
		n_design_points)
	@test !ms_res.completed_all_populations
end


text:
using IterTools
"""
`wraprange(start, stop, wrapat) = wraprange(start, 1, stop, wrapat)`
`wraprange(start, step, stop, wrapat)`

a range that wraps around back to 1, at `wrapat`.
Can be used for a sliding window index that wraps around when it reaches the size
of the dimension it's indexing

e.g.
```
julia> collect(wraprange(4,3,7))
7-element Array{Int64,1}:
 4
 5
 6
 7
 1
 2
 3

julia> collect(wraprange(4, 2, 3, 7))
4-element Array{Int64,1}:
 4
 6
 1
 3

julia> a = rand(2,5)
2×5 Array{Float64,2}:
 0.79923  0.280401  0.598733  0.213387  0.94785
 0.69775  0.246508  0.981085  0.318419  0.0183561

# return the view a[:,4:5] and a[:, 1:2]
julia> view(a, :, wraprangeidxs(4,2,5))
2×4 SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Array{Int64,1}},false}:
 0.213387  0.94785    0.79923  0.280401
 0.318419  0.0183561  0.69775  0.246508
```
"""
wraprange(start, stop, wrapat) = wraprange(start, 1, stop, wrapat)
wraprange(start, step, stop, wrapat) = begin
    stop < start || return StepRange(start, step, stop)
    start <= wrapat ||
        throw(ArgumentError(
            "start ($start) must be smaller than wrapat ($wrapat)"))
    firstidx = step - (wrapat - start)%step
    chain(start:step:wrapat, firstidx:step:stop)
end

"""
`wraprangeidxs(args...)`

convenience for `collect(wraprange(args...))`
"""
wraprangeidxs(args...) = collect(wraprange(args...))

colons(n) = fill(Colon(), n)


text:
#
# This script is used to update the Julia wrapper packages for GAP packages.
#

# input GAP package URL
# output: generate / update
using ArtifactUtils
using Downloads
using GAP
using Pkg

# GAP.jl version to require as compat
GAP_jl_version = "v0.7.2"

"""
    parse_pkginfo(pkginfopath::String)

Given the path `pkginfopath` to a GAP `PackageInfo.g` file, read it (or as
much of it as possible, anyway) into a Julia dict which is returned.
"""
function parse_pkginfo(pkginfopath::String)
    isfile(pkginfopath) || error("$pkginfopath is not a regular file")
    GAP.Globals.Read(GapObj(pkginfopath))
    r = GAP.Globals.GAPInfo.PackageInfoCurrent
    @assert GAP.Globals.IsRecord(r)

    # convert the record, but skip AvailabilityTest or anything else bound to a function
    pkginfo = Dict{Symbol,Any}()
    for key in Vector{Symbol}(GAP.Globals.RecNames(r))
        val = getproperty(r, key)
        GAP.Globals.IsFunction(val) && continue
        pkginfo[key] = GAP.gap_to_julia(val; recursive = true)
    end
    return pkginfo
end

"""
    lookup_pkg(pkgname::String)

Given the package name `pkgname`, return either its UUID (if a package with
the name exists), or `nothing`
"""
function lookup_pkg(pkgname::String)
    @static if VERSION >= v"1.7-DEV"
        uuids = Pkg.Types.registered_uuids(Pkg.Types.Context().registries, pkgname)
    else
        uuids = Pkg.Types.registered_uuids(Pkg.Types.Context(), pkgname)
    end
    isempty(uuids) && return nothing
    @assert length(uuids) == 1
    return first(uuids)
end

# given the path to a GAP package directory,
function update_pkg(pkgpath::String)

    # read PackageInfo.g
    isdir(pkgpath) || error("$path is not a directory")
    pkginfopath = joinpath(pkgpath, "PackageInfo.g")
    pkginfo = parse_pkginfo(pkginfopath)

    # get the names of the GAP package and the corresponding Julia package
    gap_pkgname = pkginfo[:PackageName]
    pkgname = lowercase(gap_pkgname)
    julia_pkgname = "GAP_pkg_$(pkgname)"

    # ensure there is a package directory
    if !isdir(julia_pkgname)
        Pkg.generate(julia_pkgname)
    end

    #
    # Update Project.toml
    #
    projpath = "$(julia_pkgname)/Project.toml"
    proj = Pkg.Types.read_project(projpath)
    proj.deps["Pkg"] = Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f")
    proj.deps["GAP"] = lookup_pkg("GAP")
    proj.compat["GAP"] = "0.7" # FIXME adjust this sometimes? but when / who does it?
    proj.compat["julia"] = "1.6"

    println("TODO: adjust version...")

    #
    # check if there is a corresponding JLL; if so, add it as a dependency
    #
    jllname = "$(julia_pkgname)_jll"
    jlluuid = lookup_pkg(jllname)
    # TODO: can we automatically "guess" the corresponding version
    # of the JLL?? maybe in some cases, but not always...

    # translate GAP package dependencies
    # TODO: beware of cycles, e.g. polycyclic vs. polenta
    deps_statements = ""
#=
    for dep in pkginfo[:Dependencies][:NeededOtherPackages]
        # vectors of the form ["GAPDoc", "1.5"]
        lowercase(dep[1]) == "gapdoc" && continue
        depname = "GAP_pkg_$(lowercase(dep[1]))"
        uuid = lookup_pkg(depname) # TODO: might not yet be in the global registry; so perhaps we should also load all UUIds from the local Project.toml files?
        uuid !== nothing || error("could not find $depname")
        proj.deps[depname] = uuid
        # TODO: compat???
        deps_statements *= "using $depname\n"
    end
=#

    # Write out Project.toml
    Pkg.Types.write_project(proj, projpath)

    #
    # generate README.md
    #
    open(joinpath(julia_pkgname, "README.md"), "w") do f
       write(f,"""
       # $(julia_pkgname)

       This is generated wrapper package for the GAP package $(gap_pkgname) $(pkginfo[:Version]).
       """)
   end

    #
    # all the Julia wrapper packages are MIT licensed, like this script
    #
    cp("LICENSE", joinpath(julia_pkgname, "LICENSE"), force=true)

    #
    # add the package tarball as an artifact
    #
    url = pkginfo[:ArchiveURL]
    types = split(pkginfo[:ArchiveFormats], " ")
    if ".tar.bz2" in types
        url *= ".tar.bz2"
    elseif ".tar.gz" in types
        url *= ".tar.gz"
    else
        url *= first(types)
    end
    artifacts_toml = joinpath(julia_pkgname, "Artifacts.toml")
    ArtifactUtils.add_artifact!(artifacts_toml, pkgname, url; force=true)

    #
    # generate source code
    #
    mkpath("$(julia_pkgname)/src")

    jlluuid !== nothing && open(joinpath(julia_pkgname, "src", "jll.jl"), "w") do f
        write(f,"""
        module JLL

        using GAP
        using $jllname

        function __init__()
            # ensure GAP finds kernel extensions or other binaries
            sopath = joinpath($jllname.find_artifact_dir(), "lib", "gap")
            @debug "GAP package '$pkgname' sopath = " * sopath
            setproperty!(GAP.Globals.DirectoriesPackageProgramsOverrides, :$pkgname, GapObj(sopath))
        end

        end # module
        """)
    end

    if jlluuid !== nothing
        deps_statements *= """

        include("jll.jl")
        """
    end

    open(joinpath(julia_pkgname, "src", "$(julia_pkgname).jl"), "w") do f
        write(f,"""
        # ATTENTION: This file was generated, do not edit by hand.
        module $julia_pkgname

        using Pkg.Artifacts
        using GAP
        $deps_statements
        function __init__()
            # ensure GAP will load the correct instance of this package
            path = artifact"$pkgname"
            @debug "GAP package '$pkgname' path = " * path
            GAP.Globals.SetPackagePath(GapObj("$pkgname"), GapObj(path))
        end

        function load()
            res = GAP.Globals.LoadPackage(GapObj("$pkgname"))
            if res != true
                error("failed to load GAP package $pkgname")
            end
        end
        end # module
        """)
    end

    #
    # generate a simple test
    #
    mkpath("$(julia_pkgname)/test")
    open(joinpath(julia_pkgname, "test", "runtests.jl"), "w") do f
        write(f,"""
        using $julia_pkgname
        using Test

        import $julia_pkgname
        println("TODO")
        """)
    end

    return nothing
end


text:
##=================####==============##

# Different update schemes

##=================####==============##

"""
    update_state!(agents)
Finds agent's next state and updates it, it just packages
`get_next_state!` and `update_state!`
"""
function update_single_state!(ag)
    push!(ag.previous, ag.state)
    #  if (ag.state == "S" && ag.new_state == "I")
    if ag.new_state == "I"
        ag.counter = ag.counter + 1
    end
    ag.state = ag.new_state
end

function update_state!(agents)
    Threads.@threads for ag in agents
        update_single_state!(ag)
    end
end

##=================####==============##

function update_coop!(agents,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            if length(ag.coopf)/ag.degree_t > threshold
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

##=================####==============##

function update_single_effect_distance_coop_stochastic!(agents,g,d,threshold,steps;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    the_final = 0
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = abs(curr_effect - mean_effect)
        the_change = the_diff*steps[dist]
        the_sign = rand(di) |> rand_sign
#         the_sign = rand() |> rand_sign
        if rand() <= 1/dist
            if agents[v].attitude == "ra"
                the_final = curr_effect - the_change * the_sign
            elseif agents[v].attitude == "rt"
                the_final = curr_effect + the_change * the_sign
            end
        end
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    agents[v].new_coop_effect = the_final
end

##=================####==============##

function update_effect_given_distance_coop_stochastic!(agents,g,d,threshold,step)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_coop_stochastic!(agents,g,d,threshold,step;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

function update_single_effect_distance_coop_prior!(agents,g,d,threshold,steps,sd_change;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    the_final = 0
    di = truncated(Normal(agents[v].prior_bel,sd_change),0,1)
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = abs(curr_effect - mean_effect)
        the_change = the_diff*steps[dist]
        the_sign = rand(di) |> rand_sign
        if rand() <= 1/dist
            if agents[v].attitude == "ra"
                the_final = curr_effect - (the_change * the_sign)
            elseif agents[v].attitude == "rt"
                the_final = curr_effect + (the_change * the_sign)
            end
        end
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    agents[v].new_coop_effect = the_final
end

##=================####==============##

function update_effect_given_distance_coop_prior!(agents,g,d,threshold,step,sd_change)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_coop_prior!(agents,g,d,threshold,step,sd_change;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

function update_coop_distance!(agents,g,d,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            status = agents[neighborhood(g,ag.id,d)|> unique] |> f -> map(x -> x.at_home,f)
            thecoop = sum(status)
            tot = length(status)
            if thecoop/tot >= threshold
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

##=================####==============##

function prev_infected(agent,threshold)
    historic = agent.previous[end-threshold:end]
    Infected = filter(x->x=="I",historic)
    Susceptible = filter(x->x=="S",historic)
    Recovered = filter(x->x=="R",historic)
    Infected, Susceptible, Recovered
end

##=================####==============##

function update_coop_infections!(agents,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            Infected, Susceptible, Recovered = prev_infected(ag, threshold)
            if length(Infected) >= 1
                ag.at_home = true
                #  ag.at_home = !(ag.at_home)
            else
                if lrt
                    ag.at_home = false
                    #  ag.at_home = !(ag.at_home)
                end
            end
        end
    end
end

##=================####==============##

function update_coop_distance_inf!(agents,g,d,threshold;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            status = agents[neighborhood(g,ag.id,d)|> unique] |> f -> map(x -> prev_infected(x,threshold)[1],f)
            theinfected = sum(map(x -> "I" in x, status))
            tot = length(status)
            if theinfected >= 1
                ag.at_home = true
            else
                if lrt
                    ag.at_home = false
                end
            end
        end
    end
end

##=================####==============##
##### Check infections and change behavior depending on distance

function update_single_given_distance!(agents,g,v,d,threshold,probs)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)
    distances = last.(neigh)
    changed = false
    for dist in 0:d
        current = findall(x->x==dist,distances)
        infected = findall(x->x.counter >= threshold, agents[current])
        for time in 1:length(infected)
            if rand() < probs[dist+1]
                agents[v].at_home = true
                changed = true
                break
            end
        end
        if changed
            break
        end
    end
end

##=================####==============##

function update_coop_given_distance!(agents,g,d,threshold,probs;lrt=false)
    Threads.@threads for ag in agents
        if ag.adapter
            update_single_given_distance!(agents,g,ag.id,d,threshold,probs)
        end
    end
end

##=================####==============##

function update_single_effect_distance!(agents,g,d,threshold,step;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)
    distances = last.(neigh)
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        infected = findall(x->x.counter >= threshold, agents[current_nodes])
        for time in 1:length(infected)
            if agents[v].attitude == "ra"
                agents[v].coop_effect = max(agents[v].coop_effect-step[dist+1],0)
            elseif agents[v].attitude == "rt"
                agents[v].coop_effect = min(agents[v].coop_effect+step[dist+1],1)
            end
        end
    end
end

##=================####==============##

function update_effect_given_distance!(agents,g,d,threshold,step)
    Threads.@threads for ag in agents
        if ag.adapter
            update_single_effect_distance!(agents,g,d,threshold,step;v=ag.id)
        end
    end
end

##=================####==============##

function update_single_effect_distance_coop!(agents,g,d,threshold,steps;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    the_final = 0
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = abs(curr_effect - mean_effect)
        the_change = the_diff*steps[dist]
        if rand() <= 1/dist
            if agents[v].attitude == "ra"
                the_final = curr_effect - the_change
            elseif agents[v].attitude == "rt"
                the_final = curr_effect + the_change
            end
        end
    end
    if agents[v].attitude == "ra"
        agents[v].new_coop_effect = max(the_final,0)
    elseif agents[v].attitude == "rt"
        agents[v].new_coop_effect = min(the_final,1)
    end
end

##=================####==============##

function update_single_effect_distance_coop_free_dir!(agents,g,d,threshold,steps;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    the_final = 0
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = (curr_effect - mean_effect)
        the_change = the_diff*steps[dist]
        if rand() <= 1/dist
            the_final = curr_effect - the_change
        end
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    agents[v].new_coop_effect = the_final
end
##=================####==============##

function update_effect_given_distance_coop!(agents,g,d,threshold,step)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_coop!(agents,g,d,threshold,step;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

function update_effect_given_distance_coop_free_dir!(agents,g,d,threshold,step)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_coop_free_dir!(agents,g,d,threshold,step;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

#  """
#      update_state!(agent)
#  Updates `agent.state` with the one computed in `get_next_state!`
#  """
#  function update_state!(agent)
#      agent.state = agent.new_state
#  end

##=================####==============##

"""
    update_all_agents!(agents, now_t)
Finds agent's next state and updates it, it just packages
`get_next_state!` and `update_state!`
"""
function update_all_agents!(agents, params)
    # FIND AGENTS' NEXT STATE FROM INTERACTIONS
    Threads.@threads for ag_id in 1:params.num_agents
        get_next_state!(agents, ag_id, params)
    end

    #  Threads.@threads for ag in agents
    #      get_next_state!(ag, now_t)
    #  end

    # UPDATE AGENTS' STATE
    Threads.@threads for ag in agents
        update_state!(ag)
    end
end

##=================####==============##

function update_single_effect_distance_alpha!(agents,g,d,threshold,steps,sd_change,sd_alpha;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    the_final = 0
    di = truncated(Normal(agents[v].prior_bel,sd_change),0,1)
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = abs(curr_effect - mean_effect)
        if agents[v].attitude == "ra"
            side = -1
        elseif agents[v].attitude == "rt"
            side = 1
        end
        dialpha = truncated(Normal(steps[dist]*side,sd_alpha),-1,1)
        the_change = the_diff*rand(dialpha)
        the_sign = rand(di) |> rand_sign
        if rand() <= 1/dist
            the_final = curr_effect + (the_change * the_sign)
        end
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    if the_final != agents[v].coop_effect
        agents[v].new_coop_effect = the_final
    else
        agents[v].new_coop_effect = agents[v].coop_effect
    end
end

##=================####==============##

function update_effect_given_distance_coop_alpha!(agents,g,d,threshold,step,sd_change,sd_alpha)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_alpha!(agents,g,d,threshold,step,sd_change,sd_alpha;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

#  function risk_assessment!(agents, g, d; v)
#      neigh = neighborhood_dists(g,v,d)
#      nodes = last.(neigh)[2:end]
#      distances = last.(neigh)[2:end]
#      curr_effect = agents[v].coop_effect
#      the_final = 0
#      for dist in unique(distances)
#          current = findall(x->x==dist, distances)
#          current_nodes = nodes[current]
#          populations = get_populations(agents[current_nodes]);
#          I = get(populations,"I",0)
#          if I > 0
#              the_final = curr_effect + (1/dist)*(I/length(current))*0.01#(length(current)/nv(g))
#          else
#              the_final = curr_effect - (1/dist)*(length(current)/length(nodes))*0.01#(length(current)/nv(g))
#          end
#      end
#      if the_final > 1
#          the_final = 1
#      elseif the_final < 0
#          the_final = 0
#      end
#      if the_final != agents[v].coop_effect
#          agents[v].new_coop_effect = the_final
#      else
#          agents[v].new_coop_effect = agents[v].coop_effect
#      end
#  end

##=================####==============##

function update_risk_assessment!(agents,g,d)
    Threads.@threads for ag in agents#[the_adaps]
        risk_assessment(agents, g, d; v=ag.id)
    end
    Threads.@threads for ag in agents#[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

function update_past!(agents,g,d,threshold,step,sd_change,sd_alpha)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_alpha!(agents,g,d,threshold,step,sd_change,sd_alpha;v=ag.id)
        end
    end
end

##=================####==============##

function risk_up!(agents, g, d; v)
    present = agents[v].coop_effect
    future = agents[v].new_coop_effect
    neigh = neighborhood_dists(g,v,d)
    nodes = last.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    the_prob = 0
    for dist in unique(distances)
        current = findall(x->x==dist, distances)
        current_nodes = nodes[current]
        populations = get_populations(agents[current_nodes]);
        I = get(populations,"I",0)
        if I > 0
            the_prob = (the_prob + (1/dist)*(I/length(current)))
        end
    end
    the_diff = abs(present - future)
    if rand() < the_prob
        updated = agents[v].coop_effect + the_diff
    if updated > 1
        updated = 1
    elseif updated < 0
        updated = 0
    end
    agents[v].new_coop_effect = updated
    end
end

##=================####==============##

function update_single_neighbors!(agents,g,d,threshold,steps,sd_change,sd_alpha;v)
    neigh = neighborhood_dists(g,v,d)
    nodes = first.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].coop_effect
    curr_bel = agents[v].prior_bel
    the_final = 0
    final_bel = 0
    di = truncated(Normal(agents[v].prior_bel,sd_change),0,1)
    for dist in unique(distances)
        current = findall(x->x==dist,distances)
        current_nodes = nodes[current]
        mean_effect = mean([x.coop_effect for x in agents[current_nodes]])
        the_diff = abs(curr_effect - mean_effect)
        populations = get_populations(agents[current_nodes]);
        I = get(populations,"I",0)
        if agents[v].attitude == "ra"
            side = -1
        elseif agents[v].attitude == "rt"
            side = 1
        end
        if I > 0
            final_bel = curr_bel + (steps[dist])*(I/length(current))
        else
            final_bel = curr_bel - (steps[dist])*(1/length(current))
        end
        dialpha = truncated(Normal(steps[dist]*side,sd_change),-1,1)
        the_change = the_diff*rand(dialpha)
        the_sign = rand(di) |> rand_sign
        if rand() <= 1/dist
            the_final = curr_effect + (the_change * the_sign)
        end
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    if final_bel > 1
        final_bel = 1
    elseif final_bel < 0
        final_bel = 0
    end
    agents[v].new_coop_effect = the_final
    agents[v].prior_bel = final_bel
end

##=================####==============##

function update_all_neighbors!(agents,g,d,threshold,step,sd_change,sd_alpha)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            #  update_single_effect_distance_alpha!(agents,g,d,threshold,step,sd_change,sd_alpha;v=ag.id)
            update_single_neighbors!(agents,g,d,threshold,step,sd_change,sd_alpha;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end

##=================####==============##

function risk_assessment!(agents, g, d, risk, sd_risk; v)
    neigh = neighborhood_dists(g,v,d)
    nodes = last.(neigh)[2:end]
    distances = last.(neigh)[2:end]
    curr_effect = agents[v].new_coop_effect
    the_final = curr_effect
    for dist in unique(distances)
        current = findall(x->x==dist, distances)
        current_nodes = nodes[current]
        populations = get_populations(agents[current_nodes]);
        I = get(populations,"I",0)
        num = risk*(1/dist)*(I/length(current))
        di = truncated(Normal(num,sd_risk),-1,1)
        the_final = the_final + rand(di)
    end
    if the_final > 1
        the_final = 1
    elseif the_final < 0
        the_final = 0
    end
    if the_final != agents[v].new_coop_effect
        agents[v].new_coop_effect = the_final
    else
        agents[v].new_coop_effect = agents[v].new_coop_effect
    end
end

##=================####==============##

function update_effect_given_distance_coop_alpha_risk!(agents,g,d,threshold,step,sd_change,sd_alpha,risk, sd_risk)
    the_adaps = findall(x->x.adapter,agents)
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            update_single_effect_distance_alpha!(agents,g,d,threshold,step,sd_change,sd_alpha;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        if ag.adapter == true
            risk_assessment!(agents,g,d, risk, sd_risk;v=ag.id)
        end
    end
    Threads.@threads for ag in agents[the_adaps]
        ag.coop_effect = ag.new_coop_effect
    end
end


text:
# module IntegrationTests # TODO: put in module?

_to_cleanup = Array{Any,1}()
schedule_for_cleanup(item) = push!(_to_cleanup,item)

function cleanup(items)
	for i=length(items):-1:1 # iterate in reverse order
		item = items[i]

		if typeof(item) <: AbstractEntity || utils.is_synapse_id(item) != nothing ||
		   (typeof(item) <: AbstractSynapse && hasattr(item,"deleteURI"))
			try
				delete(syn, item)
			catch ex
				# TODO: rewrite if exception handling is improved
				if ex <: PyError && hasattr(ex.val, "response") && ex.val[:response][:status_code] in [404,403]
				else
					println("Error cleaning up entity: ", ex)
				end
			end
		elseif typeof(item) <: AbstractString
			if ispath(item)
				try
					rm(item, recursive=true) # works on both files and folders
				catch ex
					println(ex)
				end
			end
		else
			println(stderr, "Don't know how to clean: $item")
		end
	end
end

syn = SynapseClient.Synapse(debug=true, skip_checks=true)

print("Testing against endpoints:")
print("  " * syn.repoEndpoint)
print("  " * syn.authEndpoint)
print("  " * syn.fileHandleEndpoint)
print("  " * syn.portalEndpoint * "\n")

login(syn)

project = Project() # for scope reasons

try
    global project = store(syn, Project(name=string(UUIDs.uuid4())))
    schedule_for_cleanup(project)

    println("Created project:")
    println(project.id)

    include("integration_test.jl")

#catch ex
finally
	cleanup(_to_cleanup)
end

# end


text:
function EstimatedIonosphericModel(arg0::IonosphericMappingFunction, arg1::jdouble)
    return EstimatedIonosphericModel((IonosphericMappingFunction, jdouble), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_parameters(obj::IonosphericModel)
    return jcall(obj, "getParameters", Vector{jdouble}, ())
end

function get_parameters(obj::IonosphericModel, arg0::Field)
    return jcall(obj, "getParameters", Vector{RealFieldElement}, (Field,), arg0)
end

function get_parameters_drivers(obj::EstimatedIonosphericModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function path_delay(obj::EstimatedIonosphericModel, arg0::FieldSpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldSpacecraftState, TopocentricFrame, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::EstimatedIonosphericModel, arg0::RealFieldElement, arg1::jdouble, arg2::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (RealFieldElement, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function path_delay(obj::EstimatedIonosphericModel, arg0::SpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (SpacecraftState, TopocentricFrame, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::EstimatedIonosphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (jdouble, jdouble, Vector{jdouble}), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# Velocity's direction is alway defined relative to the +X axis.
# Default direction is +X axis.

mutable struct Velocity{T <: AbstractFloat}
    magnitude::T
    min_magnitude::T
    max_magnitude::T

    direction::Vector2D{T}

    limit_magnitude::Bool

    function Velocity{T}() where {T <: AbstractFloat}
        o = new()
        o.magnitude = 0.0
        o.min_magnitude = 0.0
        o.max_magnitude = 0.0
        o.direction = Vector2D{T}()
        o.limit_magnitude = true
        o
    end

    function Velocity{T}(velocity::Velocity{T}) where {T <: AbstractFloat}
        o = new()
        o.magnitude = velocity.magnitude
        o.min_magnitude = velocity.min_magnitude
        o.max_magnitude = velocity.max_magnitude
        Math.set!(o.direction, velocity.direction)
        o
    end

    function Velocity{T}(x::T, y::T, speed::T) where {T <: AbstractFloat}
        o = new()
        o.magnitude = speed
        Math.set!(o.direction, x, y)
        o
    end

end

function set!(v::Velocity{T}, from::Velocity{T}) where {T <: AbstractFloat}
    v.magnitude = from.magnitude
    v.min_magnitude = from.min_magnitude
    v.max_magnitude = from.max_magnitude
    Math.set!(v.direction, from.direction);
end

function set_magnitude_range!(v::Velocity{T}, min::T, max::T) where {T <: AbstractFloat}
    v.min_magnitude = min
    v.max_magnitude = max;
end

function set_magnitude!(v::Velocity{T}, magnitude::T) where {T <: AbstractFloat}
    if v.limit_magnitude
        # Push magnitude into the range
        v.magnitude = min(magnitude, v.max_magnitude)
        v.magnitude = max(v.magnitude, v.min_magnitude)
    else
        v.magnitude = magnitude
    end
end

function get_magnitude(v::Velocity{T}) where {T <: AbstractFloat}
    v.magnitude
end

function accelerate!(v::Velocity{T}, acceleration::T) where {T <: AbstractFloat}
    v.magnitude += acceleration
    set_magnitude!(v, v.magnitude);
end

function increase_speed!(v::Velocity{T}, acceleration::T) where {T <: AbstractFloat}
    v.magnitude += acceleration
    if v.limit_magnitude
        v.magnitude = min(v.magnitude, v.max_magnitude)
    end
end

function decrease_speed!(v::Velocity{T}, acceleration::T) where {T <: AbstractFloat}
    v.magnitude -= acceleration
    if v.limit_magnitude
        v.magnitude = max(v.magnitude, v.min_magnitude)
    end
end

function set_velocity!(v::Velocity{T}, degrees::T, speed::T) where {T <: AbstractFloat}
    set_direction!(v.direction, degrees)
    set_magnitude!(v, speed);
end

function set_direction!(v::Velocity{T}, degrees::T) where {T <: AbstractFloat}
    set_direction!(v.direction, degrees);
end

function set_direction!(v::Velocity{T}, direction::Vector2D{T}) where {T <: AbstractFloat}
    set_direction!(v.direction, direction);
end

function get_velocity!(v::Velocity{T}, velocity::Vector2D{T}) where {T <: AbstractFloat}
    set!(velocity, v.direction.x * v.magnitude, v.direction.y * v.magnitude);
end

# Relative to +X axis in radians
function get_angle(v::Velocity{T}) where {T <: AbstractFloat}
    atan(v.direction.x, v.direction.y)
end

# vA = vA + vB
function add!(vA::Velocity{T}, vB::Velocity{T}) where {T <: AbstractFloat}
    # Convert to vectors first
    get_velocity!(vA, v1)
    get_velocity!(vB, v2)

    add!(v1, v2, v3)
    len = length(v3)

    # limit/clamp
    vA.magnitude = min(len, vA.max_magnitude)

    set!(vA.direction, v3)
    Math.normalize!(vA.direction)
end

# vA = vB - vA
function sub!(vA::Velocity{T}, vB::Velocity{T}) where {T <: AbstractFloat}
    # Convert to vectors first
    get_velocity!(vA, v1)
    get_velocity!(vB, v2)

    sub!(v2, v1, v3)
    len = length(v3)

    # limit/clamp
    vA.magnitude = min(len, vA.max_magnitude)

    set!(vA.direction, v3)
    Math.normalize!(vA.direction)
end

# Apply velocity to vector
function apply!(velocity::Velocity{T}, vec::Vector2D{T}) where {T <: AbstractFloat}
    get_velocity!(velocity, v1)
    add!(vec, v1, vec);
end

function apply!(velocity::Velocity{T}, point_mass::Point{T}) where {T <: AbstractFloat}
    get_velocity!(velocity, v1)

    set!(v2, point_mass.x, point_mass.y)
    add!(v2, v1, v3);

    Geometry.set!(point_mass, v3.x, v3.y)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
export Circulant
"""
[`Circulant` matrix](https://en.wikipedia.org/wiki/Circulant_matrix)

```julia
julia> Circulant([1,2,3,4])
4x4 Circulant{Int64}:
 1  4  3  2
 2  1  4  3
 3  2  1  4
 4  3  2  1
```
"""
struct Circulant{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

getindex(C::Circulant, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
isassigned(C::Circulant, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
size(C::Circulant, r::Int) = (r==1 || r==2) ? length(C.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(C::Circulant) = size(C,1), size(C,2)

# Fast matrix x vector via fft()
# see Golub, van Loan, Matrix Computations, John Hopkins, Baltimore, 1996, p. 202
function *(C::Circulant{T},x::Vector{T}) where T
    xt=fft(x)
    vt=fft(C.c)
    yt=vt.*xt
    typeof(x[1])==Int ? map(Int,round.(real(ifft(yt)))) : ( (T <: Real) ? map(T,real(ifft(yt))) : ifft(yt))
end

function mul!(y::StridedVector{T},C::Circulant{T},x::StridedVector{T}) where T
    xt=fft(x)
    vt=fft(C.c)
    yt=ifft(vt.*xt)
    if T<: Int
        map!(round,y,yt)
    elseif T<: Real
        map!(real,y,yt)
    else
        copy!(y,yt)
    end
    return y
end

function Matrix(C::Circulant{T}) where T
	n=size(C, 1)
	M=Array{T}(undef,n,n)
	for i=1:n
		M[i:n,i] = C.c[1:n-i+1]
		M[1:i-1,i] = C.c[n-i+2:n]
	end
	M
end


text:
"""
    module ILog2

This module provides `ilog2`.
"""
module ILog2

export ilog2, checkispow2

const IntBits  = Union{Int8, Int16, Int32, Int64, Int128,
                       UInt8, UInt16, UInt32, UInt64, UInt128}

"""
    ilog2(x, RoundUp)

Return the smallest `m` such that `2^m >= n`.
"""
ilog2(x, ::typeof(RoundUp)) = ispow2(x) ? ilog2(x) : ilog2(x) + 1
ilog2(x, ::typeof(RoundDown)) = ilog2(x)

"""
    msbindex(::Type{T})
    T is an Integer bits type

Return one less than the number of bits in T
"""
@generated function msbindex(::Type{T}) where {T<:Integer}
    return sizeof(T)*8 - 1
end

msbindex(::Type{BigInt}) = throw(ArgumentError("Expected an integer type with fixed number of bits"))

"""
    ilog2(n::Real)

Compute the largest `m` such that `2^m <= n`.

!!! note

    `float(2^53 - 1)` is the largest `n::Float64` for which `eps(n) ≤ 1`.
    For larger `n::Float64`, "the largest `m` such that `2^m <= n`" is ambiguous, although
    `ilog2` may return a number. For large enough `n::Float64`, `ilog2` will
    throw an `InexactError`. These cautionary statements do not apply for `n::Integer`.
"""
function ilog2(n::T) where {T<:IntBits}
    n > zero(T) && return msbindex(T) - leading_zeros(n)
    throw(DomainError(n))
end
ilog2(n::BigInt) = Base.GMP.MPZ.sizeinbase(n, 2) - 1
ilog2(x::Real) = ilog2(convert(Integer, floor(x)))

# This is several times slower than the other methods. But none of the standard bitstype integers,
# nor `BigInt`, dispatch to this method.
ilog2(n::Integer) = convert(typeof(n), floor(log(2,n)))

"""
    checkispow2(n::Number)

Return base-2 logarithm of `n` if `n` is a power of two.
Otherwise throw a `DomainError`.
"""
function checkispow2(n::Number)
    if ! ispow2(n)
        throw(DomainError(n, "$n is not a power of two."))
    end
    return ilog2(n)
end

end # module ILog2


text:
"""
---------------------------------------------------
Types and Constructors
---------------------------------------------------
"""

mutable struct VERTICAL_ASCENT_SAMPLER <: SAMPLER
	z_final::Sampleable
	z_ascent::Sampleable
	vx_final::Sampleable
	θ::Sampleable
	t_extra::Sampleable
	t_ascent::Sampleable
end

function vertical_ascent_sampler(;z_final = Normal(500,20),
						z_ascent = Normal(100,8),
						vx_final = Uniform(30,45),
						θ = Uniform(3,5),
						t_extra = Uniform(0.1,0.3), #Normal(10,4),
						t_ascent = Normal(30,5))
	return VERTICAL_ASCENT_SAMPLER(z_final, z_ascent, vx_final, θ, t_extra, t_ascent)
end

mutable struct VERTICAL_ASCENT <: UAM_TRAJECTORY
	z_final::Float64
	z_ascent::Float64
	vx_final::Float64
	θ::Float64
	t_extra::Float64
	t_ascent::Float64
	sampler::SAMPLER
	dt::Float64
	vz_max::Float64
	az_max::Float64
	ax_max::Float64
	p::Array{Float64,2}
	v::Array{Float64,2}
end

function vertical_ascent(;z_final = 0.0,
						z_ascent = 0.0,
						vx_final = 0.0,
						θ = 0.0,
						t_extra = 0.0,
						t_ascent = 0.0,
						sampler = vertical_ascent_sampler(),
						dt = 1.0,
						vz_max = 550fpm2fps,
						az_max = 0.3g_ft,
						ax_max = 0.1g,
						p = zeros(1,3),
						v = zeros(1,3))
	return VERTICAL_ASCENT(z_final, z_ascent, vx_final, θ, t_extra, t_ascent, sampler, dt,
							vz_max, az_max, ax_max, p, v)
end

"""
---------------------------------------------------
Functions
---------------------------------------------------
"""

function solve_trajectory!(τ::VERTICAL_ASCENT)
	# Vertical ascent portion ######################################
	Nd = convert(Int64, floor(τ.t_ascent/τ.dt))
	# Difference matrix
	Dd = zeros(Nd-1,Nd)
	for i = 1:(Nd-1)
	    Dd[i,i] = -1
	    Dd[i,i+1] = 1
	end
	Dd = Dd./τ.dt

	pzd = Variable(Nd)
	vzd = Variable(Nd)

	objec = norm(Dd*vzd)

	constraints = [[pzd[i+1] == pzd[i] + vzd[i]*τ.dt for i = 1:Nd-1];
				   vzd[1] == 0; pzd[1] == 0;
				   pzd[Nd] == τ.z_ascent; vzd[Nd] == 0;
				   pzd ≥ 0; vzd ≥ 0; vzd ≤ τ.vz_max;
				   Dd*vzd*ft2m ≥ -τ.az_max; Dd*vzd*ft2m ≤ τ.az_max]

	prob2 = minimize(objec, constraints)
	solve!(prob2, ECOS.Optimizer(verbose=0))

	# Glide portion #################################################
	# Get some needed info
	#ttot = abs(τ.z_final/500fpm2fps) + τ.t_extra*abs(τ.z_final/500fpm2fps)
	ttot = abs(τ.z_final/τ.vz_max) + τ.t_extra*abs(τ.z_final/τ.vz_max)
	N = convert(Int64, floor(ttot/τ.dt))
	# Difference matrix
	D = zeros(N-1,N)
	for i = 1:(N-1)
	    D[i,i] = -1
	    D[i,i+1] = 1
	end
	D = D./τ.dt

	# Create and solve the optimization problem
	px = Variable(N)
	pz = Variable(N)
	vx = Variable(N)
	vz = Variable(N)

	objec = norm(D*vx) + λ*norm(vz*ft2m + vx*tand(τ.θ))

	constraints = [[px[i+1] == px[i] + vx[i]*τ.dt for i = 1:N-1];
				   [pz[i+1] == pz[i] + vz[i]*τ.dt for i = 1:N-1];
				   vx[1] == 0; px[1] == 0;
				   vz[1] == 0; pz[1] == τ.z_ascent;
				   vx[N] == τ.vx_final; pz[N] == τ.z_final; vz[N] == 0;
				   vx ≥ 0; pz ≥ 0; vz ≥ 0; vz ≤ τ.vz_max;
				   D*vx ≤ τ.ax_max; D*vz*ft2m ≤ τ.az_max;
				   D*vx ≥ -τ.ax_max; D*vz*ft2m ≥ -τ.az_max]

	prob1 = minimize(objec, constraints)
	solve!(prob1, ECOS.Optimizer(verbose=0))

	# Update the trajectory
	p = zeros(Nd+N,3)
	p[Nd+1:end,1] = px.value[:,1]
	p[1:Nd,3] = pzd.value[:,1]
	p[Nd+1:end,3] = pz.value
	v = zeros(Nd+N,3)
	v[Nd+1:end,1] = vx.value[:,1]
	v[1:Nd,3] = vzd.value[:,1]
	v[Nd+1:end,3] = vz.value

	τ.p = copy(p)
	τ.v = copy(v)

	return Int(prob1.status) == 1
end

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


text:
using HyperbolicPlane, Plots

"""
`regular(n,r)` draw a regular `n`-gon in the hyperbolic plane.
The vertices are distance `r` from the origin.
"""
function regular(n::Int=5, r::Real=1.0)
    plot()
    pts = [ HPoint(r, (k/n)*2pi) for k=0:n ]
    for P in pts
        set_radius(P,3)
    end
    draw(pts)

    segs = [ HSegment(pts[k],pts[k+1]) for k=1:n ]
    draw(segs)
    draw(HPlane())
    finish()
end


text:
# VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module Script

export _nullFunction
export _debug
export compile
export invoke

global const _nullFunction = function() end
global _debug = true

function compile(file::String)
	result = nothing
	try
		result = evalfile(file)
		if _debug == true
			println("file compiled.")
		end
	catch err
		println("Cannot compile file. There was an error: ", err)
	end
	return result
end

function invoke(f::Function, args...)
	result = nothing
	try
		if f == _nullFunction
			error("function is null")
		end
		result = f(args...)
		if _debug == true
			println("function invoked. result: ", result)
		end
	catch err
		println("Cannot invoke function. There was an error: ", err)
	end
	return result
end

end # module


text:
# This example is primarily for expository purposes. Do not copy
# the practices demosntrated here!!! Instead refer to some of the
# other examples for best practices.

using Weber
sid,skip = @read_args("A simple frequency discrimination experiment.")

const low = ramp(tone(1kHz,0.5s))
const high = ramp(tone(1.1kHz,0.5s))

function one_trial()
  if rand(Bool)
    stim1 = moment(0.5s,play,low)
	  stim2 = moment(0.5s,play,high)
    resp = response(key"q" => "low_first", key"p" => "low_second",correct = "low_first")
  else
	  stim1 = moment(0.5s,play,high)
	  stim2 = moment(0.5s,play,low)
    resp = response(key"q" => "low_first", key"p" => "low_second",correct = "low_second")
  end
  return [show_cross(),stim1,stim2,resp,await_response(iskeydown)]
end

exp = Experiment(columns = [:sid => sid,:condition => "ConditionA",:correct],skip=skip)
setup(exp) do
  addbreak(instruct("Press 'q' when you hear the low tone first and 'p' otherwise."))
  for trial in 1:10
    addtrial(one_trial())
  end
end

run(exp)


text:
guessValidationFunction(v::Number) = isNumber
guessValidationFunction(v::Any) = (v) -> true

validationPassedIcon = ""
validationFailedIcon = "edit-delete-symbolic"

function ValidationEntry(value::Any, validation::Function=guessValidationFunction(value), placeholder::String=string(value); setPlaceholder::Bool=true)
    text = string(value)
    entry = Entry(text=text)

    if setPlaceholder
        GAccessor.placeholder_text(entry, placeholder)
    end

    icon = validation(text) ? validationPassedIcon : validationFailedIcon
    set_gtk_property!(entry, :secondary_icon_name, icon)

    @guarded signal_connect(entry, "changed") do widget
        text = Gtk.bytestring(GAccessor.text(widget))
        icon = validation(text) ? validationPassedIcon : validationFailedIcon

        set_gtk_property!(widget, :secondary_icon_name, icon)
    end

    return entry
end


text:
"""
    unsafe_append_buffer!(schunk::SChunk, src::Ptr{T}, src_size)

Compress `src_size` bytes from `src` buffer, and append it to super-chubk

The `unsafe` prefix on this function indicates that no validation is performed on the pointer `src`
to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.

Return the number of chunks in super-chunk. If some problem is detected, this number will be negative.
"""
function unsafe_append_buffer!(schunk::SChunk, src::Ptr{T}, src_size) where {T}
    return Lib.blosc2_schunk_append_buffer(schunk.ptr, src, src_size)
end

"""
    unsafe_decompress_chunk(schunk::SChunk, nchunk, dest::Ptr{T}, dest_size)

Decompress the `nchunk` (1-indexed) chunk from super-chunk into `dest`

Return count of the elements written to dest

The `unsafe` prefix on this function indicates that no validation is performed on the pointer `src`
to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_decompress_chunk(schunk::SChunk, nchunk, dest::Ptr{T}, dest_size) where {T}
    @boundscheck checkchunksbounds(schunk, nchunk)
    res = Lib.blosc2_schunk_decompress_chunk(schunk.ptr, nchunk - 1, dest, dest_size)
    return res >= 0 ?
    ceil(Int64, res / sizeof(T)) :
    res
end

"""
    unsafe_append_chunk!(schunk::SChunk, chunk::Ptr{UInt8}; copy::Bool = true)

Append an existing `chunk` o a super-chunk.
`copy` specified whether the chunk should be copied internally or can be used as-is.
Return the number of chunks in super-chunk. If some problem is detected, this number will be negative.

The `unsafe` prefix on this function indicates that no validation is performed on the pointer `chunk`
to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_append_chunk!(schunk::SChunk, chunk::Ptr{UInt8}; copy::Bool = true)
    return Lib.blosc2_schunk_append_chunk(schunk.ptr, chunk, copy)
end

"""
    unsafe_update_chunk!(schunk::SChunk, nchunk, chunk::Ptr{UInt8}; copy::Bool = true)

Update a `chunk` at a `nchunk` position (1-indexed) in a super-chunk.
`copy` specified whether the chunk should be copied internally or can be used as-is.
Return the number of chunks in super-chunk. If some problem is detected, this number will be negative.

The `unsafe` prefix on this function indicates that no validation is performed on the pointer `chunk`
to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_update_chunk!(schunk::SChunk, nchunk, chunk::Ptr{UInt8}; copy::Bool = true)
    @boundscheck checkchunksbounds(schunk, nchunk)
    return Lib.blosc2_schunk_update_chunk(schunk.ptr, nchunk - 1, chunk, copy)
end

"""
    unsafe_insert_chunk!(schunk::SChunk, nchunk, chunk::Ptr{UInt8}; copy::Bool = true)

Insert a `chunk` at a `nchunk` position (1-indexed) in a super-chunk.
`copy` specified whether the chunk should be copied internally or can be used as-is.
Return the number of chunks in super-chunk. If some problem is detected, this number will be negative.

The `unsafe` prefix on this function indicates that no validation is performed on the pointer `chunk`
to ensure that they are valid. Incorrect usage may corrupt or segfault your program, in the same manner as C.
"""
function unsafe_insert_chunk!(schunk::SChunk, nchunk, chunk::Ptr{UInt8}; copy::Bool = true)
    @boundscheck checkchunksbounds(schunk, nchunk)
    return Lib.blosc2_schunk_insert_chunk(schunk.ptr, nchunk - 1, chunk, copy)
end

"""
    unsafe_get_chunk(schunk::SChunk, nchunk)

Return a compressed chunk that is part of a super-chunk as `Vector{UInt8}`

The unsafe prefix on this function indicates that using the result of this function
after the schunk argument to this function is no longer accessible to the program may cause undefined behavior,
including program corruption or segfaults, at any later time.
"""
function unsafe_get_chunk(schunk::SChunk, nchunk)
    @boundscheck checkchunksbounds(schunk, nchunk)
    result = Ref{Ptr{UInt8}}()
    needs_free = Ref{Bool}()
    size = Lib.blosc2_schunk_get_chunk(schunk.ptr, nchunk - 1, result, needs_free)
    size < 0 && error("decompression error")
    return unsafe_wrap(Array, result[], size, own = needs_free[])
end


text:
mutable struct Measure
    mean::Float64
    error::Float64
end

Measure() = Measure(0,-1)
Measure(a::Observable) = Measure(mean(a), error(a))
Measure(a::Measure) = Measure(a.mean, a.error)

function measure_binomial(nsucc::Integer, trials::Integer)
    p = nsucc /trials
    e = sqrt(p*(1-p) / trials)
    return Measure(p, e)
end

function (&)(a::Measure, b::Measure)
	if a.error <= 0
        a.mean = b.mean
        a.error = b.error
        return
    end
    if b.error <= 0
        return
    end

    w1 = 1 / a.error^2
    w2 = 1 / b.error^2
    m = (a.mean * w1 + b.mean * w2 ) / (w1+w2)
    e = 1 / sqrt(w1+w2)
    return Measure(m, e)
end

function +(a::Measure, b::Measure)
	m = a.mean + b.mean
    e = sqrt(a.error^2 + b.error^2)
    Measure(m,e)
end

function *(a::Measure, b::Measure)
	m = a.mean*b.mean
    e = sqrt(b.mean^2*a.error^2 + a.mean^2*b.error^2)
    Measure(m,e)
end

+(a::Measure, val::Number) = a + Measure(val,0)

function *(a::Measure, val::Number)
	m = a.mean * val
    e = a.error * val
    Measure(m,e)
end

function shortshow(io::IO, a::Measure)
    print(io,a.mean, " ", a.error)
end

Base.show(io::IO, a::Measure) = shortshow(io, a)

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
####
#### Defines list of tutorials given `GENERATED_DIR`
####

generate_tutorials =
    parse(Bool, get(ENV, "CLIMATEMACHINE_DOCS_GENERATE_TUTORIALS", "true"))

tutorials = []

# Allow flag to skip generated
# tutorials since this is by
# far the slowest part of the
# docs build.
if generate_tutorials

    # generate tutorials

    include("pages_helper.jl")

    tutorials = [
        "Home" => "TutorialList.jl",
        "Balance Law" => "how_to_make_a_balance_law.jl",
        "Atmos" => [
            "Dry Idealized GCM (Held-Suarez)" => "Atmos/heldsuarez.jl",
            "Single Element Stack Experiment (Burgers Equation)" =>
                "Atmos/burgers_single_stack.jl",
            "LES Experiment (Density Current)" => "Atmos/densitycurrent.jl",
            "LES Experiment (Rising Thermal Bubble)" =>
                "Atmos/risingbubble.jl",
            "Linear Hydrostatic Mountain (Topography)" =>
                "Atmos/agnesi_hs_lin.jl",
            "Linear Non-Hydrostatic Mountain (Topography)" =>
                "Atmos/agnesi_nh_lin.jl",
        ],
        "Land" => [
            "Heat" => ["Heat Equation" => "Land/Heat/heat_equation.jl"],
            "Soil" => [
                "Hydraulic Functions" =>
                    "Land/Soil/Water/hydraulic_functions.jl",
                "Soil Heat Equation" =>
                    "Land/Soil/Heat/bonan_heat_tutorial.jl",
                "Richards Equation" =>
                    "Land/Soil/Water/equilibrium_test.jl",
                "Coupled Water and Heat" =>
                    "Land/Soil/Coupled/equilibrium_test.jl",
            ],
        ],
        "Ocean" => [
            "One-dimensional geostrophic adjustment" =>
                "Ocean/geostrophic_adjustment.jl",
            "Propagating mode-1 internal wave" => "Ocean/internal_wave.jl",
            "Shear instability" => "Ocean/shear_instability.jl",
        ],
        "Numerics" => [
            "System Solvers" => [
                "Conjugate Gradient" => "Numerics/SystemSolvers/cg.jl",
                "Batched Generalized Minimal Residual" =>
                    "Numerics/SystemSolvers/bgmres.jl",
            ],
            "DG Methods" =>
                ["Filters" => "Numerics/DGMethods/showcase_filters.jl"],
        ],
        "Diagnostics" => [
            "Debug" => [
                "State Statistics Regression" =>
                    "Diagnostics/Debug/StateCheck.jl",
            ],
        ],
    ]

    # Prepend tutorials_dir
    tutorials_jl = flatten_to_array_of_strings(get_second(tutorials))
    println("Building literate tutorials...")

    @everywhere function generate_tutorial(tutorials_dir, tutorial)
        gen_dir =
            joinpath(GENERATED_DIR, relpath(dirname(tutorial), tutorials_dir))
        input = abspath(tutorial)
        script = Literate.script(input, gen_dir)
        code = strip(read(script, String))
        mdpost(str) = replace(str, "@__CODE__" => code)
        Literate.markdown(input, gen_dir, postprocess = mdpost)
        Literate.notebook(input, gen_dir, execute = true)
    end

    tutorials_dir = joinpath(@__DIR__, "..", "tutorials")
    tutorials_jl = map(x -> joinpath(tutorials_dir, x), tutorials_jl)
    pmap(t -> generate_tutorial(tutorials_dir, t), tutorials_jl)

    # update list of rendered markdown tutorial output for mkdocs
    ext_jl2md(x) = joinpath(basename(GENERATED_DIR), replace(x, ".jl" => ".md"))
    tutorials = transform_second(x -> ext_jl2md(x), tutorials)
end

max_forks_repo_name: jm-c/CLIMA


text:
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
using Optim
f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)   # random starting value
result = optimize(minusf, startval, BFGS())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
using DataFrames
using CSV
using HTTP
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)

using GLM
bols = inv(X'*X)*X'*y
df.white = df.race.==1
bols_lm = lm(@formula(married ~ age + white + collgrad), df)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function logit(alpha, X, d)

    # your turn

    return loglike
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# see Lecture 3 slides for example

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==10,:occupation] .= 9
df[df.occupation.==11,:occupation] .= 9
df[df.occupation.==12,:occupation] .= 9
df[df.occupation.==13,:occupation] .= 9
freqtable(df, :occupation) # problem solved

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mlogit(alpha, X, d)

    # your turn

    return loglike
end


text:
@generated function Base.promote_rule{M1, N1, R1, T1, M2, N2, R2, T2}(::Type{SIMDVector{M1, N1, R1, T1}},
                                                                      ::Type{SIMDVector{M2, N2, R2, T2}})
    Tr = promote_type(T1, T2)

    simd_len, rest, buckets = compute_lengths(M1*N1 + R1, Tr)
    return quote
        $(Expr(:meta, :inline))
        SIMDVector{$buckets, $simd_len, $rest, $Tr}
    end
end

function generate_conversion_expr(simd_len, N, rest, T)
    simd_array_create_expr = Expr(:tuple)
    if simd_len != 0
        for i in 1:simd_len:N-rest
            exp_simd_ele = Expr(:call, :VecRegister)
            push!(exp_simd_ele.args, Expr(:tuple, [:(VecElement(convert($T, v[$j]))) for j in i:simd_len+i-1]...))
            push!(simd_array_create_expr.args, exp_simd_ele)
        end
    end

    rest_array_create_expr = Expr(:tuple, [:(convert($T, v[$j])) for j in (N-rest+1):N]...)

    return simd_array_create_expr, rest_array_create_expr
end

# Identity conversion
@inline Base.convert{M, N, R, T}(::Type{SIMDVector{M, N, R, T}}, vec::SIMDVector{M, N, R, T}) = vec

@generated function Base.convert{M1, N1, R1, T1, M2, N2, R2, T2}(::Type{SIMDVector{M1, N1, R1, T1}},
                                                                 v::SIMDVector{M2, N2, R2, T2})
    Tr = promote_type(T1, T2)
    simd_len, rest, buckets = compute_lengths(M1*N1 + R1, Tr)

    simd_array_create_expr, rest_array_create_expr = generate_conversion_expr(simd_len, M1*N1 + R1, rest, Tr)
    return quote
        $(Expr(:meta, :inline))
        @assert M1 * N1 + R1 == M2 * N2 + R2
        @inbounds simd_tup = $simd_array_create_expr
        @inbounds rest = $rest_array_create_expr
        SIMDVector{$buckets, $simd_len, $rest, $Tr}(simd_tup, rest)
    end
end

promote_eltype{M, N, R, T <: Number}(v::SIMDVector{M, N, R, T}, n::T) = v, n
promote_eltype{M, N, R, T <: Number}(n::T, v::SIMDVector{M, N, R, T}) = n, v

@generated function promote_eltype{M, N, R, T1 <: Number, T2 <: Number}(v::SIMDVector{M, N, R, T1},
                                                                        n::T2)
    Tr = promote_type(T1, T2)
    if Tr == T1
        return quote
            return v, convert(T1, n)
        end
    end

    simd_len, rest, buckets = compute_lengths(M*N + R, Tr)
    simd_array_create_expr, rest_array_create_expr = generate_conversion_expr(simd_len, M*N + R, rest, Tr)
    return quote
        $(Expr(:meta, :inline))
        @inbounds simd_tup = $simd_array_create_expr
        @inbounds rest = $rest_array_create_expr
        return SIMDVector{$buckets, $simd_len, $rest, $Tr}(simd_tup, rest), convert($Tr, n)
    end
end

function promote_eltype{M, N, R, T1 <: Number, T2 <: Number}(n::T1, v::SIMDVector{M, N, R, T2})
    vr, nr = promote_eltype(v, n)
    return nr, vr
end


text:
export
        LongitudinalDriverModel,
        ProportionalSpeedTracker,
        IntelligentDriverModel,
        StaticLongitudinalDriver,
        PrincetonLongitudinalDriver,
        track_longitudinal!

abstract LongitudinalDriverModel
get_name(::LongitudinalDriverModel) = "???"
set_desired_speed!(::LongitudinalDriverModel, v_des::Float64) = model # # do nothing by default
reset_hidden_state!(model::LongitudinalDriverModel) = model # do nothing by default
track_longitudinal!(model::LongitudinalDriverModel, scene::Scene, roadway::Roadway, ego_index::Int, target_index::Int) = model # do nothing by default
observe!(model::LongitudinalDriverModel, scene::Scene, roadway::Roadway, egoid::Int) = model  # do nothing by default
Base.rand(model::LongitudinalDriverModel) = error("rand not implemented for model $model")
Distributions.pdf(model::LongitudinalDriverModel, a_lon::Float64) = error("pdf not implemented for model $model")
Distributions.logpdf(model::LongitudinalDriverModel, a_lon::Float64) = error("logpdf not implemented for model $model")

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

type StaticLongitudinalDriver <: LongitudinalDriverModel
    a::Float64
    StaticLongitudinalDriver(a::Float64=0.0) = new(a)
end
get_name(::StaticLongitudinalDriver) = "ProportionalSpeedTracker"
Base.rand(model::StaticLongitudinalDriver) = model.a
Distributions.pdf(model::StaticLongitudinalDriver, a_lon::Float64) = a_lon == model.a ? Inf : 0.0
Distributions.logpdf(model::StaticLongitudinalDriver, a_lon::Float64) = a_lon == model.a ? Inf : -Inf

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

type ProportionalSpeedTracker <: LongitudinalDriverModel
    a::Float64 # predicted acceleration
    σ::Float64 # optional stdev on top of the model, set to zero or NaN for deterministic behavior
    k::Float64 # proportional constant for speed tracking [s⁻¹]
    v_des::Float64 # desired speed [m/s]

    function ProportionalSpeedTracker(;
        σ::Float64 = NaN,
        k::Float64 = 1.0,
        v_des::Float64 = 29.0,
        )

        retval = new()
        retval.a = NaN
        retval.σ = σ
        retval.k = k
        retval.v_des = v_des
        retval
    end
end
get_name(::ProportionalSpeedTracker) = "ProportionalSpeedTracker"
function set_desired_speed!(model::ProportionalSpeedTracker, v_des::Float64)
    model.v_des = v_des
    model
end
function track_longitudinal!(model::ProportionalSpeedTracker, scene::Scene, roadway::Roadway, ego_index::Int, target_index::Int)
    veh_ego = scene[ego_index]
    v = veh_ego.state.v

    Δv = model.v_des - v
    model.a = Δv*model.k # predicted accel to match target speed

    model
end
function observe!(model::ProportionalSpeedTracker, scene::Scene, roadway::Roadway, egoid::Int)
    ego_index = get_index_of_first_vehicle_with_id(scene, egoid)
    track_longitudinal!(model, scene, roadway, ego_index, 0)
    model
end
function Base.rand(model::ProportionalSpeedTracker)
    if isnan(model.σ) || model.σ ≤ 0.0
        model.a
    else
        rand(Normal(model.a, model.σ))
    end
end
function Distributions.pdf(model::ProportionalSpeedTracker, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a_lon)
    end
end
function Distributions.logpdf(model::ProportionalSpeedTracker, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a_lon)
    end
end

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

type PrincetonLongitudinalDriver <: LongitudinalDriverModel
    a::Float64 # predicted acceleration
    σ::Float64 # optional stdev on top of the model, set to zero or NaN for deterministic behavior
    k::Float64 # proportional constant for speed tracking [s⁻¹]
    v_des::Float64 # desired speed [m/s]

    function PrincetonLongitudinalDriver(;
        σ::Float64 = NaN,
        k::Float64 = 1.0,
        v_des::Float64 = 29.0,
        )

        retval = new()
        retval.a = NaN
        retval.σ = σ
        retval.k = k
        retval.v_des = v_des
        retval
    end
end
get_name(::PrincetonLongitudinalDriver) = "PrincetonLongitudinalDriver"
function set_desired_speed!(model::PrincetonLongitudinalDriver, v_des::Float64)
    model.v_des = v_des
    model
end
function track_longitudinal!(model::PrincetonLongitudinalDriver, scene::Scene, roadway::Roadway, ego_index::Int, target_index::Int)

    veh_ego = scene[ego_index]
    v_des = model.v_des

    if target_index > 0
        veh_target = scene[target_index]
        speed_M = veh_target.state.v
        dist_M = abs(veh_target.state.posG - veh_ego.state.posG) - veh_target.def.length/2 - veh_ego.def.length/2
        v_des = min(speed_M*(1-exp(-model.k*dist_M/speed_M - 1)), model.v_des)
    end

    Δv = v_des - veh_ego.state.v
    model.a = Δv*model.k # predicted accel to match target speed

    model
end
function observe!(model::PrincetonLongitudinalDriver, scene::Scene, roadway::Roadway, egoid::Int)
    ego_index = get_index_of_first_vehicle_with_id(scene, egoid)
    track_longitudinal!(model, scene, roadway, ego_index, 0)
    model
end
function Base.rand(model::PrincetonLongitudinalDriver)
    if isnan(model.σ) || model.σ ≤ 0.0
        model.a
    else
        rand(Normal(model.a, model.σ))
    end
end
function Distributions.pdf(model::PrincetonLongitudinalDriver, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a_lon)
    end
end
function Distributions.logpdf(model::PrincetonLongitudinalDriver, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a_lon)
    end
end

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

"""
Commonly referred to as IDM
"""
type IntelligentDriverModel <: LongitudinalDriverModel
    a::Float64 # predicted acceleration
    σ::Float64 # optional stdev on top of the model, set to zero or NaN for deterministic behavior

    k_spd::Float64 # proportional constant for speed tracking when in freeflow [s⁻¹]

    δ::Float64 # acceleration exponent [-]
    T::Float64 # desired time headway [s]
    v_des::Float64 # desired speed [m/s]
    s_min::Float64 # minimum acceptable gap [m]
    a_max::Float64 # maximum acceleration ability [m/s²]
    d_cmf::Float64 # comfortable deceleration [m/s²] (positive)
    d_max::Float64 # maximum decelleration [m/s²] (positive)

    function IntelligentDriverModel(;
        σ::Float64     =   NaN,
        k_spd::Float64 =   1.0,
        δ::Float64     =   4.0,
        T::Float64     =   1.5,
        v_des::Float64 =  29.0, # typically overwritten
        s_min::Float64 =   5.0,
        a_max::Float64 =   3.0,
        d_cmf::Float64 =   2.0,
        d_max::Float64 =   9.0,
        )

        retval = new()
        retval.a     = NaN
        retval.σ     = σ
        retval.k_spd = k_spd
        retval.δ     = δ
        retval.T     = T
        retval.v_des = v_des
        retval.s_min = s_min
        retval.a_max = a_max
        retval.d_cmf = d_cmf
        retval.d_max = d_max
        retval
    end
end
get_name(::IntelligentDriverModel) = "IDM"
function set_desired_speed!(model::IntelligentDriverModel, v_des::Float64)
    model.v_des = v_des
    model
end
function track_longitudinal!(model::IntelligentDriverModel, scene::Scene, roadway::Roadway, ego_index::Int, target_index::Int)
    veh_ego = scene[ego_index]
    v = veh_ego.state.v

    if target_index > 0
        veh_target = scene[target_index]

        s_gap = get_frenet_relative_position(get_rear_center(veh_target),
                                             veh_ego.state.posF.roadind, roadway).Δs

        if s_gap > 0.0
            Δv = veh_target.state.v - v
            s_des = model.s_min + v*model.T - v*Δv / (2*sqrt(model.a_max*model.d_cmf))
            v_ratio = model.v_des > 0.0 ? (v/model.v_des) : 1.0
            model.a = model.a_max * (1.0 - v_ratio^model.δ - (s_des/s_gap)^2)
        elseif s_gap > -veh_ego.def.length
            model.a = -model.d_max
        else
            Δv = model.v_des - v
            model.a = Δv*model.k_spd
        end

        if isnan(model.a)

            warn("IDM acceleration was NaN!")
            if s_gap > 0.0
                Δv = veh_target.state.v - v
                s_des = model.s_min + v*model.T - v*Δv / (2*sqrt(model.a_max*model.d_cmf))
                println("\tΔv: ", Δv)
                println("\ts_des: ", s_des)
                println("\tv_des: ", model.v_des)
                println("\tδ: ", model.δ)
                println("\ts_gap: ", s_gap)
            elseif s_gap > -veh_ego.def.length
                println("\td_max: ", model.d_max)
            end

            model.a = 0.0
        end
    else
        # no lead vehicle, just drive to match desired speed
        Δv = model.v_des - v
        model.a = Δv*model.k_spd # predicted accel to match target speed
    end

    model.a = clamp(model.a, -model.d_max, model.a_max)

    model
end
function observe!(model::IntelligentDriverModel, scene::Scene, roadway::Roadway, egoid::Int)

    # update the predicted accel

    vehicle_index = get_index_of_first_vehicle_with_id(scene, egoid)
    veh_ego = scene[vehicle_index]
    fore_res = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())

    track_longitudinal!(model, scene, roadway, vehicle_index, fore_res.ind)
    model
end
function Base.rand(model::IntelligentDriverModel)
    if isnan(model.σ) || model.σ ≤ 0.0
        model.a
    else
        rand(Normal(model.a, model.σ))
    end
end
function Distributions.pdf(model::IntelligentDriverModel, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a_lon)
    end
end
function Distributions.logpdf(model::IntelligentDriverModel, a_lon::Float64)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a_lon)
    end
end


text:
""" Delegates loading of code specifying types essential to the package and other packages that use it.  """

include("instruments.jl")
export AbstractInstrument, AbstractInstrument1D, AbstractInstrument2D
#export Generic1D, Generic2D   # WARNING: Might remove these in future

include("spectra.jl")
export AbstractSpectra, AbstractSpectra1D, AbstractSpectra2D
export Spectra1DBasic, Spectra2DBasic

include("chunks.jl")
export AbstractChunkOfSpectrum, ChunkOfSpectrum

include("chunk_list.jl")
export AbstractChunkList, ChunkList
export AbstractChunkListTimeseries, ChunkListTimeseries
export length, num_chunks, num_times
export extract_chunklist_timeseries_with_subset_obs

include("spectra_timeseries_common_wavelengths.jl")
export AbstractSpectralTimeSeriesCommonWavelengths
export SpectralTimeSeriesCommonWavelengths
#=
# WARNING: Still experimental
#export make_vec_metadata_from_spectral_timeseries
export make_spectral_time_series_common_wavelengths_with_selected_times
=#

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
@testset "11.container-with-most-water.jl" begin
    @test max_area_of_container([1, 8, 6, 2, 5, 4, 8, 3, 7]) == 49
    @test max_area_of_container([1, 1]) == 1
end


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, SpecialFunctions, Test
    import ApproxFunBase: testbandedoperator

@testset "Hermite and GaussWeight" begin
    @testset "Evaluation" begin
        f = Fun(x-> x + x^2, Hermite())
        @test f(1.0) ≈ 2.0
        @test values(f) ≈ points(f) + points(f).^2

        w = Fun(GaussWeight(), [1.0])
        @test w(0.1) ≈ exp(-0.1^2)
        w = Fun(GaussWeight(), Float64[])
        @test w(1) == 0

        w = Fun(GaussWeight(), Float64[1.0])
        f = Fun(x-> 1 + x + x^2, Hermite()) * w
        @test f(0.1) ≈ exp(-0.1^2) * (1+0.1+0.1^2)

        f = Fun(x-> 1 + x + x^2, Hermite(2)) * w
        @test f(0.1) ≈ exp(-0.1^2) * (1+0.1+0.1^2)

        w = Fun(GaussWeight(Hermite(2), 0), [1.0,2.0,3.0]);
        w̃ = Fun(Hermite(2), [1.0,2.0,3.0]);
        @test w(0.1) == w̃(0.1)

        @test ApproxFunOrthogonalPolynomials.Recurrence(Hermite())[1:10,1:10]/sqrt(2) ≈ ApproxFunOrthogonalPolynomials.Recurrence(Hermite(2))[1:10,1:10]

        @test points(Hermite(),10) == sqrt(2)points(Hermite(2),10)

        @test Fun(Hermite(),[1])(0.5) == 1
        @test Fun(Hermite(),[0,1])(0.5) == 2*0.5

        @test Fun(Hermite(2),[1])(0.5) == 1
        @test Fun(Hermite(2),[0,1])(0.5) == 2*sqrt(2)*0.5

        @test Fun(x-> 1 + x + x^2, Hermite())(0.5) ≈ Fun(x-> 1 + x + x^2, Hermite(2))(0.5) ≈ 1.75

        p = points(Hermite(),3)
        @test values(Fun(x-> 1 + x + x^2, Hermite())) ≈ 1 .+ p .+ p.^2

        p = points(Hermite(2),3)
        @test values(Fun(x-> 1 + x + x^2, Hermite(2))) ≈ 1 .+ p .+ p.^2

        w = Fun(GaussWeight(Hermite(2), 2), Float64[1.0])
        f = Fun(x-> 1 + x + x^2, Hermite()) * w
        @test space(f) == space(w)
        @test f(0.1) ≈ exp(-2*0.1^2) * (1+0.1+0.1^2)
        f = Fun(x-> 1 + x + x^2, Hermite(2)) * w
        @test space(f) == space(w)
        @test f(0.1) ≈ exp(-2*0.1^2) * (1+0.1+0.1^2)

        L = 1.3; x = 1.2;
        H₀ = Fun(Hermite(), [1.0])
        H̃₀ = Fun(Hermite(L), [1.0])
        @test H̃₀(x) ≈ H₀(sqrt(L) * x)

        L = 1.3; x = 1.2;
        H₁ = Fun(Hermite(), [0.0,1.0])
        H̃₁ = Fun(Hermite(L), [0.0,1.0])
        @test H̃₁(x) ≈ H₁(sqrt(L) * x)
    end

    @testset "Derivative" begin
        D = Derivative(Hermite())
        testbandedoperator(D)

        f = Fun( x-> x + x^2, Hermite())
        g = D * f
        @test g(1.) ≈ 3.
    end
end