text:
"""
    translate(i::Integer, cinfo::CodeInfo, item)

Get a wat string from item, specialized on item type.
"""
translate(i::Integer, ci::CodeInfo, item) = item
translate(i::Integer, ci::CodeInfo, item::AbstractFloat) = "(f32.const $item)"
translate(i::Integer, ci::CodeInfo, item::Number) = "(i32.const $item)"
translate(i::Integer, ci::CodeInfo, item::Bool) = item ? "(i32.const 1)" : "(i32.const 0)"
#=
function translate(i::Integer, ci::CodeInfo, item::Number, slottype=Integer)
    if ci.ssavaluetypes[i] <: AbstractFloat
        return "(f32.const $item)"
    else
        return slottype <: Integer ? "(i32.const $item)" : "(f32.const $item)"
    end
end
translate(i::Integer, ci::CodeInfo, item::Bool, slottype=Bool) = item ? "(i32.const 1)" : "(i32.const 0)"

=#
# translate(i::Integer, ci::CodeInfo, item::Nothing) = "(i32.const 0)"

translate(i::Integer, ci::CodeInfo, item::SlotNumber) = "(local.get \$$(ci.slotnames[item.id]))"
translate(i::Integer ,ci::CodeInfo, item::GlobalRef) = "call \$$(item.name)"
translate(i::Integer ,ci::CodeInfo, item::NewvarNode) = nothing
translate(i::Integer, ci::CodeInfo, item::Const) = translate(i, ci, item.val)
translate(i::Integer, ci::CodeInfo, item::ReturnNode) = "return $(translate(i, ci, item.val))"

function translate(i::Integer, ci::CodeInfo, item::TypedSlot)
    #println("translate TypedSlot, item:",item)
    if isa(item.typ, Const)
        #println("translate TypedSlot, item.typ:",item.typ)
        return translate(i, ci, item.typ.val[1])
    else
        return item
    end
end

"""
    translate(ci::CodeInfo, items::Array)

Get a wat string from items, branching to special cases based on items[1].
"""
function translate(i::Integer, ci::CodeInfo, items::Array)
    # dont translate?
    if isnothing(items)
        return nothing

    # builtin to webassembly?
    elseif is_floatop_simple(ci, items)
        return ["$(floatops[items[1].name])"; translate(i, ci, items[2:end])]
    elseif is_intop_simple(ci, items)
        return ["$(intops[items[1].name])"; translate(i, ci, items[2:end])]

    # special expression?
    elseif hasname(items[1], :(/)) && length(items)==3
        # int/int must become floats and return float. ( float/float and div(int,int) are already taken care of with builtins)
        return ["f32.div", "(f32.convert_i32_s", translate(i, ci, items[2]), ") (f32.convert_i32_s", translate(i, ci, items[3]),")"]

    elseif hasname(items[1], :(size)) && length(items)==3
        return ["call \$size$(items[3])"; translate(i, ci, items[2])]
    elseif hasname(items[1], :(tuple))
        #ignore construction of tuples...
        return [translate.((i,), (ci,), items[2:end])]
    elseif hasname(items[1], :(getindex)) && length(items)==3
        if isa(items[2],SlotNumber) && istuple(ci.slottypes[items[2].id])
            #allow indexing into tuples
            id = items[2].id
            fieldnumber = items[3]
            name = fieldnumber == 1 ? "$(ci.slotnames[id])" : "$(ci.slotnames[id])$(fieldnumber)"
            return "(local.get \$$name)"
        else
            return ["call \$getlinearindex"; translate(i, ci, items[2:end])]
        end
    #=
    elseif hasname(items[1], :(vect))
        len = length(items)-1
    =#
    elseif hasname(items[1], :(setindex!)) && length(items)==4
        #if itemtype(ci, items[3]) <: Integer
        if hasitemtype(ci, items[3], [Integer, Bool])
            return ["call \$setlinearindex_int"; translate(i, ci, items[2:end])]
        else
            return ["call \$setlinearindex"; translate(i, ci, items[2:end])]
        end
    elseif hasname(items[1], :(ifelse))
        return ["select"; translate(i, ci, items[2:end])]
    elseif hasname(items[1], Symbol("return"))
        #handle both 'return' and 'return nothing' (and simply 'nothing' at end of function) as simply return
        if hasname(items[2], Symbol("nothing")) || isnothing(items[2])
            return ["return"]
        elseif isa(items[2], Array) && hasname(items[2][1], :(tuple))
            return ["return"; translate(i, ci, items[2][2:end])]
        else
            return ["return"; translate(i, ci, items[2:end])]
        end
    elseif hasname(items[1], :(not_int))
        return ["i32.eqz", translate(i, ci, items[2])]
    elseif hasname(items[1], Symbol("==="))
        if typeof(items[3]) <: Nothing
            if isa(items[2], SlotNumber) && string(ci.slotnames[items[2].id])[1] == '_' # iterator variable begins with "_"
                return ["i32.eqz", "(local.get \$$(ci.slotnames[items[2].id])i)"] # the iteratorindex
            elseif itemtype(ci, items[2]) <: AbstractFloat
                return ["f32.eq", translate(i, ci, items[2]), "(f32.const 0.0)"]
            else
                return ["i32.eqz", translate(i, ci, items[2])]
            end
        else
            return ["i32.eq"; translate(i, ci, items[2:3])]
        end
    # elseif hasname(items[1], :(=))
    #    slotname = ci.slotnames[items[2].id]
    #    return ["local.set \$$(slotname)", translate(i, ci, items[3])]

    elseif hasname(items[1], :(=))
        slotname = ci.slotnames[items[2].id]
        slottype = ci.slottypes[items[2].id]
        if isa(items[3], Array) && is_iterate(items[3][1])
            return ["local.set \$$(slotname) ( local.set \$$(slotname)i", translate(i, ci, items[3]), ")"]
        elseif isa(items[3], Array) && istuple(slottype)
            settuples = []
            closetuples = []
            for n = 2:length(slottype.parameters)
                push!(settuples,"(local.set \$$(slotname)$(n)")
                push!(closetuples, ")")
            end
            settuples = spacedjoin(settuples)
            closetuples = spacedjoin(closetuples)
            return ["local.set \$$(slotname)",settuples, translate(i, ci, items[3]), closetuples]
        else
            #here maybe check for the type of the slotname, if its float then use f32 inside translate?
            #if (typeof(items[3]) <: Number)
            #    return ["local.set \$$(slotname)", translate(i, ci, items[3], slottype)]
            #else
                return ["local.set \$$(slotname)", translate(i, ci, items[3])]
            #end
        end
    elseif hasname(items[1], :(println))
        if itemtype(ci, items[2]) <: String
            w1,w2 = inlined_tmparray(items[2])
            return [string(w1, w2, "(call \$console_log (i32.const 4))")] #concat strings
        else
            return ["(call \$console_log", translate(i, ci, items[2]), ")"]
        end
    elseif hasname(items[1], :(error))
        if itemtype(ci, items[2]) <: String
            w1,w2 = inlined_tmparray(items[2])
            return [string(w1, w2, "(call \$console_error (i32.const 4))\n(unreachable)")] #concat strings
        else
            return ["(call \$console_error", translate(i, ci, items[2]), ")\n(unreachable)"]
        end

    # default individual translation
    else
        return translate.((i,), (ci,), items)
    end
end

function inlined_tmparray(str::String)
    #also write "unreachable" if throw() and maybe even abort function ("return" in wasm) if error() is called.
    #in julia, error is supposed to abort ececution in parent function aswell?

    v = [Int32(c) for c in str] #convert ASCII/Unicode chars to i32
    watstring1 = "(call \$setsize_tmp (i32.const $(length(v))) (i32.const 1))\n"
    watstring2 = join(["(call \$setlinearindex_int_tmp (i32.const $x) (i32.const $i)) ;;$(Char(x))\n" for (i, x) in enumerate(v)], "")

    return watstring1, watstring2
end

is_iterate(item) = isa(item, GlobalRef) && length(string(item.name)) >= 7 && string(item.name)[1:7] == "iterate"
is_floatop_simple(ci::CodeInfo, items) = isa(items[1], GlobalRef) && items[1].name in keys(floatops) && (hasitemtype(ci, items[2], AbstractFloat) || (length(items) > 2 && hasitemtype(ci, items[3], AbstractFloat)))
is_intop_simple(ci::CodeInfo, items) = isa(items[1], GlobalRef) && items[1].name in keys(intops) && (hasitemtype(ci, items[2], [Integer, Bool, Nothing]) || (length(items) > 2 && hasitemtype(ci, items[3], [Integer, Bool, Nothing])))

floatops = Dict(
:(+) => "f32.add", # these consume 2 args
:(-) => "f32.sub",
:(*) => "f32.mul",
:(/) => "f32.div",
:(==) => "f32.eq",
:(!=) => "f32.ne",
:(<) => "f32.lt",
:(>) => "f32.gt",
:(<=) => "f32.le",
:(>=) => "f32.ge",
:(min) => "f32.min",
:(max) => "f32.max",
:(copysign) => "f32.copysign",
:(abs) => "f32.abs", # these consume 1 arg
:(ceil) => "f32.ceil",
:(floor) => "f32.floor",
:(trunc) => "f32.trunc",
:(round) => "f32.nearest",
:(sqrt) => "f32.sqrt",
:(float) => "f32.convert_i32_s",
:(Int) => "i32.trunc_f32_s",
# :(^) => ["call \$pow",2],
# ""=>["f32.neg",1],
# ""=>["f32.load",1],
# ""=>["f32.store",2],
# ""=>["f32.const",1],
)

intops = Dict(
:(+) => "i32.add", # these consume 2 args
:(-) => "i32.sub",
:(*) => "i32.mul",
:(div) => "i32.div_s", # julia 4/2 will return float. but div(4,2) wont
# :(/) => "i32.div_s", # julia 4/2 will return float. but div(4,2) wont
:(%) => "i32.rem_s",
:(mod) => "i32.rem_s",
:(rem) => "i32.rem_s",
:(&&) => "i32.and",
:(||) => "i32.or",
:(⊻) => "i32.xor",
:(xor) => "i32.xor",
:(<<) => "i32.shl",
:(>>) => "i32.shr_s",
:(==) => "i32.eq",
# :(===) => "i32.eq", #handle this in a special way.
:(!=) => "i32.ne",
:(<) => "i32.lt_s",
:(>) => "i32.gt_s",
:(<=) => "i32.le_s",
:(>=) => "i32.ge_s",
:(!) => "i32.eqz", # these consume 1 arg
:(not_int) => "i32.eqz",
:(leading_zeros) => "i32.clz",
:(trailing_zeros) => "i32.ctz",
:(count_ones) => "i32.popcnt",
:(float) => "f32.convert_i32_s",
:(Int) => "i32.trunc_f32_s",
# :(^) => ["call \$powi",2],
# ""=>["i32.load",1]
# ""=>["i32.store",2]
# ""=>["i32.const",1]
)


text:
struct BaseVertex{T}
    id::Int64

    edges::Vector{T}
    requireds::Vector{T}
end

ext: jl


text:
# SteinProbitRegressionGaussianPrior
#
# We assume we have the following model:
#
# y_i = \Psi (beta_1 x_1 + ... + beta_p x_p)
#
# where \Psi is assumed to be the noraml CDF (and
# a Gaussian prior for beta).

using Distributions: Normal
import Distributions

type SteinProbitRegressionGaussianPrior <: SteinPosterior
    # the X data
    X::Array{Float64,2}
    # the 0/1 y values
    y::Array{Float64,1}
    # priorsigma2 is prior variance
    priorsigma2::Float64
end

defaultsigma2 = 0.1

SteinProbitRegressionGaussianPrior(X, y) =
    SteinProbitRegressionGaussianPrior(X, y, defaultsigma2)

# Returns lower bound of support
function supportlowerbound(d::SteinProbitRegressionGaussianPrior, i::Int)
    -Inf
end

# Returns upper bound of support
function supportupperbound(d::SteinProbitRegressionGaussianPrior, i::Int)
    Inf
end

# Returns the log prior log pi(beta)
function logprior(d::SteinProbitRegressionGaussianPrior,
                  beta::Array{Float64, 1})
    logpriorvalue =
        -0.5 * sum((beta .^ 2) ./ d.priorsigma2) -
        -0.5 * sum(log(2 * pi .* d.priorsigma2))
    logpriorvalue
end

# Returns log pi(d.y[idx] | beta)
function loglikelihood(d::SteinProbitRegressionGaussianPrior,
                       beta::Array{Float64, 1};
                       idx=1:length(d.y))
    y = d.y[idx]
    X = d.X[idx,:]
    ypos = find(y .> 0.0)
    yneg = find(y .<= 0.0)
    u = vec(X * beta)
    gsn = Normal()

    logprobpos = Distributions.logcdf(gsn, u[ypos])
    logprobneg = log(1.0 - Distributions.cdf(gsn, u[yneg]))
    sum(logprobpos) + sum(logprobneg)
end

# Returns the log prior log pi(beta)
function gradlogprior(d::SteinProbitRegressionGaussianPrior,
                      beta::Array{Float64, 1})
    -beta ./ d.priorsigma2
end

# Returns grad_{beta} log pi(d.y[idx] | beta)
function gradloglikelihood(d::SteinProbitRegressionGaussianPrior,
                           beta::Array{Float64,1};
                           idx=1:length(d.y))
    y = d.y[idx]
    X = d.X[idx,:]
    ypos = find(y .> 0.0)
    yneg = find(y .<= 0.0)
    u = vec(X * beta)
    gsn = Normal()

    posweights = Distributions.pdf(gsn, u[ypos]) ./ Distributions.cdf(gsn, u[ypos])
    negweights = -Distributions.pdf(gsn, u[yneg]) ./ (1.0 - Distributions.cdf(gsn, u[yneg]))

    xpos = vec(sum(broadcast(*, X[ypos,:], posweights), 1))
    xneg = vec(sum(broadcast(*, X[yneg,:], negweights), 1))

    xneg + xpos
end

function numdimensions(d::SteinProbitRegressionGaussianPrior)
    size(d.X, 2)
end


text:
using Test
using UnPack
using MOONS.Grids
using MOONS.Fields
using MOONS.FFTs

@testset "FFTs" begin

    # Solve ∇² u = f
    # where f = -2*(3*π)^2*cos(3*π*x)*cos(3*π*y)
    #       u = cos(3*π*x)*cos(3*π*y))

    FT = Float64
    n = 128
    cx = Coordinates(FT(0),FT(1),n)
    cy = Coordinates(FT(0),FT(1),n)
    cz = Coordinates(FT(0),FT(1),2)
    g = Grid((cx,cy,cz))

    f_exact = CellCenter(g)
    u_exact = CellCenter(g)
    f = CellCenter(g)
    u = CellCenter(g)
    u_err = CellCenter(g)

    @inbounds for local_grid in GridField(g,f_exact)
        @unpack ijk,x,y = local_grid
        u_exact[ijk] = cos(3*π*x)*cos(3*π*y)
        f_exact[ijk] = -2*(3*π)^2*cos(3*π*x)*cos(3*π*y)
    end

    ∇²!(f, u_exact, g)
    extrap!(f, g)
    solve!(u, f_exact, g, 3)

    @test max(average(u, 3)...) ≈ 1.0003012584260522
    @test max(abs.(diff(u; dims=3))...) ≈ 0

    u_err .= abs.(u_exact - u);
    assign!(u_err, 0, Ghost())
    @test max(u_err...) ≈ 0.0004518490779501505
end


text:
module Renderer

import ..@fullLibraryPath
import ..Vector2
import ..Vector3
import ..Vector4
import ..Matrix4x4

import ..Electron, ..Proton, ..Neutron
import ..particleRadius

import ..Menu

# Bindings

const Window = Ptr{Cvoid}
const RenderInfo = Ptr{Cvoid}

# Particle drawing
function drawSphere(renderer, radius, color, position)
  ccall((:renderSphere, @fullLibraryPath), Cvoid, (RenderInfo, Cfloat, Vector3, Vector3), renderer, radius, color, position)
end

function color(::Electron)
  return Vector3((0.0, 0.0, 0.8))
end

function color(::Proton)
  return Vector3((0.8, 0.0, 0.0))
end

function color(::Neutron)
  return Vector3((0.2, 0.2, 0.2))
end

function drawParticle(renderer, particle::Electron)
  drawSphere(renderer, particleRadius(particle), color(particle), particle.position)
end

function drawParticle(renderer, particle::Proton)
  drawSphere(renderer, particleRadius(particle), color(particle), particle.position)
end

function drawParticle(renderer, particle::Neutron)
  drawSphere(renderer, particleRadius(particle), color(particle), particle.position)
end

function drawMenu(renderer, menu)
  return ccall((:renderMenu, @fullLibraryPath), Cvoid, (RenderInfo, Menu), renderer, menu)
end

function init()
  return ccall((:init, @fullLibraryPath), Int32, ())
end

function createWindow(width::Int, height::Int, title::String)
  return ccall((:createWindow, @fullLibraryPath), Window, (Int32, Int32, Cstring), width, height, title)
end

function windowCloseStatus(window::Window)
  return ccall((:windowCloseStatus, @fullLibraryPath), Bool, (Window,), window)
end

function closeWindow(window::Window)
  return ccall((:closeWindow, @fullLibraryPath), Cvoid, (Window,), window)
end

function terminate(renderer::RenderInfo)
  ccall((:terminate, @fullLibraryPath), Cvoid, (RenderInfo,), renderer)
end

function createRenderer(window::Window)
  ccall((:createRenderer, @fullLibraryPath), RenderInfo, (Window,), window)
end

function render(renderer::RenderInfo, simulation, menus)
  # Begin Render
  ccall((:beginRender, @fullLibraryPath), Cvoid, (RenderInfo,), renderer)

  # Iterate over particles
  for particle in simulation.objects
    drawParticle(renderer, particle)
  end

  # Iterate over menus
  for menu in menus
    ccall((:updateMenu, @fullLibraryPath), Cvoid, (Menu, RenderInfo), menu, renderer)

    drawMenu(renderer, menu)
  end

  ccall((:renderOrientation, @fullLibraryPath), Cvoid, (RenderInfo,), renderer)

  ccall((:renderGrid, @fullLibraryPath), Cvoid, (RenderInfo,), renderer)

  # End Render
  ccall((:endRender, @fullLibraryPath), Cvoid, (RenderInfo,), renderer)
end

function clear(window::Window)
  ccall((:clear, @fullLibraryPath), Cvoid, (Window,), window)
end

function swapBuffers(window::Window)
  ccall((:swapBuffers, @fullLibraryPath), Cvoid, (Window,), window)
end

function setClearColor(window::Window, red::Float64, green::Float64, blue::Float64, alpha::Float64)
  ccall((:setClearColor, @fullLibraryPath), Cvoid, (Window, Cfloat, Cfloat, Cfloat, Cfloat), window, red, green, blue, alpha)
end

function addKeyCallback(renderer::RenderInfo, callback, key, action)
  ccall((:addKeyCallback, @fullLibraryPath), Cvoid, (RenderInfo, Cint, Cint, Ptr{Cvoid}), renderer, key, action, @cfunction($callback, Cvoid, (RenderInfo, Cint)))
end

function setResizeCallback(renderer, callback)
  ccall((:setResizeCallback, @fullLibraryPath), Cvoid, (RenderInfo, Ptr{Cvoid}), renderer, @cfunction($callback, Cvoid, (RenderInfo, Cint, Cint)))
end

function setClickCallback(renderer, callback)
  ccall((:setClickCallback, @fullLibraryPath), Cvoid, (RenderInfo, Ptr{Cvoid}), renderer, @cfunction($callback, Cvoid, (RenderInfo, Cint, Cint, Cint)))
end

function getCursorPosition(renderer)
  return ccall((:getCursorPosition, @fullLibraryPath), Vector2, (RenderInfo,), renderer)
end

function getWindowSize(renderer)
  return ccall((:getWindowSize, @fullLibraryPath), Vector2, (RenderInfo,), renderer)
end

@enum Key begin
  Unknown = -1
  Space = 32
  Apostrophe = 39
  Comma = 44
  Minus = 45
  Period = 46
  Slash = 47
  Key0 = 48
  Key1 = 49
  Key2 = 50
  Key3 = 51
  Key4 = 52
  Key5 = 53
  Key6 = 54
  Key7 = 55
  Key8 = 56
  Key9 = 57
  Semicolon = 59
  Equal = 61
  KeyA = 65
  KeyB = 66
  KeyC = 67
  KeyD = 68
  KeyE = 69
  KeyF = 70
  KeyG = 71
  KeyH = 72
  KeyI = 73
  KeyJ = 74
  KeyK = 75
  KeyL = 76
  KeyM = 77
  KeyN = 78
  KeyO = 79
  KeyP = 80
  KeyQ = 81
  KeyR = 82
  KeyS = 83
  KeyT = 84
  KeyU = 85
  KeyV = 86
  KeyW = 87
  KeyX = 88
  KeyY = 89
  KeyZ = 90
  LeftBracket = 91
  Backslash = 92
  RightBracket = 93
  GraveAccent = 96
  World1 = 161
  World2 = 162
  Escape = 256
  Enter = 257
  Tab = 258
  Backspace = 259
  Insert = 260
  Delete = 261
  KeyRight = 262
  KeyLeft = 263
  KeyDown = 264
  KeyUp = 265
  PageUp = 266
  PageDown = 267
  KeyHome = 268
  KeyEnd = 269
  KeyCapsLock = 280
  KeyScrollLock = 281
  KeyNumLock = 282
  KeyPrintScreen = 283
  KeyPause = 284
  KeyF1 = 290
  KeyF2 = 291
  KeyF3 = 292
  KeyF4 = 293
  KeyF5 = 294
  KeyF6 = 295
  KeyF7 = 296
  KeyF8 = 297
  KeyF9 = 298
  KeyF10 = 299
  KeyF11 = 300
  KeyF12 = 301
  KeyF13 = 302
  KeyF14 = 303
  KeyF15 = 304
  KeyF16 = 305
  KeyF17 = 306
  KeyF18 = 307
  KeyF19 = 308
  KeyF20 = 209
  KeyF21 = 310
  KeyF22 = 311
  KeyF23 = 312
  KeyF24 = 313
  KeyF25 = 314
  KeyPad0 = 320
  KeyPad1 = 321
  KeyPad2 = 322
  KeyPad3 = 323
  KeyPad4 = 324
  KeyPad5 = 325
  KeyPad6 = 326
  KeyPad7 = 327
  KeyPad8 = 328
  KeyPad9 = 329
  KeyPadDecimal = 330
  KeyPadDivide = 331
  KeyPadMultiply = 332
  KeyPadSubtract = 333
  KeyPadAdd = 334
  KeyPadEnter = 335
  KeyPadEqual = 336
  LeftShift = 340
  LeftControl = 341
  LeftAlt = 342
  LeftSuper = 343
  RightShift = 344
  RightControl = 345
  RightAlt = 346
  RightSuper = 347
  KeyMenu = 348
end

@enum KeyAction begin
  Release = 0
  Press = 1
  Repeat = 2
end

@enum MouseButton begin
  MouseButtonLast = 7
  MouseButtonLeft = 0
  MouseButtonMiddle = 2
  MouseButtonRight = 1
end

const Camera = Ptr{Cvoid}

function getCamera(renderer::RenderInfo)
  return ccall((:getCamera, @fullLibraryPath), Camera, (RenderInfo,), renderer)
end

function cameraGetPosition(camera::Camera)
  return ccall((:cameraGetPosition, @fullLibraryPath), Vector3, (Camera,), camera)
end

function cameraGetTarget(camera::Camera)
  return ccall((:cameraGetTarget, @fullLibraryPath), Vector3, (Camera,), camera)
end

function cameraGetUp(camera::Camera)
  return ccall((:cameraGetUp, @fullLibraryPath), Vector3, (Camera,), camera)
end

function cameraSetPosition(camera::Camera, position::Vector3)
  return ccall((:cameraSetPosition, @fullLibraryPath), Cvoid, (Camera, Vector3), camera, position)
end

function cameraSetTarget(camera::Camera, target::Vector3)
  return ccall((:cameraSetTarget, @fullLibraryPath), Cvoid, (Camera, Vector3), camera, target)
end

function cameraSetUp(camera::Camera, up::Vector3)
  return ccall((:cameraSetUp, @fullLibraryPath), Cvoid, (Camera, Vector3), camera, up)
end

function projectionMatrix(renderer)
  return ccall((:projectionMatrix, @fullLibraryPath), Matrix4x4, (RenderInfo,), renderer)
end

function viewMatrix(camera)
  return ccall((:getViewMatrix, @fullLibraryPath), Matrix4x4, (Camera,), camera)
end

import ..invertMatrix
import ..normalize
import ..dot

function pickObject(renderer::RenderInfo, simulation)
  cursor = getCursorPosition(renderer)
  window = getWindowSize(renderer)
  camera = getCamera(renderer)

  ray = Vector3(((2.0 * cursor.x) / window.x) - 1.0,
    ((2.0 * cursor.y) / window.y) - 1.0,
    1.0)
  clipRay = Vector4(ray.x, ray.y, -1.0, 1.0)

  projection = projectionMatrix(renderer)
  inverseProjection = invertMatrix(projection)
  eyeRay = inverseProjection * clipRay
  eyeRay = Vector4(eyeRay.x, eyeRay.y, -1.0, 0.0)

  view = viewMatrix(camera)
  inverseView = invertMatrix(view)
  worldRay = inverseView * eyeRay
  worldRay = normalize(Vector3(worldRay.x, worldRay.y, worldRay.z))

  for i in 1:length(simulation.objects)
    object = simulation.objects[i]

    # Ray
    origin = cameraGetPosition(camera)
    direction = worldRay

    # Particle
    center = object.position
    radius = particleRadius(object)

    b = dot(direction, origin - center)
    c = dot(origin - center, origin - center) - (radius * radius)

    discriminant = (b * b) - c
    if discriminant >= 0
      return i
    end
  end
  return -1
end

end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
# implementation of the StandardLinearSolver
# This is the interface algorithms should use for doing linear solves

"""
  Calculates the preconditioner for the linear solver.  Thsi preconditioner
  will be used for all linear solves until this function is called again.

  For direct solvers, this function calculates the linear operator itself.
  Prefer [`calcPCandLO`](@ref) whenever possible.

  **Inputs**

   * ls: StandardLinearSolver
   * mesh
   * sbp
   * eqn
   * opts
   * ctx_residual: the ctx required by [`physicsRhs`](@ref) like functions
   * t: current time

  **Keyword Arguments**

   * start_comm: start parallel communication (if required by the PC), default
                  false.  This means the user is generally required to make sure
                  parallel communication is started before calling this
                  function.
"""
function calcPC(ls::StandardLinearSolver, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t; start_comm=false)

  if typeof(ls.pc) <: PCNone
    if start_comm && needParallelData(ls.lo)
      startSolutionExchange(mesh, sbp, eqn, opts)
    end
    calcLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t)
  else
    if start_comm && needParallelData(ls.pc)
      startSolutionExchange(mesh, sbp, eqn, opts)
    end
    calcPC(ls.pc, mesh, sbp, eqn, opts, ctx_residual, t)
  end

  return nothing
end

"""
  Calculates the linear operator.  Use this function only if you want to
  calculate the linear operator and not the preconditioner.
  Prefer [`calcPCandLO`](@ref), which avoids calculating the matrix twice if
  the preconditioner and linear operator share the same matrix

  **Inputs**

   * ls: StandardLinearSolver
   * mesh
   * sbp
   * eqn
   * opts
   * ctx_residual: the ctx required by [`physicsRhs`](@ref) like functions
   * t: current time

  **Keyword Arguments**

   * start_comm: start parallel communication (if required by the lo), default
                  false.  This means the user is generally required to make sure
                  parallel communication is started before calling this
                  function.

"""
function calcLinearOperator(ls::StandardLinearSolver, mesh::AbstractMesh,
                         sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts::Dict, ctx_residual, t;
                         start_comm=false)

  if start_comm && needParallelData(ls.lo)
    startSolutionExchange(mesh, sbp, eqn, opts)
  end
  calcLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t)

end

"""
  Calculates both the preconditioner and linear operator.  In the case where
  they share the matrix, the calculation is only performed once.  This function
  should be preferred to calling [`calcPC`](@ref) and [`calcLinearOperator`](@ref) one
  after the other

  **Inputs**

   * ls: StandardLinearSolver
   * mesh
   * sbp
   * eqn
   * opts
   * ctx_residual: the ctx required by [`physicsRhs`](@ref) like functions
   * t: current time

  **Keyword Arguments**

   * start_comm: start parallel communication (if required by the lo), default
                  false.  This means the user is generally required to make sure
                  parallel communication is started before calling this
                  function.  Note that it is not possible to interleave
                  communication and computation in this case, so performing
                  communication will be very expensive.

"""
function calcPCandLO(ls::StandardLinearSolver, mesh::AbstractMesh,
                     sbp::AbstractOperator,
                     eqn::AbstractSolutionData, opts::Dict, ctx_residual, t;
                     start_comm=false)

  #TODO: why is wait=true here? shouldn't the physics module wait if it hasn't
  #      been done yet?
  if typeof(ls.pc) <: PCNone
    if start_comm && needParallelData(ls.lo)
      startSolutionExchange(mesh, sbp, eqn, opts, wait=true)
    end
    calcLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t)
  elseif ls.shared_mat
    if start_comm && needParallelData(ls.lo)
      startSolutionExchange(mesh, sbp, eqn, opts, wait=true)
    end
    calcLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t)
    pc2 = getBasePC(ls.pc)
    pc2.is_setup = false  # we don't call calcPC but this flag still needs to
                          # be set
  else
    if start_comm && ( needParallelData(ls.lo) || needParallelData(ls.pc))
      startSolutionExchange(mesh, sbp, eqn, opts, wait=true)
    end
    calcPC(ls.pc, mesh, sbp, eqn, opts, ctx_residual, t)
    calcLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t)
  end

  return nothing
end

"""
  Apply the preconditioner, ie. x = inv(Ap)*b, where Ap is an approximation
  to the linear operator used for preconditioner.  This function also works
  for matrix-free methods.

  Note that [`calcPC`](@ref) or [`calcPCandLO`](@ref) must be called before
  this function can be used.

  For direct methods, a linear solve is performed because there is no
  preconditioner.
"""
function applyPC(ls::StandardLinearSolver, mesh::AbstractMesh, sbp::AbstractOperator,
                 eqn::AbstractSolutionData, opts::Dict, t, b::AbstractVector,
                 x::AbstractVector)

  if typeof(ls.pc) <: PCNone
    # apply the most recent linear operator
    linearSolve(ls, b, x)
  else
    applyPC(ls.pc, mesh, sbp, eqn, opts, t, b, x)
  end

  return nothing
end

"""
  Like [`applyPC`](@ref), but applies the transpose (if possible, otherwise
  throws an error).
"""
function applyPCTranspose(ls::StandardLinearSolver, mesh::AbstractMesh,
                          sbp::AbstractOperator,
                          eqn::AbstractSolutionData, opts::Dict, t,
                          b::AbstractVector, x::AbstractVector)
  if typeof(ls.pc) <: PCNone
    linearSolveTranspose(ls, b, x)
  else
    applyPCTranspose(ls.pc, mesh, sbp, eqn, opts, t, b, x)
  end

  return nothing
end

"""
  Solves the linear system Ax = b for x. This function does not recompute the
  precondtioner or linear operator,  The preconditioner and linear operator
  used are the ones calculated by the most recent call to
  [`calcPC`](@ref), [`calcLinearOperator`](@ref) (or [`calcPCandLO`](@ref).

  For Petsc matrices, this function does the final matrix assembly.

  **Inputs**

   * ls: the [`StandardLinearSolver`](@ref) object.
   * b: the right hand side vector (local process portion only)
   * verbose: verbosity level, default 5, < 4 indicates no output

  **Inputs/Outputs**

   * x: vector overwritten with result (local process portion only)

  Implementation Notes:

  The preconditioner and matrix factorization (for direct solves) might
  be computed during this function, but they will be computed at the state
  correspondong to the last call to the functions that calculate them.
"""
function linearSolve(ls::StandardLinearSolver, b::AbstractVector,
                     x::AbstractVector, verbose=5)

  # call the specific solver

  #TODO: have time_all return a struct
  tmp, t_solve, t_gc, alloc = @time_all _linearSolve(ls, b, x)

  @verbose5 if ls.myrank == 0
    println(BSTDOUT, "matrix solve:")
    print_time_all(BSTDOUT, t_solve, t_gc, alloc)
  end

  return nothing
end

"""
  Similar to [`linearSolve`](@ref), but solves A.'x = b.  See that function
  for details.
"""
function linearSolveTranspose(ls::StandardLinearSolver, b::AbstractVector,
                     x::AbstractVector, verbose=5)

  # call the specific solver

  #TODO: have time_all return a struct
  tmp, t_solve, t_gc, alloc = @time_all _linearSolve(ls, b, x, trans=true)

  @verbose5 if ls.myrank == 0
    println(BSTDOUT, "matrix transpose solve:")
    print_time_all(BSTDOUT, t_solve, t_gc, alloc)
  end

  return nothing
end

function _linearSolve(
ls::StandardLinearSolver{Tpc, Tlo},
b::AbstractVector, x::AbstractVector; trans=false) where {Tlo <: AbstractDenseLO, Tpc}

  @assert typeof(ls.pc) <: PCNone

  lo2 = getBaseLO(ls.lo)

  # factor if needed
  if !getIsSetup(lo2)
    getrf!(lo2.A, lo2.ipiv)
    setIsSetup(lo2, true)
    lo2.nfactorizations += 1
  end

  tchar = trans ? 'T' : 'N'

  # solve
  copy!(x, b)
  info = getrs2!(tchar, lo2.A, lo2.ipiv, x)
  @assert info == 0

  if trans
    lo2.ntsolves += 1
  else
    lo2.nsolves += 1
  end

  return nothing
end

"""
  On Skylake, there is a problem with reusing the symbolic factorization
  with UMFPACK.  If this constant is true, it will recalculate the the
  symbolic factorization each time.  On other systems, this constant should
  be set to false.
"""
global const SKYLAKE_STACKSMASH=true

function _linearSolve(ls::StandardLinearSolver{Tpc, Tlo},
                      b::AbstractVector, x::AbstractVector;
                      trans=false) where {Tlo <: AbstractSparseDirectLO, Tpc}

  @assert typeof(ls.pc) <: PCNone

  lo2 = getBaseLO(ls.lo)

  make_zerobased(lo2.A)
  # compute factorization if needed
  if !getIsSetup(lo2)
    umfpack_free_numeric(lo2.fac)  # free old factorization
    if SKYLAKE_STACKSMASH
      umfpack_free_symbolic(lo2.fac)  # DEBUGGING
    end
    # note: the matrix stored in the factorization object must alias. lo2.A
    umfpack_numeric!(lo2.fac)
    setIsSetup(lo2, true)
    lo2.nfactorizations += 1
  end

  tchar = trans ? UMFPACK_At : UMFPACK_A

  solve_suitesparse(lo2.fac, b, tchar, x)
  make_onebased(lo2.A)
  if trans
    lo2.ntsolves += 1
  else
    lo2.nsolves += 1
  end

  return nothing
end

function _linearSolve(
ls::StandardLinearSolver{Tpc, Tlo},
b::AbstractVector, x::AbstractVector; trans=false) where {Tlo <: PetscLO , Tpc}

  myrank = ls.myrank
  pc2 = getBasePC(ls.pc)
  lo2 = getBaseLO(ls.lo)
  @assert !(typeof(pc2) <: PCNone)

  # prepare the data structures
  #TODO: copy x into xtmp to support initial guess non-zero
  t_assem = @elapsed assemblePetscData(ls, b, lo2.btmp, x, lo2.xtmp)
  @mpi_master println(BSTDOUT, "Final matrix assembly time = ", t_assem)

  if !isPCMatFree(ls) && !pc2.is_setup
    setupPC(pc2)
  end

  # do the solve
  ksp = ls.ksp
  SetTolerances(ksp, ls.reltol, ls.abstol, ls.dtol, PetscInt(ls.itermax))

  if trans
    KSPSolveTranspose(ksp, lo2.btmp, lo2.xtmp)
    lo2.ntsolves += 1
  else
    KSPSolve(ksp, lo2.btmp, lo2.xtmp)
    lo2.nsolves += 1
  end

  # print convergence info
  @mpi_master begin
    reason = GetConvergedReason(ksp)
    println(BSTDOUT, "KSP converged reason = ", KSPConvergedReasonDict[reason])
    rnorm = GetResidualNorm(ksp)
    println("Linear residual = ", rnorm)
  end

  # copy result back to x
  xtmp = VecGetArrayRead(lo2.xtmp)
  copy!(x, xtmp)
  VecRestoreArrayRead(lo2.xtmp, xtmp)

  # Reuse preconditionre until next time setupPC() is called
  PCSetReusePreconditioner(pc2.pc, PETSC_TRUE)

  return nothing
end

"""
  Helper function for doing linear solves with Petsc matrices.  It assembles
  A, and Ap, and copies b into the local part of the petsc vector for b.
"""
function assemblePetscData(ls::StandardLinearSolver, b::AbstractVector,
                           b_petsc::PetscVec, x::AbstractVector,
                           x_petsc::PetscVec)

  lo_matfree = isLOMatFree(ls)
  pc_matfree = isPCMatFree(ls)

  pc2 = getBasePC(ls.pc)
  lo2 = getBaseLO(ls.lo)
  myrank = ls.myrank

  # assemble things
  if !lo_matfree && !getIsSetup(lo2)
    MatAssemblyBegin(lo2.A, MAT_FINAL_ASSEMBLY)
  end

  if !pc_matfree && !ls.shared_mat
    MatAssemblyBegin(pc2.A, MAT_FINAL_ASSEMBLY)
  end

  # copy values into the vector
  btmp = VecGetArray(b_petsc)
  copy!(btmp, b)
  VecRestoreArray(b_petsc, btmp)

  xtmp = VecGetArray(x_petsc)
  copy!(xtmp, x)
  VecRestoreArray(x_petsc, xtmp)

  # end assembly
  if !lo_matfree && !getIsSetup(lo2)
    MatAssemblyEnd(lo2.A, MAT_FINAL_ASSEMBLY)
    setIsSetup(lo2, true)
    lo2.nassemblies[1] += 1
    matinfo = MatGetInfo(lo2.A, PETSc2.MAT_LOCAL)
    if matinfo.mallocs > 0.5  # if any mallocs
      println(BSTDERR, "Warning: non-zero number of mallocs for A on process $myrank: $(matinfo.mallocs) mallocs")
    end
  end

  if !pc_matfree && !ls.shared_mat
    MatAssemblyEnd(pc2.A, MAT_FINAL_ASSEMBLY)
    setIsAssembled(pc2, true)
    pc2.nassemblies[1] += 1
    matinfo = MatGetInfo(pc2.A, PETSc2.MAT_LOCAL)
    if matinfo.mallocs > 0.5  # if any mallocs
      println(BSTDERR, "Warning: non-zero number of mallocs for Ap on process $myrank: $(matinfo.mallocs) mallocs")
    end
  end

#  pc2.is_setup = true

  return nothing
end

"""
  Returns true if the linear operator is matrix free, false otherwise
"""
function isLOMatFree(ls::StandardLinearSolver)

  return typeof(ls.lo) <: AbstractPetscMatFreeLO
end

"""
  Returns true if the preconditioner is matrix free, false otherwise
"""
function isPCMatFree(ls::StandardLinearSolver)

  return typeof(ls.pc) <: AbstractPetscMatFreePC
end

"""
  Set the tolerances for iterative solves.  Has no effect for direct solves.
  Supplying a negative value results in retaining the original value.

  **Inputs**

   * ls: StandardLinearSolver
   * reltol: relative residual tolerance
   * abstol: absolute residual tolerance
   * dtol: divergence tolerance
   * itermax: maximum number of iterations
"""
function setTolerances(ls::StandardLinearSolver, reltol::Number, abstol::Number,
                       dtol::Number, itermax::Integer)

  if reltol > 0
    ls.reltol = reltol
  end

  if abstol > 0
    ls.abstol = abstol
  end

  if dtol > 0
    ls.dtol = dtol
  end

  if itermax > 0
    ls.itermax = itermax
  end

  return nothing
end

"""
  This function frees any memory owned by external libraries, both in the
  StandardLinearSolver object itself and in the pc and lo objects.
  Therefore, at the end of a run, all you need to do is free the
  StandardLinearSolver and everything will be taken care of.

  It is safe to call this function multiple times.
"""
function free(ls::StandardLinearSolver)

  if !ls.is_finalized
    if ls.ksp.pobj != C_NULL
      PetscDestroy(ls.ksp)
      ls.ksp.pobj = C_NULL
    end

    free(ls.pc)
    free(ls.lo)

  end

  ls.is_finalized = true

  return nothing
end


text:
export Config
module Config

println("Configuring Mocha...")
const use_cuda_key             = "MOCHA_USE_CUDA"
const cuda_dev_key             = "MOCHA_CUDA_DEVICE"

const use_native_extension_key = "MOCHA_USE_NATIVE_EXT"

parseEnvInt(key; dflt=0) = begin
  try
    parseint(ENV[key])
  catch
    dflt
  end
end

parseEnvBool(key; dflt=false) = begin
  try
    lowercase(ENV[key]) != "false"
  catch
    dflt
  end
end

const use_cuda = parseEnvBool(use_cuda_key)
const cuda_dev_id = parseEnvInt(cuda_dev_key)

const use_native_extension = parseEnvBool(use_native_extension_key)

use_cuda && println(" * CUDA       enabled [DEV=$cuda_dev_id] ($use_cuda_key environment variable detected)")
use_cuda || println(" * CUDA       disabled by default")

use_native_extension && println(" * Native Ext enabled ($use_native_extension_key environment variable detected)")
use_native_extension || println(" * Native Ext disabled by default")

println("Mocha configured, continue loading module...")

end


text:
using Documenter, PkgBenchmark

makedocs(
    modules = [PkgBenchmark],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "PkgBenchmark.jl",
    pages = Any[
        "Home" => "index.md",
        "define_benchmarks.md",
        "run_benchmarks.md",
        "comparing_commits.md",
        "export_markdown.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaCI/PkgBenchmark.jl.git",
)


text:
# See https://github.com/tkoolen/VariableHeightInvertedPendulum
function cubic_orbital_energy_control_law(g, zf, x, z, xd, zd)
    a = xd / x
    b = zd - a * z
    u = -7 * a^2 + (3 * zf * a^3 - g * a) / b - (10 * a^3 * b) / g
end

mutable struct VariableHeightMomentumController
    τ::Vector{Float64}
    dynamicsResult::DynamicsResult{Float64}
    positionControlledJointReferences::Dict{GenericJoint{Float64}, Vector{Float64}}
    centroidalFrame::CartesianFrame3D
    soleFrame::CartesianFrame3D

    function VariableHeightMomentumController(
            mechanism::Mechanism, positionControlledJoints::Vector{GenericJoint{Float64}},
            centroidalFrame::CartesianFrame3D, soleFrame::CartesianFrame3D)
        τ = zeros(num_velocities(mechanism))
        positionControlledJointReferences = Dict(j => zeros(num_positions(j)) for j in positionControlledJoints)
        new(τ, DynamicsResult{Float64}(mechanism), positionControlledJointReferences, centroidalFrame, soleFrame)
    end
end

function set_desired_configuration!(controller::VariableHeightMomentumController, joint::Joint, qdes::AbstractVector)
    controller.positionControlledJointReferences[joint] .= qdes
end

function control(controller::VariableHeightMomentumController, t, state)
    # create model
    solver = Gurobi.GurobiSolver() # currently overkill; it's just a linear solve right now.
    Gurobi.setparameters!(solver, Silent = true)
    model = Model(solver = solver)
    nv = num_velocities(state)
    @variable(model, v̇[1 : nv])

    # Compute desired rate of change of linear momentum
    centroidalFrame = controller.centroidalFrame
    soleFrame = controller.soleFrame
    mechanism = state.mechanism
    com = center_of_mass(state) - transform(state, Point3D(Float64, soleFrame), root_frame(mechanism))
    v = velocity(state)
    A = momentum_matrix(state)
    m = mass(state.mechanism)
    comdot = FreeVector3D(A.frame, A.linear * v) / m
    g = -RigidBodyDynamics.gravitational_spatial_acceleration(mechanism).linear[3]
    zf = 0.95
    x, z = com.v[1], com.v[3]
    xd, zd = comdot.v[1], comdot.v[3]
    u = cubic_orbital_energy_control_law(g, zf, x, z, xd, zd)
    u = max(u, zero(u))
    fx = m * x * u
    fz = m * (z * u - g)
    ḣdes = Wrench{Float64}(centroidalFrame, zero(SVector{3, eltype(u)}), SVector(fx, zero(u), fz))
    centroidalToWorld = Transform3D(centroidalFrame, root_frame(mechanism), com.v)
    ḣdes = transform(ḣdes, centroidalToWorld)

    # Set up QP
    Ȧv = momentum_rate_bias(state)
    @framecheck A.frame ḣdes.frame
    @framecheck A.frame Ȧv.frame
    @constraint(model, Array(A.linear) * v̇ + Array(Ȧv.linear) .== Array(ḣdes.linear))
    @constraint(model, Array(A.angular) * v̇ + Array(Ȧv.angular) .== Array(ḣdes.angular))

    # PD control for upper body and free leg
    kp = 100.
    kd = 20.
    for (joint, q_desired) in controller.positionControlledJointReferences
        v̇pd = kp .* (q_desired .- configuration(state, joint)) - kd .* velocity(state, joint)
        @constraint(model, v̇[velocity_range(state, joint)] .== v̇pd)
    end

    # minimize squared joint accelerations
    @objective(model, Min, sum(v̇[i]^2 for i = 1 : nv))

    # solve
    status = solve(model)

    # Compute CoP
    gravitationalWrench = transform(Wrench(centroidalFrame, zero(SVector{3}), SVector(0, 0, -m * g)), centroidalToWorld)
    groundReactionWrench = ḣdes - gravitationalWrench
    groundReactionWrench = transform(state, groundReactionWrench, soleFrame)
    cop = center_of_pressure(groundReactionWrench)
    cop = transform(state, cop, root_frame(mechanism))
    @assert all(isnan(cop.v)) || isapprox(cop.v[1], 0., atol = 1e-8)

    # inverse dynamics
    result = controller.dynamicsResult
    result.v̇ = getvalue(v̇)
    inverse_dynamics!(controller.τ, result.jointwrenches, result.accelerations, state, result.v̇)
    controller.τ
end

max_stars_repo_name: tkoolen/VariableHeightMomentumControl

max_issues_repo_name: tkoolen/VariableHeightMomentumControl

max_forks_repo_name: tkoolen/VariableHeightMomentumControl


text:
@testset "benchmark_shape" begin
    Z = ground_truth(SynthSphere(50))
    Z2 = ground_truth(SynthSphere(51))
    error, a = benchmark_shape(Z, Z2)
    @test error ≈ 0.0255463213362199
end


text:
__precompile__()

module IntelVectorMath

export IVM
const IVM = IntelVectorMath

# import Base: .^, ./
# using Libdl
include("../deps/deps.jl")

include("setup.jl")

for t in (Float32, Float64, ComplexF32, ComplexF64)
    # Unary, real or complex
    def_unary_op(t, t, :acos, :acos!, :Acos)
    def_unary_op(t, t, :asin, :asin!, :Asin)
    def_unary_op(t, t, :acosh, :acosh!, :Acosh)
    def_unary_op(t, t, :asinh, :asinh!, :Asinh)
    def_unary_op(t, t, :sqrt, :sqrt!, :Sqrt)
    def_unary_op(t, t, :exp, :exp!, :Exp)
    def_unary_op(t, t, :log, :log!, :Ln)

    # # Binary, real or complex
    def_binary_op(t, t, :pow, :pow!, :Pow, true)
    def_binary_op(t, t, :divide, :divide!, :Div, true)
end

for t in (Float32, Float64)
    # Unary, real-only
    def_unary_op(t, t, :cbrt, :cbrt!, :Cbrt)
    def_unary_op(t, t, :expm1, :expm1!, :Expm1)
    def_unary_op(t, t, :log1p, :log1p!, :Log1p)
    def_unary_op(t, t, :abs, :abs!, :Abs)
    def_unary_op(t, t, :abs2, :abs2!, :Sqr)
    def_unary_op(t, t, :ceil, :ceil!, :Ceil)
    def_unary_op(t, t, :floor, :floor!, :Floor)
    def_unary_op(t, t, :round, :round!, :Round)
    def_unary_op(t, t, :trunc, :trunc!, :Trunc)

    # now in SpecialFunctions (make smart, maybe?)
    def_unary_op(t, t, :erf, :erf!, :Erf)
    def_unary_op(t, t, :erfc, :erfc!, :Erfc)
    def_unary_op(t, t, :erfinv, :erfinv!, :ErfInv)
    def_unary_op(t, t, :erfcinv, :erfcinv!, :ErfcInv)
    def_unary_op(t, t, :lgamma, :lgamma!, :LGamma)
    def_unary_op(t, t, :gamma, :gamma!, :TGamma)
    # Not in Base
    def_unary_op(t, t, :inv_cbrt, :inv_cbrt!, :InvCbrt)
    def_unary_op(t, t, :inv_sqrt, :inv_sqrt!, :InvSqrt)
    def_unary_op(t, t, :pow2o3, :pow2o3!, :Pow2o3)
    def_unary_op(t, t, :pow3o2, :pow3o2!, :Pow3o2)

    # Enabled only for Real. MKL guarantees higher accuracy, but at a
    # substantial performance cost.
    def_unary_op(t, t, :atan, :atan!, :Atan)
    def_unary_op(t, t, :cos, :cos!, :Cos)
    def_unary_op(t, t, :sin, :sin!, :Sin)
    def_unary_op(t, t, :tan, :tan!, :Tan)
    def_unary_op(t, t, :atanh, :atanh!, :Atanh)
    def_unary_op(t, t, :cosh, :cosh!, :Cosh)
    def_unary_op(t, t, :sinh, :sinh!, :Sinh)
    def_unary_op(t, t, :tanh, :tanh!, :Tanh)
    def_unary_op(t, t, :log10, :log10!, :Log10)

    # # .^ to scalar power
    # mklfn = Base.Meta.quot(Symbol("$(vml_prefix(t))Powx"))
    # @eval begin
    #     export pow!
    #     function pow!{N}(out::Array{$t,N}, A::Array{$t,N}, b::$t)
    #         size(out) == size(A) || throw(DimensionMismatch())
    #         ccall(($mklfn, lib), Nothing, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
    #         vml_check_error()
    #         out
    #     end
    #     function (.^){N}(A::Array{$t,N}, b::$t)
    #         out = similar(A)
    #         ccall(($mklfn, lib), Nothing, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
    #         vml_check_error()
    #         out
    #     end
    # end

    # # Binary, real-only
    def_binary_op(t, t, :atan, :atan!, :Atan2, false)
    def_binary_op(t, t, :hypot, :hypot!, :Hypot, false)

    # Unary, complex-only
    def_unary_op(Complex{t}, Complex{t}, :conj, :conj!, :Conj)
    def_unary_op(Complex{t}, t, :abs, :abs!, :Abs)
    def_unary_op(Complex{t}, t, :angle, :angle!, :Arg)

    ### cis is special, IntelVectorMath function is based on output
    def_unary_op(t, Complex{t}, :cis, :cis!, :CIS; vmltype = Complex{t})

    # Binary, complex-only. These are more accurate but performance is
    # either equivalent to Base or slower.
    # def_binary_op(Complex{t}, Complex{t}, (:+), :add!, :Add, false)
    # def_binary_op(Complex{t}, Complex{t}, (:.+), :add!, :Add, true)
    # def_binary_op(Complex{t}, Complex{t}, (:.*), :multiply!, :Mul, true)
    # def_binary_op(Complex{t}, Complex{t}, (:-), :subtract!, :Sub, false)
    # def_binary_op(Complex{t}, Complex{t}, (:.-), :subtract!, :Sub, true)
    # def_binary_op(Complex{t}, Complex{t}, :multiply_conj, :multiply_conj!, :Mul, false)
end

export VML_LA, VML_HA, VML_EP, vml_set_accuracy, vml_get_accuracy

end


text:
@testset "CSV" begin
    function Models.HKAgent(id, op)
        return Models.HKAgent(id, op, op, -1)
    end
    function Models.HKAgent(; id = -1, op1 = -1, op2 = -1, op3 = -1)
        return Models.HKAgent(id, op1, op2, op3)
    end

    function Models.Fighter(id, p1, p2, p3, has, cap, shape)
        return Models.Fighter(id, (p1, p2, p3), has, cap, shape)
    end

    @agent Foo GridAgent{2} begin end
    @agent Bar GridAgent{2} begin end

    model = ABM(Union{Foo,Bar}, GridSpace((5,5)); warn = false)

    @test_throws AssertionError AgentsIO.populate_from_csv!(model, "test.csv")

    model, _ = Models.hk(numagents = 10)
    empty_model, _ = Models.hk(numagents = 0)

    AgentsIO.dump_to_csv("test.csv", allagents(model))

    open("test.csv", "r") do f
        @test length(split(readline(f), ',')) == 4
    end

    AgentsIO.populate_from_csv!(empty_model, "test.csv")

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].old_opinion == empty_model[i].old_opinion for i in allids(model))
    @test all(model[i].new_opinion == empty_model[i].new_opinion for i in allids(model))
    @test all(model[i].previous_opinion == empty_model[i].previous_opinion for i in allids(model))

    genocide!(empty_model)
    AgentsIO.populate_from_csv!(empty_model, "test.csv", Models.HKAgent, Dict(:id => 1, :op1 => 3, :op2 => 2))

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].old_opinion == empty_model[i].old_opinion for i in allids(model))
    @test all(model[i].new_opinion == empty_model[i].new_opinion for i in allids(model))
    @test all(model[i].previous_opinion == empty_model[i].previous_opinion for i in allids(model))

    AgentsIO.dump_to_csv("test.csv", Models.HKAgent[model[i] for i in 1:nagents(model)], [:old_opinion])

    open("test.csv", "r") do f
        @test length(split(readline(f), ',')) == 1
    end

    genocide!(empty_model)
    AgentsIO.populate_from_csv!(empty_model, "test.csv"; row_number_is_id = true)

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].old_opinion == empty_model[i].old_opinion for i in allids(model))
    @test all(model[i].new_opinion == empty_model[i].new_opinion for i in allids(model))
    @test all(model[i].previous_opinion == empty_model[i].previous_opinion for i in allids(model))

    genocide!(empty_model)
    AgentsIO.populate_from_csv!(empty_model, "test.csv", Models.HKAgent, Dict(:op1 => 1, :op2 => 1); row_number_is_id = true)

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].old_opinion == empty_model[i].old_opinion for i in allids(model))
    @test all(model[i].new_opinion == empty_model[i].new_opinion for i in allids(model))
    @test all(model[i].previous_opinion == empty_model[i].previous_opinion for i in allids(model))

    AgentsIO.dump_to_csv("test.csv", allagents(model), [:id, :old_opinion])

    open("test.csv", "r") do f
        @test length(split(readline(f), ',')) == 2
    end

    genocide!(empty_model)
    AgentsIO.populate_from_csv!(empty_model, "test.csv")

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].old_opinion == empty_model[i].old_opinion for i in allids(model))
    @test all(model[i].new_opinion == empty_model[i].new_opinion for i in allids(model))
    @test all(model[i].previous_opinion == empty_model[i].previous_opinion for i in allids(model))

    model, _ = Models.battle(; fighters = 10)
    empty_model, _ = Models.battle(; fighters = 0)

    AgentsIO.dump_to_csv("test.csv", allagents(model))
    AgentsIO.populate_from_csv!(empty_model, "test.csv")

    @test nagents(empty_model) == nagents(model)
    @test all(haskey(empty_model.agents, i) for i in allids(model))
    @test all(model[i].pos == empty_model[i].pos for i in allids(model))
    @test all(model[i].has_prisoner == empty_model[i].has_prisoner for i in allids(model))
    @test all(model[i].capture_time == empty_model[i].capture_time for i in allids(model))
    @test all(model[i].shape == empty_model[i].shape for i in allids(model))

    rm("test.csv")
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 3b47d4a0-f7fa-11eb-0798-d3633587bd07
md"# Empezando

Temas:

1. Cómo imprimir
2. Cómo asignar variables
3. Cómo comentar
4. Sintaxis para matemáticas básicas"

# ╔═╡ 5a6a11d5-651b-4c26-accf-50a2387a13c8
md"## Cómo imprimir

En Julia, comúnmente se usa `println()` para imprimir, pero `print()` también funciona. ¿Cuál es la diferencia?"

# ╔═╡ 86f39c0a-5521-4cc3-aa4e-d7a3651d8022
println("¡Hola, CORE 2021!")

# ╔═╡ 77ef2843-9190-4767-9031-040bf71b9c6b
print("¡Hola, CORE 2021!")

# ╔═╡ 8186de26-f7be-479d-9af5-6c213e0338dc
md"## Cómo asignar variables

Solo se necesita un nombre de una variable, un valor y un signo de igual.

Julia se encarga de los _tipos_ por nosotros."

# ╔═╡ 104f1b92-873b-453c-93d1-21e7c62dfcb0
mi_respuesta = 42

# ╔═╡ 4bce4673-a035-47eb-9518-3002b7815a83
typeof(mi_respuesta)

# ╔═╡ c141093f-65b7-4eec-beeb-da2de8b1c24e
mi_pi = 3.14159

# ╔═╡ bef34fa9-72ae-45a1-b53f-72dacc493044
typeof(mi_pi)

# ╔═╡ ed58ead3-b719-4bb9-a8ac-8677b9b86404
md"Después de asignar una variable, se puede reasignar un valor de un tipo diferente a esa variable sin ningún problema (porque Julia es un lenguaje dinámico), pero no en Pluto. Usa el REPL para estas líneas."

# ╔═╡ e116f6b8-7ce7-430f-9fca-8182f76496cd
md"## Cómo comentar"

# ╔═╡ ef1b1f95-237f-4122-961a-0ff25114ed37
# Se pueden dejar comentarios de una sola línea usando el signo de "gato"

# ╔═╡ 33386794-6953-4d59-9906-a69d690be41d
#=

Para comentarios multilínea,
se utiliza la secuencia " #= ... =#".

=#

# ╔═╡ cabd36ab-10c2-4d47-b7aa-eeddfb7ebf73
md"## Sintaxis para matemáticas básicas"

# ╔═╡ e7df69fc-635a-4508-8bbb-b8fd67c938c5
suma = 3 + 7

# ╔═╡ 26844704-c42a-4f7f-b129-19cc012c13fb
diferencia = 10 - 3

# ╔═╡ c618a243-9494-4963-acfe-e9393ec30f45
producto = 20 * 5

# ╔═╡ 7e652f15-a2bd-4aec-b5cc-abcc188a6b78
cociente = 100 / 10

# ╔═╡ 3dde5f23-383e-482c-a703-7ba655328c36
potencia = 10 ^2

# ╔═╡ 074a1980-5ee0-4af1-a357-38ba0ba43bed
módulo = 101 % 2

# ╔═╡ 642f5b99-fd8e-4a55-b9b1-580cba30653c
md"## Ejercicios

### 1.1

Busca la documentación para la función `convert`."

# ╔═╡ 7686b892-9cb6-44ba-9f10-20a3fa984793
md"### 1.2

Asigna `365` a una variable llamada `días`."

# ╔═╡ 52ba4bde-9043-43c4-8445-d422d073dfc6
# Escribe tu respuesta

# ╔═╡ 9ea03840-a5bf-4dfd-b6e5-3a8e2f3d8ad4
if !@isdefined días
	md"""
	!!! danger "Incorrecto"
		La variable `días` no ha sido definida.
	"""
else
	if !(días isa Int && días == 365)
		md"""
		!!! warning "Incorrecto"
			La variable ya existe, pero no le asignaste `365`. Recuerda que debe ser un número entero.
		"""
	else
		md"""
		!!! correct "Correcto"
			¡Bien hecho!
		"""
	end
end

# ╔═╡ 2108d797-5d83-43ee-9fbc-0d2ad16c5aed
md"### 1.3
Convierte `días` a un flotante y asígnala a la variable `días_flotante`."

# ╔═╡ 7d29b778-cc6b-446a-afe8-018265f49785
# Escribe tu respuesta

# ╔═╡ 91219fe7-8f18-43fd-88cf-051b97364172
if !@isdefined días_flotante
	md"""
	!!! danger "Incorrecto"
		La variable `días_flotante` no ha sido definida.
	"""
else
	if !(días_flotante isa Float64 && días_flotante == 365.0)
		md"""
		!!! warning "Ya casi"
			La variable ya existe, pero no es un número de tipo flotante o no tiene el valor solicitado.
		"""
	else
		md"""
		!!! correct "Correcto"
			¡Bien hecho!
		"""
	end
end

# ╔═╡ 77cc11b3-e1b8-477e-8e4f-cc20230f3ebd
md"""### 1.4

Observa qué sucede cuando se ejecuta lo siguiente:

```julia
convert(Int64, "1")
```
y

```julia
parse(Int64, "1")
```
"""

# ╔═╡ e3964dc0-8979-4748-8e31-bfc88a620578

# ╔═╡ 36c014a3-5558-4a75-93fc-2aefd4eb61ca

# ╔═╡ Cell order:
# ╟─3b47d4a0-f7fa-11eb-0798-d3633587bd07
# ╟─5a6a11d5-651b-4c26-accf-50a2387a13c8
# ╠═86f39c0a-5521-4cc3-aa4e-d7a3651d8022
# ╠═77ef2843-9190-4767-9031-040bf71b9c6b
# ╟─8186de26-f7be-479d-9af5-6c213e0338dc
# ╠═104f1b92-873b-453c-93d1-21e7c62dfcb0
# ╠═4bce4673-a035-47eb-9518-3002b7815a83
# ╠═c141093f-65b7-4eec-beeb-da2de8b1c24e
# ╠═bef34fa9-72ae-45a1-b53f-72dacc493044
# ╟─ed58ead3-b719-4bb9-a8ac-8677b9b86404
# ╟─e116f6b8-7ce7-430f-9fca-8182f76496cd
# ╠═ef1b1f95-237f-4122-961a-0ff25114ed37
# ╠═33386794-6953-4d59-9906-a69d690be41d
# ╟─cabd36ab-10c2-4d47-b7aa-eeddfb7ebf73
# ╠═e7df69fc-635a-4508-8bbb-b8fd67c938c5
# ╠═26844704-c42a-4f7f-b129-19cc012c13fb
# ╠═c618a243-9494-4963-acfe-e9393ec30f45
# ╠═7e652f15-a2bd-4aec-b5cc-abcc188a6b78
# ╠═3dde5f23-383e-482c-a703-7ba655328c36
# ╠═074a1980-5ee0-4af1-a357-38ba0ba43bed
# ╟─642f5b99-fd8e-4a55-b9b1-580cba30653c
# ╟─7686b892-9cb6-44ba-9f10-20a3fa984793
# ╠═52ba4bde-9043-43c4-8445-d422d073dfc6
# ╟─9ea03840-a5bf-4dfd-b6e5-3a8e2f3d8ad4
# ╟─2108d797-5d83-43ee-9fbc-0d2ad16c5aed
# ╠═7d29b778-cc6b-446a-afe8-018265f49785
# ╟─91219fe7-8f18-43fd-88cf-051b97364172
# ╟─77cc11b3-e1b8-477e-8e4f-cc20230f3ebd
# ╠═e3964dc0-8979-4748-8e31-bfc88a620578
# ╠═36c014a3-5558-4a75-93fc-2aefd4eb61ca


text:
#= original implementation from
Single-Particle Diffusion Characterization by Deep Learning
Biophys J. 2019 Jul 23;117(2):185-192. doi: 10.1016/j.bpj.2019.06.015. Epub 2019 Jun 22.
Naor Granik et. al.
https://github.com/AnomDiffDB/DB
=#

function mittag_leffer_rand(β::Real, N::Int, γ::Real)
    t = -log.(rand(N))
    u = rand(N)
    w = sin(β * π) ./ tan.(β * π .* u) .- cos(β * π)
    t .*= w ./ β
    t .*= γ
end

function symmetric_alpha_levy(α::Real, N::Int, γ::Real)
    u = rand(N)
    v = rand(N)
    ϕ = π .* (v .- 0.5)
    w = sin.(α .* ϕ) ./ cos.(ϕ)
    z = -log.(u) .* cos.(ϕ)
    z ./= cos.((1 - α) .* ϕ)
    γ .* w .* z .^ (1 - (1 / α))
end

find_nearest(array::AbstractArray, value) = argmin(abs.(array .- value))

"""
    ctrw(α, γ, T, N) -> Vector{Float}, Vector{Float}, Vector{Float64}

Return a simulated CTRW trajectory presented by Fulger et. al., Phys. Rev. E(2008).

**Arguments**

- `α`: exponent of the waiting time distribution function.
- `β`: scale parameter for the mittag-leffler and alpha stable distributions.
- `N`: number of points to generate.
- `T`: end time.
"""
function ctrw(α::Real, γ::Real, T::Int, N::Int)
    jumpsX = mittag_leffer_rand(α, N, γ)
    rawTimeX = cumsum(jumpsX)
    tX = rawTimeX * T / maximum(rawTimeX)
    jumpsY = mittag_leffer_rand(α, N, γ)
    rawTimeY = cumsum(jumpsY)

    tY = rawTimeY * T / maximum(rawTimeY)
    x = cumsum(symmetric_alpha_levy(α, N, γ^(α / 2)))
    y = cumsum(symmetric_alpha_levy(α, N, γ^(α / 2)))
    tOut = collect(0:1:N) .* T / N
    xOut = Vector{Float64}(undef, N)
    yOut = Vector{Float64}(undef, N)

    for n = 1:N
        xOut[n] = x[find_nearest(tX, tOut[n]), 1]
        yOut[n] = y[find_nearest(tY, tOut[n]), 1]
    end
    return xOut, yOut, tOut
end

"""
Permute the states of `hmm` according to `perm`.

**Arguments**

- `perm::Vector{<:Integer}`: permutation of the states.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.8 0.2; 0.1 0.9], [Normal(0,1), Normal(10,1)])
hmm = permute(hmm, [2, 1])
hmm.A # [0.9 0.1; 0.2 0.8]
hmm.B # [Normal(10,1), Normal(0,1)]
```
"""

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
# journalentry.jl
# ===============
#
# A type and method used for constructing journaled sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

const SRC = UInt8
const SRC_NULL = 0b0000
const SRC_REFERENCE = 0b0001
const SRC_INSERTION = 0b0010

"""
A JournalEntry is a struct of values which are used to represent some sub-range
of either the reference or the insertion buffer in a Journal.
"""
immutable JournalEntry
    "A flag for where the segment maps (reference or insertion buffer)."
    source::SRC

    "Position in the reference string, or insertion buffer."
    physicalpos::Int

    "The physical position in the virtual string."
    virtualpos::Int

    "The physical position in the reference string."
    physicaloriginpos::Int

    "The length of the segment."
    len::Int
end

const JE = JournalEntry

JournalEntry() = JournalEntry(SRC_NULL, 0, 0, 0, 0)

source(entry::JE) = entry.source

"Test if the segment represents a region of the reference."
isorigional(entry::JE) = source(entry) == SRC_REFERENCE

"Test if the segment represents a region of the insertion buffer."
isinsertion(entry::JE) = source(entry) == SRC_INSERTION

"Get the physical location in the target."
vposition(entry::JE) = entry.virtualpos

"Get the position in the origonal string, or insertion buffer."
pposition(entry::JE) = entry.physicalpos

"Get the physical origin position in the reference string."
poposition(entry::JE) = entry.phyoriginpos

ltvirtual(entrya::JE, entryb::JE) = isless(vposition(entrya), vposition(entryb))

Base.isless(entrya::JE, entryb::JE) = ltvirtual(entrya, entryb)

Base.length(entry::JE) = entry.len

Base.isnull(entry::JE) = source(entry) == SRC_NULL

function Base.show(io::IO, ::MIME"text/plain", entry::JE)
    println(io, "Source of Segment: $(source(entry))")
    println(io, "Virtual Position: $(vposition(entry))")
    println(io, "Physical Position: $(pposition(entry))")
    println(io, "Physical Origin Position: $(poposition(entry))")
    println(io, "Length of Segment: $(length(entry))")
end


text:
module Data

include("model.jl")
export SummaryPair
include("datasets.jl")
export xsum, cnn_dm
include("utils.jl")
export out_dir, snapshot_file, snapshot_files

function __init__()
    register(xsum_dependency)
    register(cnn_dm_dependency)
end

end # module


text:
function replace_tables(s::String)
    lines = split(s,"\r\n")
    out=String[]
    intable=false
    for line in lines

        if startswith(line, r"\s*\|\|") && !intable  # first row of table
            println("found table")
            push!(out,"")   # empty line needed
            push!(out,replace(line, "||" => "|" ))
            push!(out, rstrip(replace(line, r"\s*(\|\|[^\|]*)" => "| --- ") , ['-', ' ']))
            intable=true
        elseif startswith(line, r"\s*\|\|") && intable # non first row of table
            push!(out,replace(line, "||" => "|"))
        else # normal text
            push!(out,line)
            intable=false
        end
    end
    return join(out,"\r\n")
end


text:
using Pkg
 if pwd() == pkgdir(RvSpectML) Pkg.activate("examples") end

verbose = true
 if verbose   println("# Loading RvSpecML")    end
 using RvSpectML
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_expres_data_101501.jl"))
continuum_normalize_spectra!(all_spectra)

order_list_timeseries = extract_orders(all_spectra,pipeline_plan)

linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",linelist_for_ccf_filename)
line_list_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true)

mask_scale_factors = [  1.0, 2, 4.0, 6.0, 8, 10, 12 ]
 rms_rvs = zeros(length(mask_scale_factors))
 rms_nightly_rvs = zeros(length(mask_scale_factors))
 rms_within_night_rvs = zeros(length(mask_scale_factors))
 println("Starting Tophat CCFs w/o variances...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors)
   println("# mask_scale_factor = ", mask_scale_factor)
   local (ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_scale_factor=mask_scale_factor, mask_type=:tophat,
   ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=false, recalc=true)
   #local (ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_merged, pipeline_plan, mask_scale_factor=mask_scale_factor, mask_type=:tophat,
   #ccf_mid_velocity=ccf_mid_velocity, recalc=true)
   alg_fit_rv = EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=5)
   local rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
   rms_rvs[i] = std(rvs_ccf.-mean(rvs_ccf))
   rms_nightly_rvs[i] = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf.-mean(rvs_ccf)))
   rms_within_night_rvs[i] = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf.-mean(rvs_ccf))
 end

mask_scale_factors2 = [  1, 2, 4.0, 8.0, 10, 12] #, 16, 18, 20, 22, 24, 26, 28, 30, 32 ]
 rms_rvs2 = zeros(length(mask_scale_factors2))
 rms_nightly_rvs2 = zeros(length(mask_scale_factors2))
 rms_within_night_rvs2 = zeros(length(mask_scale_factors2))
 println("Starting Tophat CCFs w/ variances...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors2)
   println("# mask_scale_factor = ", mask_scale_factor)
   ((ccfs_expr, ccf_vars_expr), v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:tophat, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
   rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, ccf_vars_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
   rms_rvs2[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
   rms_nightly_rvs2[i] = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr)))
   rms_within_night_rvs2[i] = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr))
 end

mask_scale_factors3 = [  4.0, 8.0, 10, 12, 16, 18, 20] # , 22, 24 , 26, 28, 30, 32]
 rms_rvs3 = zeros(length(mask_scale_factors3))
 rms_nightly_rvs3 = zeros(length(mask_scale_factors3))
 rms_within_night_rvs3 = zeros(length(mask_scale_factors3))
  println("Starting Half Cos...")
  for (i,mask_scale_factor) in enumerate(mask_scale_factors3)
    println("# mask_scale_factor = ", mask_scale_factor)
    (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:halfcos, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
    rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
    #((ccfs_expr, ccf_vars_expr), v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:halfcos, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
    #rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, ccf_vars_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
    rms_rvs3[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
    rms_nightly_rvs3[i] = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr)))
    rms_within_night_rvs3[i] = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr))
  end

mask_scale_factors4 = [  4.0, 6.0, 8.0, 9, 10 ] # ,11, 12, 15, 18 ]
 rms_rvs4 = zeros(length(mask_scale_factors4))
 rms_nightly_rvs4 = zeros(length(mask_scale_factors4))
 rms_within_night_rvs4 = zeros(length(mask_scale_factors4))
   println("Starting Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors4)
      println("# mask_scale_factor = ", mask_scale_factor)
      #alg_fit_rv = EchelleCCFs.MeasureRvFromCCFQuadratic()  # Why isn't this doing better than Gaussian?
      alg_fit_rv = EchelleCCFs.MeasureRvFromCCFGaussian()
      #(ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:gaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
      #rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
      ((ccfs_expr, ccfs_var_expr), v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:gaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true , recalc=true)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, ccfs_var_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, alg_fit_rv=alg_fit_rv, recalc=true, bin_nightly=true)
      rms_rvs4[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_nightly_rvs4[i] = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr)))
      rms_within_night_rvs4[i] = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr))

    end

mask_scale_factors5 = [  #= 6.0,8.0, =# 10.0, 12.0, 14.0, 16.0, 18.0, 20.0]
 rms_rvs5 = zeros(length(mask_scale_factors5))
  rms_nightly_rvs5 = zeros(length(mask_scale_factors5))
  rms_within_night_rvs5 = zeros(length(mask_scale_factors5))

   println("Starting Super Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors5)
      println("# mask_scale_factor = ", mask_scale_factor)
      #alg_fit_rv = EchelleCCFs.MeasureRvFromCCFQuadratic()  # Why isn't this doing better than Gaussian?
      alg_fit_rv = EchelleCCFs.MeasureRvFromCCFGaussian()
      #(ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:supergaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
      #rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
      ((ccfs_expr, ccf_vars_expr), v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:supergaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, ccf_vars_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_nightly=true)
      rms_rvs5[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_nightly_rvs5[i] = std(bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr)))
      rms_within_night_rvs5[i] = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_expr.-mean(rvs_ccf_expr))

    end

if make_plot(pipeline_plan, :ccf_total)
  using Plots

  plot(mask_scale_factors, rms_rvs, color=1, label="Tophat")
    scatter!(mask_scale_factors, rms_nightly_rvs, color=1, label=:none)
  plot!(mask_scale_factors2, rms_rvs2, color=2, label="Tophat w/ CCF Var")
    scatter!(mask_scale_factors2, rms_nightly_rvs2, color=2, label=:none)
  plot!(mask_scale_factors3, rms_rvs3, color=3, label="Half Cos")
    scatter!(mask_scale_factors3, rms_nightly_rvs3, color=3, label=:none)
  plot!(mask_scale_factors4, rms_rvs4, color=4, label="Gaussian")
    scatter!(mask_scale_factors4, rms_nightly_rvs4, color=4, label=:none)
  plot!(mask_scale_factors5, rms_rvs5, color=5, label="Super-Gaussian")
    scatter!(mask_scale_factors5, rms_nightly_rvs5, color=5, label=:none)
  ylims!(2,8)
end


text:
struct vcovBoot <: vcov
    B::Int64
    cluster::Matrix
    n_cluster_vars::Int64
    β_0::Float64
    function vcovBoot(B::Int64, cluster::Matrix; β_0 = 0.0)
        n_cluster_vars = size(cluster, 2)
        if n_cluster_vars > 2
            error("Only support one bootstrapped cluster.")
        end
        new(B, cluster, n_cluster_vars, β_0)
    end
end

function rademacher(N)
    D = rand(N) .> 0.5
    rademacher = D .* 1 .+ (1 .- D) .* -1
    return rademacher
end

function rademacher(dims::Tuple)
   hcat(rademacher.(fill(dims[1], dims[2]))...)
end

function inference(N::Int,
                   K::Int,
                   resid::Vector,
                   β::Vector,
                   XX_inv::Matrix,
                   X::Matrix,
                   vcov::vcovBoot)
    if vcov.n_cluster_vars != 1
        error("Only support one clustered variable Bootstrap atm")
    end
    B = vcov.B
    β_0 = vcov.β_0
    # if length(β_0) != K && length(β_0) == 1
    #     β_0 = fill(β_0[1], K)
    # end
    cluster_vals = vcov.cluster[:, 1]
    unique_clusters = unique(cluster_vals)
    n_clusters = length(unique_clusters)

    numeric_index_ids = Dict(unique_clusters .=> 1:n_clusters)
    cluster_ids = getindex.(Ref(numeric_index_ids), cluster_vals)

    radm_matrix = rademacher((N, B))

    # Perform inference for realised data
    se_r, pval_r, σ_sq_r, vcov_matrix_r = inference(
        N,
        K,
        resid,
        β,
        XX_inv,
        X,
        vcovCluster(vcov.cluster)
    )
    t_r = (β .- β_0) ./ se_r
    B = vcov.B
    t_b = Matrix{Float64}(undef, (B, K))
    se_b_vec = similar(t_b)
    for i in 1:B
        resid_b = resid .* radm_matrix[cluster_ids, i]
        y_b = X * β .+ resid_b
        β_b = X \ y_b
        se_b, _, σ_sq_b, _ = inference(
            N,
            K,
            resid_b,
            β_b,
            XX_inv,
            X,
            vcovCluster(vcov.cluster)
        )
        t_b[i, :] = (β_b .- β) ./ se_b
        se_b_vec[i, :] = se_b
    end
    pval_lower = (1/B) .* sum(t_b' .> t_r, dims = 2)
    pval_upper = (1/B) .* sum(t_b' .< t_r, dims = 2)
    pval = 2 .*  minimum(hcat(pval_lower, pval_upper), dims = 2)
    return mean(se_b_vec, dims = 1)[:], pval[:], σ_sq_r, vcov_matrix_r
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
function getreferencegenomes(syn, referenceFolder::AbstractString, referenceNames::AbstractArray)
    ids = map( x->getchildbyname(syn, referenceFolder, string(x,".fasta")), referenceNames )
    localPaths = map(x->localpath(syn,x), ids)
    referenceGenomes = map(loadfasta, localPaths)

    if all(x->length(x)==1, referenceGenomes)
        referenceGenomes = map(x->x[1][2], referenceGenomes) # NB: Only for genomes with 1 segment, get the sequence part
    end

    ids, referenceGenomes
end


text:
using Coconet
#=
    The core idea here is that we can derive Tonnetz very simply,
    as a small repeating tile.

    We start with a small number game, where we start from 0,
    and increase by 4 when going to the right, and 1 when going up

     ⋮  ⋮  ⋮  ⋮
    4  8 12 16 ...
    3  7 11 15 ...
    2  6 10 14 ...
    1  5  9 13 ...
    0  4  8 12 ...

    in other words, the value n of a position (x, y) can be expressed as
        n = 4x + y

    You may also notice some repetition
        starting with the 4s and going to the right.
        In fact, we can tile this pattern infinitely.

         ⋮  ⋮  ⋮  ⋮
    ... 4  8 12 16 ...
    ... 3  7 11 15 ...
    ... 2  6 10 14 ...
    ... 1  5  9 13 ...
    ... 0  4  8 12 16 ...
       ... 3  7 11 15 ...
       ... 2  6 10 14 ...
       ... 1  5  9 13 ...
       ... 0  4  8 12 ...
          ⋮   ⋮  ⋮  ⋮

    Now, since its infinitely repeating, it's not possible to have an inverse,
        but we can fix that too.
    If we just restrict our space to the first 4 rows, we can densely describe
        the infinite space pretty simply, while working our way toward
        an invertible function.

    We can describe a point's position x, y as
        x = floor(n / 4) <-- we only want integer values
        y = n % 4

    This grid, is the only way we will be interacting with notes.
    Any changes to the space will be handled by multiplying our results by
        a matrix representing a linear transformation.

    The last thing we need to do is add a y_offset, so that we can
        reach any point in the grid, which gives

        x = floor(n / 4)
        y = (n % 4) + y_offset

    We now have a geometric interpretation of all possible notes.
        You could even use negative values if you have an idea of what a
        "negative pitch" could possibly mean.
=#

#= Convert a pitch into a position in the tonnetz grid

        Pitch is just the numeric pitch (0-127 inclusive)
        y_offset lets you control where in the infinite
            grid of possible spaces you want to be.
            I intend to use it to center on the last note played,
            thus sort of minimizing the distance between consecutive notes
            in tonnetz space
=#
function pitch_to_tonnetz(pitch::Integer; y_offset = 0::Integer)

    # Subtract the vertical offset so pitches aren't changed
    pitch = pitch - y_offset

    # Since columns increase by 4,
    #   the x component is just the pitch divided by 4 with no remainder, and
    #   the y component is just the remainder
    x = floor(pitch / 4)
    y = pitch % 4

    #= Julia allows negative solutions to the mod function,
        e.g. -1 % 4 = -1
        but we want -1 % 4 = 3,
        So we have to manually correct for that =#
    y = (y < 0) ? 4 + y : y

    #= All our numbers are now correct,
    But our space is still anchored to (0, 0),
    So add the offset to y to get the correct location =#
    y = y + y_offset

    # We now have an infinite, navigatable grid of pitches
    #   that are related spacially to the tonnetz
    return x, y
end

# Convert a position on the tonnetz grid back to a pitch
function tonnetz_to_pitch(x, y)

    # Choose it based on its grid position
    pitch = x * 4 + y

    # Ta-da
    return Integer(round(pitch))
end

# Convert a song of tonnetz into an array of pitches
function song_to_pitch_array(song::AbstractArray)
    pitch_size = size(song, 1)
    sequence_length = size(song, 2)
    num_instruments = size(song, 3)

    output = zeros(sequence_length, num_instruments)

    for i in 1:num_instruments
        for j in 1:sequence_length
            pitch = NaN
            if round(song[3, j, i]) > 0
                pitch = Coconet.tonnetz_to_pitch(song[1:2, j, i])
            end

            output[i, j] = pitch
        end
    end
end

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    Aitchison

Aitchison distance between compositions.
"""
struct Aitchison <: Metric end

(::Aitchison)(c₁::Composition, c₂::Composition) = distance(c₁, c₂)

(d::Aitchison)(w₁, w₂) = d(Composition(w₁), Composition(w₂))

result_type(::Aitchison, w₁, w₂) = Float64


text:
"""
    fuzzy_scores(predicted, real; neighborhood, ignored=())

Return a fuzzy score for each cell. Compares each cell of the two arrays.
It finds the shortest distance within the neighborhood where the
value in the `real` array is present in `predict` array.

The longer the distance, the lower score for this cell. Maximum
is 1 and if not present its 0. The weight function is linear.

- `predict`: Array of predicted data
- `real`: Array of real target data

# Keywords

- `neighborhood`: any `Neighborhood` oject, `Window(1)` by default.
- `banlist`: Values in the arrays which will not be evaluated.

Adapted from: https://github.com/johanlahti/urban-lu-model/blob/master/script/Utils.py
"""
function fuzzy_scores(predicted, real; neighborhood=Window{1}(), ignored=0)
    broadcast(real, predicted, CartesianIndices(real)) do r, p, I
        r in ignored && return 0.0

        weight = 0.0
        shortestdist = Inf

        # Get neighborhood values from the `real` values array
        real_hood = _setbuffer(neighborhood, _getwindow(real, Tuple(I)))
        # Search the neighborhood
        for (hood_val, dist) in zip(neighbors(real_hood), distances(real_hood))
            if p === hood_val && dist < shortestdist
                # Store the shortest distance at which we found the value
                shortestdist = dist
                weight = 1 - dist / (length(hood) + 1)
            end
        end
        return weight
    end
end


text:
using SafeTestsets, Test

@safetestset "Differentiation Test" begin include("diff.jl") end
@safetestset "Array Test" begin include("array.jl") end


text:
import HTTP

function isdate(t::AbstractVector{T})::Bool where {T<:TimeType}
    h = Dates.hour.(t)
    m = Dates.minute.(t)
    s = Dates.second.(t)
    ms = Dates.millisecond.(t)
    return all(h.==h[1]) && all(m.==m[1]) && all(s.==s[1]) && all(ms.==ms[1])
end

function csvresp(resp::HTTP.Response; sort::Char='d')
    @assert resp.status == 200 "Error in download request."
    rowdata = Vector{String}(split(String(resp.body), '\n'))
    header = Vector{String}(split(popfirst!(rowdata), ','))
    pop!(rowdata)
    if sort == 'd'
        reverse!(rowdata)
    elseif sort != 'a'
        error("Argument `sort` must be either 'a' (ascending) or 'd' (descending)")
    end
    N = length(rowdata)
    k = length(header)
    v = map(s -> Array{String}(split(s, ',')), rowdata)
    source_is_google = (header[1] == "\ufeffDate")
    if source_is_google
        header[1] .= "Date"
        format = Dates.DateFormat("dd-uuu-yy")
        t = map(s -> Dates.DateTime(s[1], format), v)
        t .+= Dates.Year(2000)  # instead of year 0017, take year 2017
    else
        format = Dates.DateFormat("yyyy-mm-dd")
        t = map(s -> Dates.DateTime(s[1]), v)
    end
    isdate(t) ? t = Date.(t) : nothing
    data = zeros(Float64, (N,k-1))
    if length(header) == 2 && header[2] == "Stock Splits"
        # Logic to be applied for stock splits for Yahoo Finance downloads
        @inbounds for i in 1:N
            stock_split_string = split(v[i][2], '/')
            split_a = float(stock_split_string[1])
            split_b = float(stock_split_string[2])
            is_rev_split = split_a < split_b
            data[i,1] .= split_b / split_a
        end
    else
        # Standard logic
        @inbounds for i in 1:N
            @inbounds for j in 2:k
                if v[i][j] == "" || v[i][j] == "null"
                    v[i][j] = "NaN"
                end
                data[i,j-1] = parse(Float64, v[i][j])
            end
        end
    end
    return (data, t, header)
end


text:
"""
    MvCategorical

Product of independent `Categorical` distributions.
"""
const MvCategorical = ProductDistribution{Categorical}

"""
    MvCategoricalPrior

Conjugate multivariate Dirichlet prior over the parameters of a `MvCategorical` distribution.
"""
const MvCategoricalPrior = ProductDistribution{Dirichlet}

# ## MLE

# function Distributions.fit_mle(
#     ::Type{<:MvCategorical},
#     ks::AbstractVector{Integer},
#     x::AbstractMatrix{Integer},
# )
#     marginals = [fit_mle(Categorical, ks[i], x[i, :]) for i = 1:length(ks)]
#     p = [probs(dist) for dist in marginals]
#     return MvCategorical(p)
# end

# function Distributions.fit_mle(
#     ::Type{<:MvCategorical},
#     ks::AbstractVector{Integer},
#     x::AbstractMatrix{Integer},
#     w::AbstractVector{Real},
# )
#     marginals = [fit_mle(Categorical, ks[i], x[i, :], w) for i = 1:length(ks)]
#     p = [probs(dist) for dist in marginals]
#     return MvCategorical(p)
# end

# function Distributions.fit_mle(::Type{<:MvCategorical}, x::AbstractMatrix{Integer})
#     ks = vec(maximum(x, dims = 2))
#     return fit_mle(MvCategorical, ks, x)
# end

# function Distributions.fit_mle(
#     ::Type{<:MvCategorical},
#     x::AbstractMatrix{Integer},
#     w::AbstractVector{Real},
# )
#     ks = vec(maximum(x, dims = 2))
#     return fit_mle(MvCategorical, ks, x, w)
# end

# ## MAP estimation

# function fit_map(
#     ::Type{<:MvCategorical},
#     prior::MvCategoricalPrior,
#     x::AbstractMatrix{Integer},
# )
#     marginals = [
#         fit_map(Categorical, CategoricalPrior(prior.α[i]), x[i, :]) for
#         i = 1:length(prior.α)
#     ]
#     p = [probs(dist) for dist in marginals]
#     return MvCategorical(p)
# end

# function fit_map(
#     ::Type{<:MvCategorical},
#     prior::MvCategoricalPrior,
#     x::AbstractMatrix{Integer},
#     w::AbstractVector{Real},
# )
#     marginals = [
#         fit_map(Categorical, CategoricalPrior(prior.α[i]), x[i, :], w) for
#         i = 1:length(prior.α)
#     ]
#     p = [probs(dist) for dist in marginals]
#     return MvCategorical(p)
# end


text:
export OvCF, OvController, OvMsg

import Optim

struct OvMsg{X,Z}
    δx::Timed{X}
    z::Timed{Optional{Z}}
end

to_optional(x::Timed{X}) where X = convert(Timed{Optional{X}}, x)

"""
Matrices are indexed by [time, id].
"""
struct OvLog{X,Z,U}
    ũ::Matrix{U}
    x̃::Matrix{X}
    z̃::Matrix{Z}
    δx::X
    z::Z
end

function to_array(m::MMatrix{n1,n2,T}) where {n1,n2,T}
    r = Matrix{T}(undef, n1, n2)
    copyto!(r, m)
end

"""
The OneVision Controller Framework.
"""
@kwdef struct OvCF{N,X,Z,U,S,H} <: ControllerFramework{X,Z,U,OvMsg{X,Z},OvLog{X,Z,U}}
    central::CentralControl{U,S}
    world_model::WorldDynamics{N}
    delay_model::DelayModel
    x_weights::SVector{N,X}
    u_weights::SVector{N,U}
    optim_options::Optim.Options
    save_log::FuncT{Tuple{ℕ,𝕋,X,Z},Bool} = FuncT(x -> false, Tuple{ℕ,𝕋,X,Z}, Bool)
end

function OvCF(central::CentralControl{U,S}, world_model, delay_model, x_weights, u_weights;
    X, Z, N, H,
    u_tol = 1e-4, loss_tol = 1e-6, loss_rel_tol = 1e-6, max_iters = 100,
    save_log = FuncT(x -> false, Tuple{ℕ,𝕋,X,Z}, Bool)
) where {U,S}
    OvCF{N,X,Z,U,S,H}(;central, world_model, delay_model, x_weights, u_weights, save_log,
        optim_options = Optim.Options(
            x_abstol = u_tol, f_abstol = loss_tol, f_reltol = loss_rel_tol,
            iterations = max_iters))
end

function OvCF(loss_model::RegretLossModel{N,X,U,S}, delay_model;
    Z, H,
    optim_options = Optim.Options(
        x_abstol=1e-4, f_abstol=1e-6, f_reltol=1e-6, iterations=100),
    save_log = FuncT(x -> false, Tuple{ℕ,𝕋,X,Z}, Bool)
) where {N,X,U,S}
    @unpack central, world_model, x_weights, u_weights = loss_model
    OvCF{N,X,Z,U,S,H}(;central, world_model, delay_model,
        x_weights, u_weights, save_log, optim_options)
end

@kwdef mutable struct OvController{N,X,Z,U,S,H,ΔT,Hf,WDy,Ctrl} <: Controller{X,Z,U,OvMsg{X,Z},OvLog{X,Z,U}}
    id::ℕ
    cf::OvCF{N,X,Z,U,S,H}
    τ0::𝕋
    τ::𝕋                                # -- At the start of each controll step --
    u_last::U                           # t = τ+Tu-1
    u_history::TimedQueue{U}            # t ∈ [τ-Tx,τ+Tu-1]
    "predicted x from last control step, used to measure disturbance."
    pred_x::Timed{X}                    # t = τ-Tx
    "self disturbance history"
    self_δx::TimedQueue{X}              # t ∈ [τ-Tx-Tc-ΔT-1, τ-Tx-2]
    self_z::TimedQueue{Optional{Z}}     # t ∈ [τ-Tx-Tc-ΔT, τ-Tx-1]
    ideal_xz::Timed{Each{Tuple{X,Z}}}   # t = τ-Tx-Tc-ΔT
    ideal_s::S                   # t = τ-Tx-Tc-ΔT
    fp_prob::ForwardPredictProblem{N,Hf,X,Z,U,WDy,Ctrl}
    plan_prob::TrajPlanningProblem{H,ΔT,X,U}
    logs::Dict{𝕋,OvLog{X,Z,U}}
end

OneVision.write_logs(ctrl::OvController) = ctrl.logs

"""
# Arguments
- `init_status`: The initial (x,z,u) of the agent. x and z are assumed to be
measured at `t0-Tx`.
"""
function OneVision.make_controllers(
    cf::OvCF{N,X,Z,U,S,H},
    init_status::Each{Tuple{X,Z,U}},
    τ::𝕋,
)::Tuple where {N,X,Z,U,S,H}
    dm = cf.delay_model
    @unpack Tx, Tu, Tc, Ta, ΔT = short_delay_names(dm)
    function mk_controller(id)
        x0, z0, u0 = init_status[id]
        u_history = TimedQueue(u0, τ-Tx, τ+Tu-1)
        x_dy = cf.world_model.dynamics[id]
        pred_x = Timed(τ-Tx, x0)
        self_δx = TimedQueue(zero(X), τ-Tx-Tc-ΔT-1, τ-Tx-2)
        self_z = TimedQueue(missing, τ-Tx-Tc-ΔT, τ-Tx-1, eltype=Optional{Z})
        ideal_xz = Timed(τ-Tx-Tc-ΔT, [(x, z) for (x, z, _) in init_status])
        ideal_s = init_state(cf.central, τ)

        fp_prob = let Hf = H * ΔT + Ta + ΔT
            ForwardPredictProblem(cf.world_model, cf.central; X, Z, Hf)
        end
        x_weights = cf.x_weights[id]
        u_weights = cf.u_weights[id]
        plan_prob = TrajPlanningProblem(
            H, ΔT, x_dy, x_weights, u_weights, optim_options = cf.optim_options,
        )

        OvController(;
            id, cf, τ0 = τ, τ, u_last = u0, u_history, pred_x,
            self_δx, self_z, ideal_xz, ideal_s,
            fp_prob, plan_prob, logs = Dict{𝕋,OvLog{X,Z,U}}(),
        )
    end

    function mk_q(_)
        receives = [[OvMsg(Timed(t-1, zero(x)), Timed{Optional{Z}}(t, missing))
                    for (x, z, u) in init_status
                    ] for t in range(τ-Tc-Tx, length = msg_queue_length(dm))]
        FixedQueue(receives)
    end

    ctrls::NTuple{N,OvController} = ntuple(mk_controller, N)
    msg_qs::NTuple{N,MsgQueue} = ntuple(mk_q, N)
    ctrls, msg_qs
end

function OneVision.control!(
    π::OvController{N,X,Z,U,S,H,ΔT},
    x::X,
    z::Z,
    msgs::Each{OvMsg{X,Z}},
)::Tuple{U,Each{OvMsg{X,Z}}} where {N,X,Z,U,S,H,ΔT}
    @unpack id, τ, τ0 = π
    world = π.cf.world_model
    fp_prob = π.fp_prob
    @unpack Tx, Tu, Tc, Ta = short_delay_names(π.cf.delay_model)
    zt = Timed(τ-Tx, z)
    xt = Timed(τ-Tx, x)

    δx = x - π.pred_x[τ-Tx] |> attime(τ-Tx-1)   # δx: t = τ-Tx-1
    pushpop!(π.self_δx, δx)     # π.self_δx: t ∈ [τ-Tx-1-Tc, τ-Tx-1]
    pushpop!(π.self_z, zt)       # π.self_z: t ∈ [τ-Tx-Tc, τ-Tx]
    new_msg = OvMsg(δx, to_optional(zt))
    should_log = π.cf.save_log((id, π.τ, x, zt.value))

    x_self::Timed{X} = isempty(π.u_history) ? xt : self_estimate(
        world.dynamics[id],
        xt,
        π.u_history.queue
    )[end]  # x_self: t = τ+Tu

    "the time offset within the ΔT period"
    t = mod1(τ - τ0, ΔT)

    # forward predict the ideal fleet trajectory from the state at t-Tx-Tc-1
    for j in 1:N
        if j == id continue end
        @unpack δx, z = msgs[j]
        fp_prob.δx[t, j] = δx[τ-Tx-Tc-1]
        fp_prob.z_obs[t, j] = z[τ-Tx-Tc]
    end

    if t == ΔT
        fp_prob.δx[1:length(π.self_δx),id] = collect(x.value for x in π.self_δx.queue)
        fp_prob.z_obs[1:length(π.self_z),id] = collect(x.value for x in π.self_z.queue)

        # x̃[t in τ-Tx-Tc-ΔT+1: τ+Tu+HΔT], ũ[t in τ-Tx-Tc-ΔT: τ+Tu+HΔT-1]
        ptime = τ-Tx-Tc-ΔT  # the starting time of forward_predict!
        @asserteq π.self_δx[1].time ptime
        ũ, x̃, z̃ = forward_predict!(
            fp_prob, π.ideal_xz[ptime], π.ideal_s, ptime, ΔT)
        π.ideal_xz = [(x̃[ΔT,j], z̃[ΔT,j]) for j in 1:N] |> attime(ptime+ΔT)  # π.ideal_xz: t = τ-Tx-Tc

        # local planning using path following

        # u_plan: [τ+Tu:τ+Tu+HΔT-1], x_plan: [τ+Tu+1:τ+Tu+HΔT]
        u_plan, loss = let
            n_x, n_u = length(X), length(U)
            x_path = SVector{H * ΔT}(x̃[1 + ΔT + Ta : end, id])  # [τ+Tu+1:τ+Tu+HΔT]
            u_path = SVector{H * ΔT}(ũ[1 + ΔT + Ta : end, id])  # [τ+Tu:τ+Tu+HΔT-1]
            plan_trajectory(π.plan_prob, x_self[τ+Tu], x_path, u_path, τ+Tu)
        end

        π.u_last = u_plan[1]
    end
    u = π.u_last

    π.pred_x = let
        x_dy = world.dynamics[id]
        u_old = pushpop!(π.u_history, Timed(τ+Tu, u))[τ-Tx]  # u_old at t = τ-Tx
        x1 = sys_forward(x_dy, x, u_old, τ-Tx)
        Timed(τ-Tx+1, x1)
    end

    if should_log
        π.logs[π.τ] = OvLog(to_array.((ũ, x̃, z̃))..., x, z)
    end

    π.τ += 1
    u, fill(new_msg, N)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using FactCheck

facts("Testing LSTM layer") do
    n_time_steps = 5
    net = Net(2; time_steps=n_time_steps)
    data, data_value = MemoryDataLayer(net, :data, (8,))
    rand!(data_value)

    lstm = LSTMLayer(:lstm, net, data, length(data))
    init(net)
    forward(net; phase=Latte.Test)

    W_ih = net.params[1].value
    W_ix = net.params[3].value

    W_C_simh = net.params[5].value
    W_C_simx = net.params[7].value

    W_fh = net.params[9].value
    W_fx = net.params[11].value

    W_oh = net.params[13].value
    W_ox = net.params[15].value

    V = net.params[17].value

    h = zeros(Float32, 8, 2)
    C = zeros(Float32, 8, 2)
    sigmoid(x)  = 1.0f0 ./ (1.0f0 .+ exp(-x))
    for t = 1:n_time_steps
        x = get_buffer(net, :datavalue, t)
        ih = W_ih' * h
        ix = W_ix' * x
        @fact ih --> roughly(get_buffer(net, symbol(net.ensembles[2].name, :value), t))
        @fact ix --> roughly(get_buffer(net, symbol(net.ensembles[3].name, :value), t))

        C_simh = W_C_simh' * h
        C_simx = W_C_simx' * x
        @fact C_simh --> roughly(get_buffer(net, symbol(net.ensembles[4].name, :value), t))
        @fact C_simx --> roughly(get_buffer(net, symbol(net.ensembles[5].name, :value), t))

        fh = W_fh' * h
        fx = W_fx' * x
        @fact fh --> roughly(get_buffer(net, symbol(net.ensembles[6].name, :value), t))
        @fact fx --> roughly(get_buffer(net, symbol(net.ensembles[7].name, :value), t))

        oh = W_oh' * h
        ox = W_ox' * x
        @fact oh --> roughly(get_buffer(net, symbol(net.ensembles[8].name, :value), t))
        @fact ox --> roughly(get_buffer(net, symbol(net.ensembles[9].name, :value), t))

        i = sigmoid(ih .+ ix)
        @fact i --> roughly(get_buffer(net, symbol(net.ensembles[11].name, :value), t))

        C_sim = tanh(C_simh .+ C_simx)
        @fact C_sim --> roughly(get_buffer(net, symbol(net.ensembles[13].name, :value), t))

        f = sigmoid(fh .+ fx)
        @fact f --> roughly(get_buffer(net, symbol(net.ensembles[15].name, :value), t))

        f_C = f .* C
        @fact f_C --> roughly(get_buffer(net, symbol(net.ensembles[16].name, :value), t))

        i_C = i .* C_sim
        @fact i_C --> roughly(get_buffer(net, symbol(net.ensembles[17].name, :value), t))

        C = i_C .+ f_C
        @fact C --> roughly(get_buffer(net, symbol(net.ensembles[18].name, :value), t))
        o = sigmoid(V' * C .+ oh .+ ox)
        h = o .* tanh(C)
        h_actual = get_buffer(net, symbol(net.ensembles[end].name, :value), t)
        @fact h_actual --> roughly(h)
    end
end

FactCheck.exitstatus()


text:
# Some of the code in this file is taken from PredictMD:
# https://github.com/bcbi/PredictMD.jl
# https://predictmd.net

import DataFrames

const RealOrMissing = Union{DataFrames.Missing, T} where T<:Real

function factor_columns(dat)
    p = size(dat, 2)
    is_factor = falses(p)
    for j = 1:p
        if !(eltype(dat[:, j]) <: RealOrMissing)
            is_factor[j] = true
        end
    end
    indcs = findall(is_factor)

    indcs
end

# @code_warntype factor_columns(d)

function factor_to_float(v::T) where T <: AbstractArray
    unique_cats = unique(v)         # unique categories
    sort!(unique_cats)
    cat_dictionary = Dict{String, Float64}()
    val = 1.0
    for k in unique_cats
        cat_dictionary[k] = val
        val += 1.0
    end
    n = length(v)
    res = zeros(n)
    for i = 1:n
        res[i] = cat_dictionary[v[i]]
    end
    res
end

function float_to_factor(v::T, levels::S) where T <: AbstractArray where S <: AbstractVector
    if levels == DataFrames.DataFrame()
        levels = zeros(1,1)
    end
    sort!(levels)
    str_vect = map(x -> levels[convert(Int, x)], v)
    result = Array(str_vect)
    return result
end

# This function behaves a bit like R's scale()
# function when it's called with MARGIN = 2.
function rscale(X, center, scale)
    n, p = size(X)
    res = zeros(n, p)
    for i = 1:n
        for j = 1:p
            res[i, j] = (X[i, j] - center[j])/scale[j]
        end
    end
    res
end

function column_ranges(X::T) where T <: AbstractMatrix
    p = size(X, 2)
    ranges = zeros(p)

    for j = 1:p
	ranges[j] = maximum(X[:, j]) - minimum(X[:, j])
    end
    ranges
end

# The `calculate_smote_pct_under` function is taken from PredictMD:
# https://github.com/bcbi/PredictMD.jl
# https://predictmd.net

function calculate_smote_pct_under(
        ;
        pct_over::Real = 0,
        minority_to_majority_ratio::Real = 0,
        )
    if pct_over < 0
        error("pct_over must be >=0")
    end
    if minority_to_majority_ratio <= 0
        error("minority_to_majority_ratio must be >0")
    end
    result = 100*minority_to_majority_ratio*(100+pct_over)/pct_over
    return result
end


text:
## Char8 ##

bitstype 8 Char8 <: ASCIIChar

## Conversions ##

char8(x) = convert(Char8,x)
char8(x::FloatingPoint) = char8(iround(x))

char8(s::ASCIIString) = reinterpret(Char8,copy(s.data))

char8(x::AbstractArray{Char8}) = x
import Base.iround_to
char8(x::AbstractArray)= iround_to(similar(x,Char8), x)

macro c8_str(s);   ex = interp_parse_bytes(s); :(reinterpret(Char8,($ex).data)); end