text:
using CSyntax
using Quaternions
using GLTF

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

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

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

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# we use glTF file format because it's OpenGL friendly and we don't need to care about
# how to parse and extract mesh data from other mesh file format like .obj format.
# load glTF file
sphere = GLTF.load(joinpath(@__DIR__, "sphere.gltf"))
sphere_data = [read(joinpath(@__DIR__, b.uri)) for b in sphere.buffers]

# load sphere position metadata
search_name(x, keyword) = x[findfirst(x->occursin(keyword, x.name), x)]
pos = search_name(sphere.accessors, "position")
pos_bv = sphere.bufferViews[pos.bufferView]
# load sphere index metadata
indices = search_name(sphere.accessors, "indices")
idx_bv = sphere.bufferViews[indices.bufferView]

# create VBO and EBO
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)
glBindBuffer(pos_bv.target, vbo)
glBufferData(pos_bv.target, pos_bv.byteLength, C_NULL, GL_STATIC_DRAW)
pos_data = sphere_data[pos_bv.buffer]
@c glBufferSubData(pos_bv.target, 0, pos_bv.byteLength, &pos_data[pos_bv.byteOffset])

ebo = GLuint(0)
@c glGenBuffers(1, &ebo)
glBindBuffer(idx_bv.target, ebo)
glBufferData(idx_bv.target, idx_bv.byteLength, C_NULL, GL_STATIC_DRAW)
idx_data = sphere_data[idx_bv.buffer]
@c glBufferSubData(idx_bv.target, 0, idx_bv.byteLength, &idx_data[idx_bv.byteOffset])

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(pos_bv.target, vbo)
glVertexAttribPointer(0, 3, pos.componentType, pos.normalized, pos_bv.byteStride, Ptr{Cvoid}(pos.byteOffset))
glEnableVertexAttribArray(0)

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

# ray casting
function screen2world(window::GLFW.Window, mouse_x::Real, mouse_y::Real)
    # screen space (window space)
    x = 2.0 * mouse_x / width - 1.0
    y = 1.0 - (2.0 * mouse_y) / height
    z = 1.0
    # normalised device space [-1:1, -1:1, -1:1]
    ray_nds = GLfloat[x, y, z]
    # clip space [-1:1, -1:1, -1:1, -1:1]
    ray_clip = GLfloat[ray_nds[1], ray_nds[2], -1, 1]
    # eye space [-x:x, -y:y, -z:z, -w:w]
    ray_eye = inv(get_projective_matrix(window, camera)) * ray_clip
    ray_eye = GLfloat[ray_eye[1], ray_eye[2], -1, 0]
    # world space [-x:x, -y:y, -z:z, -w:w]
    ray_world_homogenous = inv(get_view_matrix(camera)) * ray_eye
    ray_world = ray_world_homogenous[1:3]
    # normalize
    normalize!(ray_world)
    return ray_world
end

function raysphere(ray_origin_world, ray_direction_world, sphere_center_world, sphere_radius)
    intersection_dist = 0
    # quadratic parameters
    dist = ray_origin_world - sphere_center_world
    b = ray_direction_world ⋅ dist
    c = dist ⋅ dist - sphere_radius * sphere_radius
    Δ = b * b - c
    # no intersection
    Δ < 0 && (flag = false;)
    # one intersection (tangent ray)
    if Δ == 0
        # if behind viewer, throw away
        t = -b + √(Δ)
        t < 0 && (flag = false;)
        intersection_dist = t
        flag = true
    end
    # two intersections (secant ray)
    if Δ > 0
        t₁ = -b + √(Δ)
        t₂ = -b - √(Δ)
        intersection_dist = t₂
        # if behind viewer, throw away
        t₁ < 0 && t₂ < 0 && (flag = false;)
        t₁ ≥ 0 && t₂ < 0 && (intersection_dist = t₁;)
        flag = true
    end
    return flag, intersection_dist
end

# spheres in world
const NUM_SPHERES = 4
sphere_world = GLfloat[-2.0 0.0  0.0;
                        2.0 0.0  0.0;
                       -2.0 0.0 -2.0;
                        1.5 1.0 -1.0]
model_mats = map(1:4) do i
    GLfloat[ 1.0 0.0 0.0 sphere_world[i,1];
             0.0 1.0 0.0 sphere_world[i,2];
             0.0 0.0 1.0 sphere_world[i,3];
             0.0 0.0 0.0       1.0        ]
end

const SPHERE_RADIUS = 1
selected = -1
function mouse_click_callback(window::GLFW.Window, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint)
    if GLFW.PRESS == action
        xpos, ypos = GLFW.GetCursorPos(window)
        ray_world = screen2world(window, xpos, ypos)
        # ray sphere
        clicked = -1
        intersect = 0
        for i = 1:NUM_SPHERES
            sp = collect(sphere_world[i,:])
            flag, dist = raysphere(vec(camera.position), ray_world, sp, SPHERE_RADIUS)
            flag || continue
            if (clicked == -1) || (dist < intersect)
                clicked = i
                intersect = dist
            end
        end
        global selected = clicked
        println("sphere ", selected, " was clicked")
    end
    return nothing
end
# set mouse click callback
GLFW.SetMouseButtonCallback(window, mouse_click_callback)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "raypick.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "raypick.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)

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

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

let
updatefps = FPSCounter()
# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glBindBuffer(idx_bv.target, ebo)
    for i = 1:NUM_SPHERES
        if i == selected
            glUniform1f(blue_loc, 1.0)
        else
            glUniform1f(blue_loc, 0.0)
        end
        glUniformMatrix4fv(model_loc, 1, GL_FALSE, model_mats[i])
        glDrawElements(GL_TRIANGLES, indices.count, indices.componentType, Ptr{Cvoid}(0))
    end
    # check and call events
    GLFW.PollEvents()
    # move camera
    updatecamera!(window, camera)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module FieldValues

using ...VersionWeights

export FieldValue, Field, validmax, secondmax

# FieldValue is a hierarchical numeric type with 5 levels.
# When summing two FieldValues, the levels are summed independently.
# When comparing them, lower levels take precedence.
# The levels are used as such:
#  l0 : for hard constraints (dependencies and requirements)
#  l1 : for favoring higher versions of the explicitly required
#       packages
#  l2 : for favoring higher versions of all other packages
#  l3 : for favoring uninstallation of non-needed packages
#
struct FieldValue
    l0::Int64
    l1::VersionWeight
    l2::VersionWeight
    l3::Int64
    FieldValue(l0::Integer = 0,
               l1::VersionWeight = zero(VersionWeight),
               l2::VersionWeight = zero(VersionWeight),
               l3::Integer = 0) = new(l0, l1, l2, l3)
end

# This isn't nice, but it's for debugging only anyway
function Base.show(io::IO, a::FieldValue)
    print(io, a.l0)
    a == FieldValue(a.l0) && return
    print(io, ".", a.l1)
    a == FieldValue(a.l0, a.l1) && return
    print(io, ".", a.l2)
    a == FieldValue(a.l0, a.l1, a.l2) && return
    print(io, ".", a.l3)
    return
end

const Field = Vector{FieldValue}

Base.zero(::Type{FieldValue}) = FieldValue()

Base.typemin(::Type{FieldValue}) = (x=typemin(Int64); y=typemin(VersionWeight); FieldValue(x, y, y, x))

Base.:-(a::FieldValue, b::FieldValue) = FieldValue(a.l0-b.l0, a.l1-b.l1, a.l2-b.l2, a.l3-b.l3)
Base.:+(a::FieldValue, b::FieldValue) = FieldValue(a.l0+b.l0, a.l1+b.l1, a.l2+b.l2, a.l3+b.l3)

function Base.isless(a::FieldValue, b::FieldValue)
    a.l0 < b.l0 && return true
    a.l0 > b.l0 && return false
    c = cmp(a.l1, b.l1)
    c < 0 && return true
    c > 0 && return false
    c = cmp(a.l2, b.l2)
    c < 0 && return true
    c > 0 && return false
    a.l3 < b.l3 && return true
    return false
end

Base.:(==)(a::FieldValue, b::FieldValue) =
    a.l0 == b.l0 && a.l1 == b.l1 && a.l2 == b.l2 && a.l3 == b.l3

Base.abs(a::FieldValue) = FieldValue(abs(a.l0), abs(a.l1), abs(a.l2), abs(a.l3))

Base.copy(a::FieldValue) = FieldValue(a.l0, copy(a.l1), copy(a.l2), a.l3)

# if the maximum field has l0 < 0, it means that
# some hard constraint is being violated
validmax(a::FieldValue) = a.l0 >= 0

# like usual indmax, but favors the highest indices
# in case of a tie
function Base.indmax(f::Field)
    m = typemin(FieldValue)
    mi = 0
    for j = length(f):-1:1
        if f[j] > m
            m = f[j]
            mi = j
        end
    end
    @assert mi != 0
    return mi
end

# secondmax returns the (normalized) value of the second maximum in a
# field. It's used to determine the most polarized field.
function secondmax(f::Field, msk::BitVector = trues(length(f)))
    m = typemin(FieldValue)
    m2 = typemin(FieldValue)
    for i = 1:length(f)
        msk[i] || continue
        a = f[i]
        if a > m
            m2 = m
            m = a
        elseif a > m2
            m2 = a
        end
    end
    return m2 - m
end

# Support broadcasting like a scalar by default
Base.Broadcast.broadcastable(a::FieldValue) = Ref(a)

end

ext: jl

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using RecipesBase
using AstroImages: pix2world_xformatter, pix2world_yformatter

@testset "Plot recipes" begin
    data = randn(10, 10)
    img = AstroImage(data)
    wcs1 = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"])
    wcs2 = WCSTransform(2; ctype = ["GLON--", "GLAT--"])
    wcs3 = WCSTransform(2; ctype = ["TLON--", "TLAT--"])
    wcs4 = WCSTransform(2; ctype = ["UNK---", "UNK---"])
    wcs5 = WCSTransform(2;)

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, 1)
    @test getfield(rec[1], 1) == Dict{Symbol, Any}(:seriestype   => :heatmap,
                                                   :aspect_ratio => :equal,
                                                   :color        => :grays)

                                                   rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img)
    @test getfield(rec[1], 1) == Dict{Symbol, Any}(:seriestype   => :heatmap,
                                                   :aspect_ratio => :equal,
                                                   :color        => :grays)
    @test rec[1].args == (img.data[1],)

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, wcs1)
    @test rec[1].plotattributes[:xlabel] == "Right Ascension (ICRS)" && rec[1].plotattributes[:ylabel] == "Declination (ICRS)"

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, wcs2)
    @test rec[1].plotattributes[:xlabel] == "Galactic Coordinate" && rec[1].plotattributes[:ylabel] == "Galactic Coordinate"

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, wcs3)
    @test rec[1].plotattributes[:xlabel] == "ITRS" && rec[1].plotattributes[:ylabel] == "ITRS"

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, wcs4)
    @test rec[1].plotattributes[:xlabel] == "UNK---" && rec[1].plotattributes[:ylabel] == "UNK---"

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), img, wcs5)
    @test rec[1].plotattributes[:xlabel] == "" && rec[1].plotattributes[:ylabel] == ""

end

@testset "formatters" begin
    wcs1 = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"])
    wcs2 = WCSTransform(2; ctype = ["GLON--", "GLAT--"])

    @test pix2world_xformatter(255, wcs1) == "165.58°"
    @test pix2world_xformatter(255, wcs2) == 255.0

    @test pix2world_yformatter(255, wcs1) == "13.98°"
    @test pix2world_yformatter(255, wcs2) == 255.0
end


text:
export perlineiterator

function eachstep(x::StepValues{Nstep}, order=ntuple(i->i,Nstep)) where Nstep
  length(order) != Nstep && throw(ArgumentError())
  if order==(1,2,3)
    return ((i,j,k) for i=1:length(x.values[1]),
                        j=1:length(x.values[2]),
                        k=1:length(x.values[3]))
  elseif order==(1,3,2)
    return ((i,j,k) for i=1:length(x.values[1]),
                        k=1:length(x.values[3]),
                        j=1:length(x.values[2]))
  elseif order==(2,1,3)
    return ((i,j,k) for j=1:length(x.values[2]),
                        i=1:length(x.values[1]),
                        k=1:length(x.values[3]))
  elseif order==(2,3,1)
    return ((i,j,k) for j=1:length(x.values[2]),
                        k=1:length(x.values[3]),
                        i=1:length(x.values[1]))
  elseif order==(3,1,2)
    return ((i,j,k) for k=1:length(x.values[3]),
                        i=1:length(x.values[1]),
                        j=1:length(x.values[2]))
  elseif order==(3,2,1)
    return ((i,j,k) for k=1:length(x.values[3]),
                        j=1:length(x.values[2]),
                        i=1:length(x.values[1]))
  elseif order==(1,2)
    return ((i,j) for   i=1:length(x.values[1]),
                        j=1:length(x.values[2]))
  elseif order==(2,1)
    return ((i,j) for   j=1:length(x.values[2]),
                        i=1:length(x.values[1]))
  elseif order==(1,)
    return ((i) for     i=1:length(x.values[1]))
  else
    throw(ArgumentError())
  end
end

function eachstep(x::LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                         order=ntuple(i->i,Nstep)) where {Nparam,Nmeas,Nmdim,Nstep}
  eachstep(x.stepvalues,order)
end
function eachstep(x::LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                         order::NTuple{Nstep,AbstractString}) where {Nparam,Nmeas,Nmdim,Nstep}
  eachstep(x,ntuple(i->findfirst(x.stepnames,order[i]),Nstep))
end

struct ResultNamesIndices
  isparameter :: Array{Bool,1}
  parametervalue :: Array{Float64,1}
  measurementindex :: Array{Int,1}
end
ResultNamesIndices(n::Int) = ResultNamesIndices(Array{Bool}(undef, n),Array{Float64}(undef, n),Array{Int}(undef, n))

"""
```julia
perlineiterator(simulation, <keyword arguments>)
```
Retruns iterator that flattens multidimensional data.

**Keyword Arguments**

- `steporder`     -- specify order of steps
- `resultnames`   -- specify parameters and measurements for output
- `header`        -- `true` to make first line header

The step order defaults to the order the step values appear in the circuit file.
Step order can be specified by passing an array of step names.  By default
there is one column for each step, measurement, and parameter.  The desired
measurements and parameters can be set by passing an array of names to
resultnames.

```julia
# write CSV with headers
using DelimitedFiles
writedlm("test.csv",perlineiterator(circuit2,header=true),',')
```
"""
function perlineiterator(x :: LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep};
                         steporder = ntuple(i->i,Nstep),
                         resultnames = (x.parameternames...,
                                        x.measurementnames...,),
                         header::Bool = false) where {Nparam,Nmeas,Nmdim,Nstep}
  _perlineiterator(x,steporder,resultnames,header)
end
perlineiterator(::NonSteppedSimulation;
                steporder=nothing,
                resultnames=nothing,
                header=nothing) = ()
function _perlineiterator(x :: LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                         steporder,
                         resultnames,
                         header::Bool) where {Nparam,Nmeas,Nmdim,Nstep}
  if header
    return Iterators.flatten(([headerline(x,steporder,resultnames)],
          _perlineiterator(x,steporder,resultnames)))
  else
    return _perlineiterator(x,steporder,resultnames)
  end
end
function _perlineiterator(x :: LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                         steporder,
                         resultnames) where {Nparam,Nmeas,Nmdim,Nstep}
  length(steporder)!=Nstep && throw(ArgumentError("must include all steps"))
  _perlineiterator(x,
                  ntuple(i->something(findfirst(isequal(steporder[i]), x.stepnames), 0),Nstep),
                  resultnames)
end

function _perlineiterator(x :: LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                         steporder::NTuple{Nstep,Int},
                         resultnames) where {Nparam,Nmeas,Nmdim,Nstep}
  allunique(steporder) || throw(ArgumentError("steps must be unique"))
  run!(x)
  resultnameslength = length(resultnames)
  rni = ResultNamesIndices(resultnameslength)
  for i in eachindex(resultnames)
    if haskey(x.parameterdict,resultnames[i])
      rni.isparameter[i] = true
      rni.parametervalue[i] = x[resultnames[i]]
    elseif haskey(x.measurementdict,resultnames[i])
      rni.isparameter[i] = false
      rni.measurementindex[i] = x.measurementdict[resultnames[i]]
    else
      throw(ArgumentError("result name \"$(resultnames[i])\" not found"))
    end
  end
  n = Nstep + resultnameslength
  return (ntuple(j->pliresultvalue(x,rni,i,j,steporder),n) for i in eachstep(x,steporder))
end
function pliresultvalue(x::LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
                        rni,i,j,steporder) where {Nparam,Nmeas,Nmdim,Nstep}
  if j<=Nstep
    return x.stepvalues.values[steporder[j]][i[steporder[j]]]
  elseif rni.isparameter[j-Nstep]
    return rni.parametervalue[j-Nstep]
  else
    return x.measurementvalues[i...,rni.measurementindex[j-Nstep]]
  end
end

function headerline(x,steporder::NTuple{N,Int},resultnames) where N
  (ntuple(i->x.stepnames[steporder[i]],N)...,resultnames...)
end
function headerline(x,steporder,resultnames)
  (steporder...,resultnames...)
end


text:
@testset "empty matrix" begin
    A = InteractionMatrix(Int[], Int[], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (0, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Int}(Int[], Int[], +)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (0, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Float64}(Int[], Int[], /)
    @test A isa AbstractArray{Float64, 2}
    @test size(A) == (0, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]
end

@testset "no rowelems" begin
    A = InteractionMatrix(Int[], Int[42, 40, 47], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (0, 3)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Int}(Int[], Int[42, 40, 47], +)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (0, 3)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Float64}(Int[], Int[42, 40, 47], /)
    @test A isa AbstractArray{Float64, 2}
    @test size(A) == (0, 3)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]
end

@testset "no colelems" begin
    A = InteractionMatrix(Int[42, 40, 47], Int[], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (3, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Int}(Int[42, 40, 47], Int[], +)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (3, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]

    A = InteractionMatrix{Float64}(Int[42, 40, 47], Int[], /)
    @test A isa AbstractArray{Float64, 2}
    @test size(A) == (3, 0)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[1]
end

@testset "single value" begin
    A = InteractionMatrix(Int[2], Int[42], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (1, 1)
    @test A[] isa Int
    @test A[] == 11
    @test A[1] isa Int
    @test A[1] == 11
    @test A[1, 1] isa Int
    @test A[1, 1] == 11
    @test_throws BoundsError A[2]
    @test_throws BoundsError A[1, 2]
    @test_throws BoundsError A[2, 1]
    @test_throws BoundsError A[2, 2]
    @test_throws CanonicalIndexError A[] = 13
    @test_throws CanonicalIndexError A[1] = 13
    @test_throws CanonicalIndexError A[1, 1] = 13

    A = InteractionMatrix{Int}(Int[2], Int[42], -)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (1, 1)
    @test A[] isa Int
    @test A[] == -40
    @test A[1] isa Int
    @test A[1] == -40
    @test A[1, 1] isa Int
    @test A[1, 1] == -40
    @test_throws BoundsError A[2]
    @test_throws BoundsError A[1, 2]
    @test_throws BoundsError A[2, 1]
    @test_throws BoundsError A[2, 2]
    @test_throws CanonicalIndexError A[] = 13
    @test_throws CanonicalIndexError A[1] = 13
    @test_throws CanonicalIndexError A[1, 1] = 13

    A = InteractionMatrix{Int}(Int[42], Int[2], -)
    @test A[1] == 40

    A = InteractionMatrix{Int}(Int[42], Int[2], /)
    @test_throws TypeError A[1]

    A = InteractionMatrix{Float64}(Int[42], Int[2], /)
    @test A isa AbstractArray{Float64, 2}
    @test A[1] isa Float64
    @test A[1] == 21.0
end

@testset "single rowelem" begin
    A = InteractionMatrix(Int[2], Int[42, 40, 47], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (1, 3)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == 11
    @test A[2] isa Int
    @test A[2] == 11
    @test A[3] isa Int
    @test A[3] == 11
    @test A[1, 1] isa Int
    @test A[1, 1] == 11
    @test A[1, 2] isa Int
    @test A[1, 2] == 11
    @test A[1, 3] isa Int
    @test A[1, 3] == 11
    @test_throws BoundsError A[4]
    @test_throws BoundsError A[1, 4]
    @test_throws BoundsError A[2, 1]

    A = InteractionMatrix{Int}(Int[2], Int[42, 40, 47], -)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (1, 3)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == -40
    @test A[2] isa Int
    @test A[2] == -38
    @test A[3] isa Int
    @test A[3] == -45
    @test A[1, 1] isa Int
    @test A[1, 1] == -40
    @test A[1, 2] isa Int
    @test A[1, 2] == -38
    @test A[1, 3] isa Int
    @test A[1, 3] == -45
    @test_throws BoundsError A[4]
    @test_throws BoundsError A[1, 4]
    @test_throws BoundsError A[2, 1]

    A = InteractionMatrix{Int}(Int[2], Int[42, 40, 47], \)
    @test_throws TypeError A[1]
    @test_throws TypeError A[2]
    @test_throws TypeError A[3]
    @test_throws TypeError A[1, 1]
    @test_throws TypeError A[1, 2]
    @test_throws TypeError A[1, 3]

    A = InteractionMatrix{Float64}(Int[2], Int[42, 40,47], \)
    @test A[1] isa Float64
    @test A[1] == 21.0
    @test A[2] isa Float64
    @test A[2] == 20.0
    @test A[3] isa Float64
    @test A[3] == 23.5
    @test A[1, 1] isa Float64
    @test A[1, 1] == 21.0
    @test A[1, 2] isa Float64
    @test A[1, 2] == 20.0
    @test A[1, 3] isa Float64
    @test A[1, 3] == 23.5
end

@testset "single colelem" begin
    A = InteractionMatrix(Int[42, 40, 47], Int[2], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (3, 1)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == 11
    @test A[2] isa Int
    @test A[2] == 11
    @test A[3] isa Int
    @test A[3] == 11
    @test A[1, 1] isa Int
    @test A[1, 1] == 11
    @test A[2, 1] isa Int
    @test A[2, 1] == 11
    @test A[3, 1] isa Int
    @test A[3, 1] == 11
    @test_throws BoundsError A[4]
    @test_throws BoundsError A[4, 1]
    @test_throws BoundsError A[1, 2]

    A = InteractionMatrix{Int}(Int[42, 40, 47], Int[2], -)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (3, 1)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == 40
    @test A[2] isa Int
    @test A[2] == 38
    @test A[3] isa Int
    @test A[3] == 45
    @test A[1, 1] isa Int
    @test A[1, 1] == 40
    @test A[2, 1] isa Int
    @test A[2, 1] == 38
    @test A[3, 1] isa Int
    @test A[3, 1] == 45
    @test_throws BoundsError A[4]
    @test_throws BoundsError A[4, 1]
    @test_throws BoundsError A[1, 2]

    A = InteractionMatrix{Int}(Int[42, 40, 47], Int[2], /)
    @test_throws TypeError A[1]
    @test_throws TypeError A[2]
    @test_throws TypeError A[3]
    @test_throws TypeError A[1, 1]
    @test_throws TypeError A[2, 1]
    @test_throws TypeError A[3, 1]

    A = InteractionMatrix{Float64}(Int[42, 40,47], Int[2], /)
    @test A[1] isa Float64
    @test A[1] == 21.0
    @test A[2] isa Float64
    @test A[2] == 20.0
    @test A[3] isa Float64
    @test A[3] == 23.5
    @test A[1, 1] isa Float64
    @test A[1, 1] == 21.0
    @test A[2, 1] isa Float64
    @test A[2, 1] == 20.0
    @test A[3, 1] isa Float64
    @test A[3, 1] == 23.5
end

@testset "ordinary" begin
    A = InteractionMatrix(Int[2, 5], Int[42, 40], 11)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (2, 2)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == 11
    @test A[2] isa Int
    @test A[2] == 11
    @test A[3] isa Int
    @test A[3] == 11
    @test A[4] isa Int
    @test A[4] == 11
    @test A[1, 1] isa Int
    @test A[1, 1] == 11
    @test A[2, 1] isa Int
    @test A[2, 1] == 11
    @test A[1, 2] isa Int
    @test A[1, 2] == 11
    @test A[2, 2] isa Int
    @test A[2, 2] == 11

    A = InteractionMatrix{Int}(Int[2, 5], Int[42, 40], -)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (2, 2)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == -40
    @test A[2] isa Int
    @test A[2] == -37
    @test A[3] isa Int
    @test A[3] == -38
    @test A[4] isa Int
    @test A[4] == -35
    @test A[1, 1] isa Int
    @test A[1, 1] == -40
    @test A[2, 1] isa Int
    @test A[2, 1] == -37
    @test A[1, 2] isa Int
    @test A[1, 2] == -38
    @test A[2, 2] isa Int
    @test A[2, 2] == -35

    A = InteractionMatrix{Int}(Int[42, 40], Int[2, 5], -)
    @test A isa AbstractArray{Int, 2}
    @test size(A) == (2, 2)
    @test_throws BoundsError A[]
    @test A[1] isa Int
    @test A[1] == 40
    @test A[2] isa Int
    @test A[2] == 38
    @test A[3] isa Int
    @test A[3] == 37
    @test A[4] isa Int
    @test A[4] == 35
    @test A[1, 1] isa Int
    @test A[1, 1] == 40
    @test A[2, 1] isa Int
    @test A[2, 1] == 38
    @test A[1, 2] isa Int
    @test A[1, 2] == 37
    @test A[2, 2] isa Int
    @test A[2, 2] == 35

    A = InteractionMatrix{Int}(Int[42, 40], Int[2, 5], /)
    @test_throws TypeError A[1]
    @test_throws TypeError A[2]
    @test_throws TypeError A[3]
    @test_throws TypeError A[4]
    @test_throws TypeError A[1, 1]
    @test_throws TypeError A[2, 1]
    @test_throws TypeError A[1, 2]
    @test_throws TypeError A[2, 2]

    A = InteractionMatrix{Float64}(Int[42, 40], Int[2, 5], /)
    @test A[1] isa Float64
    @test A[1] == 21.0
    @test A[2] isa Float64
    @test A[2] == 20.0
    @test A[3] isa Float64
    @test A[3] == 8.4
    @test A[4] isa Float64
    @test A[4] == 8.0
    @test A[1, 1] isa Float64
    @test A[1, 1] == 21.0
    @test A[2, 1] isa Float64
    @test A[2, 1] == 20.0
    @test A[1, 2] isa Float64
    @test A[1, 2] == 8.4
    @test A[2, 2] isa Float64
    @test A[2, 2] == 8.0
end


text:
include("../btg.jl")
using LinearAlgebra
"""
Precompute some quantities
diffs - array of differences between x_i and x_j, used to compute partial_Σθ_theta for various theta
"""
function globalvars(btg::btg)
    n = getNumPts(btg.trainingData)
    d = getDimension(btg.trainingData)
    x = getPosition(btg.trainingData)
    diffs = Array{Array{T, 2} where T, 2}(undef, n, n)
    for i = 1:n
        for j = 1:n
            diffs[i, j] = x[i:i, :] - x[j:j, :]
        end
    end
    return diffs
end

"""
Takes dot product of vectir a with each entry in matrix of vectors b and stores
result in corresponding entry in b
"""
function entrywise_dot(a, b)
    #@assert size(a, 2) == size(b[1, 1], 2)
    #@assert size(a, 1) == size(b[1, 1], 1)
    a = reshape(a, size(b[1, 1], 1), size(b[1, 1], 2))
    c = copy(b)
    for i = 1:size(c, 1)
        for j = 1:size(c, 2)
            c[i,j] = a .* c[i, j]
        end
    end
    return c
end

function hadamard(a, b)
    @assert size(a, 1) == size(b, 1)
    @assert size(a, 2) == size(b, 2)
    m = size(a, 1)
    n = size(a, 2)
    c = similar(a)
    for i =1:m
        for j = 1:n
            c[i, j] = a[i, j] .* b[i, j]
        end
    end
    return
end

"""
Squares each entry of vector v
"""
function vector_square(v)
    w = copy(v)
    for i = 1:length(v)
        w[i] = v[i]^2
    end
    return w
end

"""
log(p(theta, lambda|z))
"""
function tdist_mle(btg::btg, theta, lambda)
    (x, Fx, y, _, n, p) = unpack(btg.trainingData) #unpack training data
    function Σθ(theta)
        return correlation(btg.k, theta, getPosition(btg.trainingData); jitter = 1e-6)
    end
    g = btg.g
    dg = (y, λ) -> partialx(btg.g, y, λ)
    try
        global choleskyΣθ = cholesky(Hermitian(Σθ(theta)))
    catch e
        println("problematic kernel matrix")
        display(minimum(eigvals(Hermitian(Σθ(theta)))))
        println("theta")
        display(theta)
    end
    gλz = g(y, lambda)
    Σθ_inv_y = choleskyΣθ\gλz
    Σθ_inv_X = choleskyΣθ\Fx
    choleskyXΣX = cholesky(Hermitian(Fx'*Σθ_inv_X))
    pθ = btg.priorθ
    pλ = btg.priorλ
    function βhat()
        return choleskyXΣX\(Fx'*Σθ_inv_y)
    end
    βhat = βhat()
    function qtilde()
        qtilde =  gλz'*Σθ_inv_y  - 2*gλz'*Σθ_inv_X*βhat + βhat'*Fx'*Σθ_inv_X*βhat
    end
    qtilde = qtilde()
    logdetΣθ = logdet(choleskyΣθ)
    logdetXΣX = logdet(choleskyXΣX)

    logJ = z -> reduce(+, map(x -> log(dg(x, lambda)), z))
    logJ_ret = logJ(y)

    logpθ = logProb(pθ, theta)
    logpλ = logProb(pλ, lambda)

    logprob = logpθ + logpλ + (1-p/n) * logJ_ret - (n-p)/2 * log(qtilde) - 0.5 * logdetXΣX -0.5 * logdetΣθ
    @info "logpθ", logpθ
    @info "logpλ", logpλ
    @info "(1-p/n) * logJ_ret", (1-p/n) * logJ_ret
    @info "- (n-p)/2 * log(qtilde)", - (n-p)/2 * log(qtilde)
    @info "- 0.5 * logdetXΣX", - 0.5 * logdetXΣX
    @info "-0.5 * logdetΣθ", -0.5 * logdetΣθ
    @info "logprob", logprob
    return -logprob
end

"""
Maximizes p(theta, lambda|z)
"""
function partial_tdist_mle(btg::btg, theta, lambda)
    ##############################
    ######### derivatives #########
    ##############################
    function partial_Σθ_theta(theta)
        diffs = globalvars(btg);
        square = vector_square.(diffs);
        return -0.5 * map(x ->exp.(x),- 0.5 * sum.(entrywise_dot(theta, square))) .* square;
    end

    function partial_detΣθ_theta()

    end

    function partial_det_XΣX_theta()

    end

    function partial_q_theta()
        g = btg.g
        dg = (y, λ) -> partialx(btg.g, y, λ)
    end

    function partial_q_lambda()
    end

    function partial_J_lambda()
    end

    function partial_p_theta()
        return partialx(btg.priorθ, theta)
    end

    function partial_p_lambda()
        return partialx(btg.priorλ, lambda)
    end

    function compute_betahat_prime_theta(choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}, expr_mid, Σθ_prime, X::Array{Float64,2}, gλz, Σθ_inv_X)
        AA = choleskyXΣX\(expr_mid)*(choleskyXΣX\(X'*(choleskyΣθ\gλz)))
        BB = - (choleskyXΣX\(X'*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\gλz)))))
        βhat_prime_theta = AA + BB
        βhat_prime_theta = reshape(βhat_prime_theta, size(βhat_prime_theta, 1), size(βhat_prime_theta, 2)) #turn 1D array into 2D array
    end

    function compute_qtilde_prime_theta(gλz, X::Array{Float64,2}, βhat, βhat_prime_theta, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Σθ_prime::Array{Float64,2})
        meanvv = gλz - X*βhat
        rr = X*βhat_prime_theta
        AA = (-rr)' * (choleskyΣθ \ meanvv)
        BB = - meanvv' * (choleskyΣθ \ (Σθ_prime * (choleskyΣθ \ meanvv)))
        CC =  meanvv' * (choleskyΣθ \ (-rr))
        qtilde_prime_theta = AA .+ BB .+ CC
    end

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
#
# General utils
#

#
# Utils for Akl Toussaint
#
"""
Creates the "T inverse" matrix described by wikipedia and
above that will right multiply points to get the Barycentric
coordinates
"""
function create_t_inv(a::Point{2}, b::Point{2}, c::Point{2})
    # Pull out points
    a1, a2 = a
    b1, b2 = b
    c1, c2 = c

    # Build t inverse matrix
    temp = (a1*b2 - a1*c2 - c1*b2) - (b1*a2 - b1*c2 - c1*a2)
    return [(b2-c2)/temp (c1-b1)/temp; (c2-a2)/temp (a1-c1)/temp]
end

# One liner to check if λ values are in (0, 1) and sum to less than 1
@inline check_λ(λ1, λ2) = (0 < λ1 < 1) && (0 < λ2 < 1) && (λ1 + λ2 < 1) ? false : true

"""
Determines which points in a vector of points lie inside or outside the
triangle determined by the triangle abc.
"""
function pt_out_tri{T}(x::Vector{Point{2, T}}, a::Point{2}, b::Point{2}, c::Point{2})
    # Get number of points in x
    npts = length(x)

    # Allocate space
    out = Array(Bool, npts)

    # Create our "T inverse" matrixs
    tinv = create_t_inv(a, b, c)

    @inbounds for n=1:npts
        # Get the vector we need to multiply by T inverse to get
        # the Barycentric coordinates
        temp = x[n] - c

        # Get Barycentric coordinates
        λ1, λ2 = tinv*temp
        out[n] = check_λ(λ1, λ2)
    end

    return out
end

#
# Wrap around set of points
#
"""
Takes a vector of points (which must be sorted either by
relative angle or x-value) and wraps around them by assuring
that they only make counter-clockwise turns.
"""
function wrappoints{T<:Real}(p::Vector{Point{2, T}})

    # Get number of points
    npts = length(p)

    # Put first two points inside the convex hull
    ch = [p[1], p[2]]
    sizehint!(ch, npts)

    # Starting value for number of extreme points
    m = 2

    # Will check each point
    @inbounds for i=3:npts
        # If it makes a clockwise turn then the
        # last element of the convex hull doesn't
        # belong, so pop it out. Make sure that at
        # least two elements are always in it
        while m>1 && cw(ch[m-1], ch[m], p[i])
            pop!(ch)
            m = m-1
        end
        push!(ch, p[i])
        m = m+1
    end

    return ch
end

#
# Utils for Monotone Chain
#
function split_xL_xU{T<:Real}(psort::Vector{Point{2, T}})
    # Allocate space to identify where points belong
    # We make endpoints true because we want them in both sets
    npts = length(psort)
    xL_bool = Array(Bool, npts)
    xL_bool[1] = true; xL_bool[end] = true
    xU_bool = Array(Bool, npts)
    xU_bool[1] = true; xU_bool[end] = true

    # Get the endpoints and slope
    pmin = psort[1]; pmax = psort[end]
    ls = LineSegment(pmin, pmax)

    # We check whether it is above or below the line between
    # the points, but could also check whether it was clockwise
    # or counter-clockwise (above/below seems like less operations)
    # Iterate through all points except endpoints
    @inbounds for i=2:npts-1
        # Pull out x and y values for ith point
        p_i = psort[i]
        xi, yi = p_i[1], p_i[2]

        # Compute where y on line would be
        yi_line = evaluatey(ls, xi)

        if yi_line <= yi
            xL_bool[i] = false
            xU_bool[i] = true
        else
            xL_bool[i] = true
            xU_bool[i] = false
        end
    end
    xL = psort[xL_bool]
    xU = psort[xU_bool]

    return xL, xU
end

#
# Additional functions for Point{2, T} that will be useful
#

# Matrix times point
function _unsafe_mult_2by2_point{T}(A::Matrix{T}, x::Point{2, T})
    x1 = x[1]
    x2 = x[2]
    out = (A[1, 1]*x1 + A[1, 2]*x2, A[2, 1]*x1 + A[2, 2]*x2)
    return out
end

function _unsafe_mult_Nby2_point{T}(A::Matrix{T}, x::Point{2, T})
    nr, nc = size(A)
    x1 = x[1]
    x2 = x[2]

    out = Vector(T, nr)
    for i=1:nr
        out[i] = A[i, 1]*x1+A[i, 2]*x2
    end
    return out
end

function *(A::Matrix, x::Point{2})
    nr, nc = size(A)

    # Depending on size call the right function
    if (nc==nr==2)
        out = _unsafe_mult_2by2_point(A, x)
    elseif (nc==2)
        out = _unsafe_mult_nby2(A, x)
    else
        error("Size is wrong. Number of cols should be 2 (currently is $nc)")
    end

    return out
end

function extremax{T}(p::Vector{Point{2, T}})
    # Get number of points
    npts = length(p)

    # Set initial values
    minval, minind = Inf, 0
    maxval, maxind = -Inf, 0

    for i=1:npts
        cx, cy = p[i]
        minval, minind = cx < minval ? (cx, i) : (minval, minind)
        maxval, maxind = cx > maxval ? (cx, i) : (maxval, maxind)
    end

    return minval, minind, maxval, maxind
end

function extremay{T}(p::Vector{Point{2, T}})
    # Get number of points
    npts = length(p)

    # Set initial values
    minval, minind = Inf, 0
    maxval, maxind = -Inf, 0

    for i=1:npts
        cx, cy = p[i]
        minval, minind = cy < minval ? (cy, i) : (minval, minind)
        maxval, maxind = cy > maxval ? (cy, i) : (maxval, maxind)
    end

    return minval, minind, maxval, maxind
end

@inline distance{T}(a::Point{2, T}, b::Point{2, T}) = sqrt((b[1] - a[1])^2 - (b[2] - a[2])^2)

# Comparisons
isless{T}(x::Point{2, T}, y::Point{2, T}) = x[1] < y[1] ? true : (x[1]==y[1] && x[2]<y[2]) ? true : false
isequal{T}(x::Point{2, T}, y::Point{2, T}) = (x[1]==y[1]) && (x[2]==y[2]) ? true : false

# Orientation
@inline orientation(a::Point{2}, b::Point{2}, c::Point{2}) =
    (b[1]-a[1])*(c[2]-a[2]) - (b[2]-a[2])*(c[1]-a[1])
cw(a::Point{2}, b::Point{2}, c::Point{2}) = orientation(a, b, c) < 0.0 ? true : false
ccw(a::Point{2}, b::Point{2}, c::Point{2}) = orientation(a, b, c) >= 0.0 ? true : false


text:
# Nocedal (1980) based on
# Specialized S and Y matrices are constructed and updated
struct TwoLoop end
struct CompactLimited end
struct LBFGS{TA, F, T, TP}
  approx::TA
	type::F
  memory::T
  P::TP
end
summary(lbfgs::LBFGS{Inverse}) = "Inverse LBFGS"
summary(lbfgs::LBFGS{Direct}) = "Direct LBFGS"

hasprecon(::LBFGS{<:Inverse, <:Any, <:Any, <:Nothing}) = NoPrecon()
hasprecon(::LBFGS{<:Inverse, <:Any, <:Any, <:Any}) = HasPrecon()

LBFGS(m::Int=5) = LBFGS(Inverse(), TwoLoop(), m, nothing)
LBFGS(approx, m=5) = LBFGS(approx, TwoLoop(), m, nothing)
"""
	q holds gradient at current state
	history is a named tuple of S and Y
	memory is the number of past (s,y) pairs we have
"""
function find_direction!(scheme::LBFGS{<:Inverse, <:TwoLoop}, q,
                  qnvars,
                  memory,
                  scaling,
                  precon)

    S, Y = qnvars.S, qnvars.Y
	  α, ρ = qnvars.α, qnvars.ρ
    d    = qnvars.d

    # Backward pass
    @inbounds for i in memory:-1:1
        α[i] = ρ[i] * real(dot(S[i], q))
        q .-= α[i] .* Y[i]
    end
    # Copy scaled or preconditioned q into s for forward pass
    if memory > 0
        if scaling isa InitialScaling{<:ShannoPhua} && precon isa Nothing # we need a pair to scale
            k = scaling(S[memory], Y[memory])
            @. d = k*q
        elseif !(precon isa Nothing)
            mul!(d, precon, q)
        end
    else
        d .= q
    end
    # Forward pass
    @inbounds for i in 1:memory
        β = ρ[i] * real(dot(Y[i], d))
        d .+= S[i] .* (α[i] - β)
    end
    # Negate search direction
    negate(InPlace(), d)
end
function update_obj!(problem, qnvars, α, x, ∇fx, z, ∇fz, current_memory, scheme::LBFGS{<:Inverse, <:TwoLoop}, scale=nothing)
    # Calculate final step vector and update the state
    fz, ∇fz = upto_gradient(problem, ∇fz, z)
    # add Project gradient

    # Quasi-Newton update
    qnvars = update!(scheme, qnvars, ∇fx, ∇fz, current_memory)

    return fz, ∇fz, qnvars
end
@inbounds function update!(scheme::LBFGS{<:Inverse, <:TwoLoop}, qnvars, ∇fx, ∇fz, current_memory)
  S, Y, ρ = qnvars.S, qnvars.Y, qnvars.ρ
  n = length(S)
  m = min(n, 1+current_memory)
  if current_memory == n
    s1, y1 = S[1], Y[1] # hoisting, no allocation
    for i = 2:n
      S[i-1] = S[i]
      Y[i-1] = Y[i]
      ρ[i-1] = ρ[i]
    end
    S[end] = s1
    Y[end] = y1
  end
  @. Y[m] = ∇fz - ∇fx
  ρ[m] = 1/dot(S[m], Y[m])
  TwoLoopVars(qnvars.d, S, Y, qnvars.α, ρ)
end


text:
# Model/chain/residue and model/fragment hierarchies

module Hierarchies

using ..Dorothy
using ..Dorothy.Graphs
using ..Dorothy.Multicollections
using ..Dorothy.Utils

export
		mcrptree, mcrp, chains, residues, mfptree, mfp, fragments,
		chainat, residueat, fragmentat

function mcrptree(model::ParticleCollection)
	n = length(model)
	chainids = get(model, :chainids, Repeated("", n))
	resids = get(model, :resids, Repeated(0, n))
	mcrptree(chainids, resids)
end

function mcrptree(chainids::AbstractVector{<:AbstractString},
		resids::AbstractVector{<:Integer})
	n = length(chainids)
	@boundscheck length(resids) == n ||
			error("size mismatch between chainid and resid arrays")
	model = Vector{Vector{Int}}[]
	if n > 0
		residue = [1]
		chain = [residue]
		push!(model, chain)
		lastchainid = chainids[1]
		lastresid = resids[1]
		for i = 2:n
			thischainid = chainids[i]
			thisresid = resids[i]
			if thischainid != lastchainid
				residue = [i]
				chain = [residue]
				push!(model, chain)
			elseif thisresid != lastresid
				residue = [i]
				push!(chain, residue)
			else
				push!(residue, i)
			end
			lastchainid = thischainid
			lastresid = thisresid
		end
	end
	model
end

mcrp(model::ParticleCollection) = mcrp(mcrptree(model), length(model))

mcrp(tree::Vector{Vector{Vector{Int}}}, n::Integer) = H3Hierarchy(tree, n)

chains(model::ParticleCollection) = chains(model, mcrp(model))

chains(model::ParticleCollection, mcrp::H3Hierarchy) =
		H3IteratorH2(model, mcrp)

residues(model::ParticleCollection) = residues(model, mcrp(model))

residues(model::ParticleCollection, mcrp::H3Hierarchy) =
		H2Iterator(model, mcrp.flath2)

residues(model::ParticleCollection, chainindex::Integer) =
		residues(model, chainindex, mcrp(model))

residues(model::ParticleCollection, ichain::Integer, mcrp::H3Hierarchy) =
		H3IteratorH1(model, ichain, mcrp)

mfptree(model::ParticleCollection) =
		mfptree(get(model, :topology, Graph(length(model))))

function mfptree(G::AbstractGraph)
	n = length(G)
	model = Vector{Int}[]
	if n > 0
		reachable = trues(n)
		i = 1
		while true
			fragment = sort!(connected!(Int[], G, i, reachable))
			push!(model, fragment)
			i = findnext(reachable, i+1)
			if i == nothing
				break
			end
		end
	end
	model
end

mfp(model::ParticleCollection) = mfp(mfptree(model), length(model))

mfp(tree::Vector{Vector{Int}}, n::Integer) = H2Hierarchy(tree, n)

fragments(model::ParticleCollection) = fragments(model, mfp(model))

fragments(model::ParticleCollection, mfp::H2Hierarchy) = H2Iterator(model, mfp)

function chainat(model::ParticleCollection, i::Integer)
	@boundscheck checkbounds(model, i)
	chainids = get(model, :chainids, Repeated("", length(model)))
	view(model, chainat(chainids, i))
end

function chainat(chainids::AbstractArray{<:AbstractString}, i::Integer)
	@boundscheck checkbounds(chainids, i)
	n = length(chainids)
	thischainid = chainids[i]
	firsti = i - 1
	while firsti > 0 && chainids[firsti] == thischainid
		firsti -= 1
	end
	firsti += 1
	lasti = i + 1
	while lasti <= n && chainids[lasti] == thischainid
		lasti += 1
	end
	lasti -= 1
	firsti:lasti
end

function residueat(model::ParticleCollection, i::Integer)
	@boundscheck checkbounds(model, i)
	n = length(model)
	chainids = get(model, :chainids, Repeated("", n))
	resids = get(model, :resids, Repeated(0, n))
	view(model, residueat(chainids, resids, i))
end

function residueat(chainids::AbstractArray{<:AbstractString},
		resids::AbstractArray{<:Integer}, i::Integer)
	n = length(chainids)
	@boundscheck begin
		length(resids) == n ||
				error("size mismatch between chainid and resid arrays")
		checkbounds(chainids, i)
	end
	thischainid = chainids[i]
	thisresid = resids[i]
	firsti = i - 1
	while firsti > 0 && chainids[firsti] == thischainid &&
			resids[firsti] == thisresid
		firsti -= 1
	end
	firsti += 1
	lasti = i + 1
	while lasti <= n && chainids[lasti] == thischainid &&
			resids[lasti] == thisresid
		lasti += 1
	end
	lasti -= 1
	firsti:lasti
end

function fragmentat(model::ParticleCollection, i::Integer)
	@boundscheck checkbounds(model, i)
	topology = get(model, :topology, Graph(length(model)))
	view(model, fragmentat(topology, i))
end

fragmentat(G::AbstractGraph, i::Integer) = sort!(connected(G, i))

end # module


text:
using ChaosTools
using Test, StaticArrays

println("\nTesting orbit diagrams...")

@testset "Discrete Orbit Diagram" begin
    @testset "1D" begin
        ds = Systems.logistic(0.5112661)
        i = 1
        parameter = 1
        pvalues = 2:0.01:4
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test output[1][1] ≈ 0.5
        @test output[2][1] ≈ 0.502487562
        @test output[end][1] != output[end][2] != output[end][3]

        @testset "specific range" begin
        ulims = (0.0, 0.5)
        ics = zeros(length(pvalues)); ics[end] = 0.5112661
        output = orbitdiagram(
            ds, i, parameter, pvalues;
            n = n, Ttr = Ttr, u0 = ics, ulims = ulims
        )
        @test all(iszero, output[1])
        @test !all(iszero, output[end])
        @test all(≤(0.5), output[end])
        end
    end

    @testset "2D" begin
        ds = Systems.standardmap()
        i = 2
        parameter = 1
        pvalues = 0:0.05:2
        ics = [0.001rand(2) for m in pvalues]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr, u0=ics)
        @test length(output[1]) == n
    end

    @testset "Vector values collection" begin
        ds = Systems.coupledstandardmaps(4)
        i = SVector(2, 3)
        parameter = 1
        pvalues = 0:0.05:2
        ics = [0.001rand(8) for m in pvalues]
        n = 50
        Ttr = 500
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output) == length(pvalues)
        @test length(output[1]) == n
        @test length(output[1][1]) == 2
    end

    @testset "IIP" begin
        ds = Systems.henon_iip()
        i = 2
        parameter = 1
        pvalues = [1.4, 1.5]
        ics = [rand(2), rand(2)]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
    end
end


text:
abstract type Kernel end

mutable struct Linear <: Kernel
    σ::Real
    function Linear(σ::Real)
        @assert σ > 0 "The output variance, σ, must be positive."
        new(σ)
    end
end

mutable struct RBF <: Kernel
    l::Real
    σ::Real
    function RBF(l::Real, σ::Real)
        @assert σ > 0 "The output variance, σ, must be positive."
        @assert l > 0 "The lengthscale, l, must be positive."
        new(l, σ)
    end
end

mutable struct RationalQuadratic <: Kernel
    a::Real
    l::Real
    σ::Real
    function RationalQuadratic(a::Real, l::Real, σ::Real)
        @assert a > 0 "The relative weighting of different lengthscales, a, must be positive."
        @assert σ > 0 "The output variance, σ, must be positive."
        @assert l > 0 "The lengthscale, l, must be positive."
        new(a, l, σ)
    end
end

mutable struct Periodic <: Kernel
    p::Real
    l::Real
    σ::Real
    function Periodic(p::Real, l::Real, σ::Real)
        @assert p > 0 "The period, p, must be positive."
        @assert σ > 0 "The output variance, σ, must be positive."
        @assert l > 0 "The lengthscale, l, must be positive."
        new(p, l, σ)
    end
end

mutable struct LocallyPeriodic <: Kernel
    p::Real
    l::Real
    σ::Real
    function LocallyPeriodic(p::Real, l::Real, σ::Real)
        @assert p > 0 "The period, p, must be positive."
        @assert σ > 0 "The output variance, σ, must be positive."
        @assert l > 0 "The lengthscale, l, must be positive."
        new(p, l, σ)
    end
end

abstract type Model end

mutable struct RLS{T <: AbstractFloat} <: Model
    y::Vector{T}
    w::Vector{T}
    λ::Real
    ŷ::Vector{T}
    σ::Vector{T}
end

mutable struct KRLS{T <: AbstractFloat} <: Model
    y::REALVEC
    x::REALMAT
    kernel::Kernel
    λ::Real
    budget::Integer
    mode::String
    μ::Vector{T}
    Σ::Matrix{T}
    Q::Matrix{T}
    ŷ::Vector{T}
    σ::Vector{T}
    #function KRLS(λ::Real, budget::Integer, forgetting::String)
    #    @assert λ <= 1 && λ > 0 "Forgetting factor, λ, must be ]0;1]."
    #    @assert budget > 0 "Budget size must be a positive integer."
    #    @assert (forgetting == "B2P" || forgetting == "UI") "Forgetting mode must be either B2P (back 2 prior) or UI (uncertainty injection)."
    #    new(λ, budget, forgetting)
    #end
end


text:
using Propagators: Cell, propagator, constraint, block,
  constant, quadratic, propagate!
using IntervalArithmetic

# h = gt²/2

function fallduration(t, h)
  block(t, h) do
    g = Cell()
    onehalf = Cell()
    t² = Cell()
    gt² = Cell()
    constant(9.789..9.832)(g)
    constant(1/2..1/2)(onehalf)
    quadratic(t, t²)
    constraint(*)(g, t², gt²)
    constraint(*)(onehalf, gt², h)
  end
end

# s_ba/h_ba == s/h

function similartriangles(s_ba, h_ba, s, h)
  block(s_ba, h_ba, s, h) do
    ratio = Cell()
    constraint(*)(s_ba, ratio, h_ba)
    constraint(*)(s, ratio, h)
  end
end

barometerheight = Cell()
barometershadow = Cell()
buildingheight = Cell()
buildingshadow = Cell()
similartriangles(barometershadow, barometerheight, buildingshadow, buildingheight)

buildingshadow[] = 54.9..55.1
barometerheight[] = 0.3..0.32
barometershadow[] = 0.36..0.37
propagate!()

# Initial guess based on barometer height/shadow
buildingheight[]

falltime = Cell()
fallduration(falltime, buildingheight)
propagate!()

# We have an idea of falltime even before measuring it
falltime[]
# Add our measurement
falltime[] = 2.9..3.1

propagate!()

# Building height with shadow + fall time
buildingheight[]

# We refined both of our measurements
falltime[]
barometerheight[]

# If we know the building height, it refines our measurements further
buildingheight[] = 45..45
propagate!()

falltime[]
barometerheight[]
barometershadow[]
buildingshadow[]


text:
# test to see that deleteHybridizationUpdate undoes all attributes
# prompted by Cecile finding cases when containRoot was not updated
# Claudia December 2015

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
    const DEBUG = true
end

if isdefined(:PhyloNetworks)
    PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")
else
    CHECKNET || error("need CHECKNET==true to test snaq in test_correctLik.jl")
end

tree = "(((((((1,2),3),4),5),(6,7)),(8,9)),10);"

#seed = 2738
seed = 56326

currT0 = readTopologyLevel1(tree);
## printEdges(currT0)
## printNodes(currT0)
## writeTopologyLevel1(currT0)
checkNet(currT0)
srand(seed);
besttree = deepcopy(currT0);

# ===== first hybridization ==========================
println("starting first hybridization")
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success || error("not able to place first hybridization")
#printEdges(besttree)
#printNodes(besttree)
println(writeTopologyLevel1(besttree))
net = deepcopy(besttree);
# test contain root
!net.edge[15].containRoot || error("edge 15 wrong contain Root")
!net.edge[19].containRoot || error("edge 19 wrong contain Root")
!net.edge[20].containRoot || error("edge 20 wrong contain Root")

# test inCycle
net.edge[12].inCycle == 11 || error("wrong incycle")
net.edge[13].inCycle == 11 || error("wrong incycle")
net.edge[16].inCycle == 11 || error("wrong incycle")
net.edge[18].inCycle == 11 || error("wrong incycle")
net.edge[19].inCycle == 11 || error("wrong incycle")
net.edge[20].inCycle == 11 || error("wrong incycle")

net.node[12].inCycle == 11 || error("wrong incycle")
net.node[13].inCycle == 11 || error("wrong incycle")
net.node[16].inCycle == 11 || error("wrong incycle")
net.node[17].inCycle == 11 || error("wrong incycle")
net.node[19].inCycle == 11 || error("wrong incycle")
net.node[20].inCycle == 11 || error("wrong incycle")

# test partition
length(net.partition) == 6 || error("wrong partition")
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [9,7,5,3,1,2,4,6,8] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [14] || error("wrong partition")

# ===== second hybridization ==========================
println("starting second hybridization")
success = false
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success || error("could not add second hybridization")
#printEdges(besttree)
#printNodes(besttree)
println(writeTopologyLevel1(besttree,true))
net = deepcopy(besttree);

# test contain root
!net.edge[15].containRoot || error("edge 15 wrong contain Root")
!net.edge[19].containRoot || error("edge 19 wrong contain Root")
!net.edge[20].containRoot || error("edge 20 wrong contain Root")
!net.edge[3].containRoot || error("edge 15 wrong contain Root")
!net.edge[1].containRoot || error("edge 19 wrong contain Root")
!net.edge[2].containRoot || error("edge 20 wrong contain Root")
!net.edge[23].containRoot || error("edge 15 wrong contain Root")
!net.edge[22].containRoot || error("edge 19 wrong contain Root")

# test inCycle
net.edge[12].inCycle == 11 || error("wrong incycle")
net.edge[13].inCycle == 11 || error("wrong incycle")
net.edge[16].inCycle == 11 || error("wrong incycle")
net.edge[18].inCycle == 11 || error("wrong incycle")
net.edge[19].inCycle == 11 || error("wrong incycle")
net.edge[20].inCycle == 11 || error("wrong incycle")

net.node[12].inCycle == 11 || error("wrong incycle")
net.node[13].inCycle == 11 || error("wrong incycle")
net.node[16].inCycle == 11 || error("wrong incycle")
net.node[17].inCycle == 11 || error("wrong incycle")
net.node[19].inCycle == 11 || error("wrong incycle")
net.node[20].inCycle == 11 || error("wrong incycle")

net.edge[9].inCycle == 13 || error("wrong incycle")
net.edge[7].inCycle == 13 || error("wrong incycle")
net.edge[5].inCycle == 13 || error("wrong incycle")
net.edge[22].inCycle == 13 || error("wrong incycle")
net.edge[23].inCycle == 13 || error("wrong incycle")

net.node[5].inCycle == 13 || error("wrong incycle")
net.node[7].inCycle == 13 || error("wrong incycle")
net.node[9].inCycle == 13 || error("wrong incycle")
net.node[21].inCycle == 13 || error("wrong incycle")
net.node[22].inCycle == 13 || error("wrong incycle")

# test partition
length(net.partition) == 10 || error("wrong partition")
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [14] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [3,1,2] || error("wrong partition")
[n.number for n in net.partition[7].edges] == [21] || error("wrong partition")
[n.number for n in net.partition[8].edges] == [8] || error("wrong partition")
[n.number for n in net.partition[9].edges] == [6] || error("wrong partition")
[n.number for n in net.partition[10].edges] == [4] || error("wrong partition")

## # ===== identify containRoot for net.node[21]
## net0=deepcopy(net);
## hybrid = net.node[21];
## hybrid.isBadDiamondI
## nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,hybrid);
## [n.number for n in edgesInCycle] == [23,9,7,5,22] || error("wrong identified cycle")
## [n.number for n in nodesInCycle] == [13,14,-4,-5,-6] || error("wrong identified cycle")
## edgesRoot = identifyContainRoot(net,hybrid);
## [n.number for n in edgesRoot] == [3,1,2] || error("wrong identified contain root")
## [e.containRoot for e in edgesRoot] == [false,false,false] || error("edges root should be false contain root")
## edges = hybridEdges(hybrid);
## [e.number for e in edges] == [22,23,3] || error("wrong identified edges")
## undoGammaz!(hybrid,net);
## othermaj = getOtherNode(edges[1],hybrid);
## othermaj.number == -6 || error("wrong othermaj")
## edgesmaj = hybridEdges(othermaj);
## [e.number for e in edgesmaj] == [22,5,4] || error("wrong identified edges")
## edgesmaj[3].containRoot || error("wrong edgesmaj[3] contain root")
## undoContainRoot!(edgesRoot);
## [e.containRoot for e in edgesRoot] == [true,true,true] || error("edges root should be false contain root")
## deleteHybrid!(hybrid,net,true, false)
## [e.containRoot for e in edgesRoot] == [true,true,true] || error("edges root should be false contain root")
## printEdges(net)

# ================= delete second hybridization =============================
println("starting deletion of second hybridization")
deleteHybridizationUpdate!(net,net.node[21], false,false);
length(net.partition) == 6 || error("wrong partition")
# 15,11,10,[9,7,5,3,1,2,4,6,8],17,14
[n.number for n in net.partition[1].edges] == [15] || error("wrong partition")
[n.number for n in net.partition[2].edges] == [11] || error("wrong partition")
[n.number for n in net.partition[3].edges] == [10] || error("wrong partition")
[n.number for n in net.partition[4].edges] == [17] || error("wrong partition")
[n.number for n in net.partition[5].edges] == [14] || error("wrong partition")
[n.number for n in net.partition[6].edges] == [3,1,2,8,6,4,9,7,5] || error("wrong partition")
#printNodes(net)
#printEdges(net)
#printPartitions(net)
println(writeTopologyLevel1(net))

# test contain root
!net.edge[15].containRoot || error("edge 15 wrong contain Root")
!net.edge[19].containRoot || error("edge 19 wrong contain Root")
!net.edge[20].containRoot || error("edge 20 wrong contain Root")
net.edge[1].containRoot || error("edge 1 wrong contain Root")
net.edge[2].containRoot || error("edge 2 wrong contain Root")
net.edge[3].containRoot || error("edge 3 wrong contain Root")

# test inCycle
net.edge[12].inCycle == 11 || error("wrong incycle")
net.edge[13].inCycle == 11 || error("wrong incycle")
net.edge[16].inCycle == 11 || error("wrong incycle")
net.edge[18].inCycle == 11 || error("wrong incycle")
net.edge[19].inCycle == 11 || error("wrong incycle")
net.edge[20].inCycle == 11 || error("wrong incycle")

net.node[12].inCycle == 11 || error("wrong incycle")
net.node[13].inCycle == 11 || error("wrong incycle")
net.node[16].inCycle == 11 || error("wrong incycle")
net.node[17].inCycle == 11 || error("wrong incycle")
net.node[19].inCycle == 11 || error("wrong incycle")
net.node[20].inCycle == 11 || error("wrong incycle")

# =============== delete first hybridization ===================
println("starting deletion of first hybridization")
deleteHybridizationUpdate!(net,net.node[19]);
checkNet(net)
length(net.partition) == 0 || error("wrong partition")
#printEdges(net)


text:
using Dates

@testset "Mot.jl" begin
    filename = "data\\test.mot"
    dc = DDR2import.Mot.read(filename)
    @test Dates.hour(dc["NAT1"].onlytimes[1].begintime) == 0
    @test Dates.second(dc["NAT1"].onlytimes[1].endtime) == 59
    #LAST
    @test dc["UKOV_FRA"].timetypes[1].type == "EX"
end


text:
using Test
using LocalProjections

using CSV
using CodecZlib: GzipDecompressorStream
using DataFrames
using LinearAlgebra: I
using LocalProjections: kron_fastl, kron_fastr, getscore, _geto,
    OLS, VarName, _makeYX, _firststage, _lp, _toint, _toname,
    Ridge, _basismatrix, _makeYSr, _makeP
using Tables: getcolumn

function exampledata(name::Union{Symbol,String})
    # Need to allow CSV v0.8 to work
    f = open(datafile(name)) |> GzipDecompressorStream |> read |> CSV.File
    return DataFrame(f, copycols=true)
end

const tests = [
    "utils",
    "lp",
    "slp",
    "vce",
    "irf"
]

printstyled("Running tests:\n", color=:blue, bold=true)

@time for test in tests
    include("$test.jl")
    println("\033[1m\033[32mPASSED\033[0m: $(test)")
end


text:
using MathOptInterface
const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

struct Solution
    ret_val::Int
    primal::Vector{Float64}
    dual_eq::Vector{Float64}
    dual_ineq::Vector{Float64}
    slack::Vector{Float64}
    objective_value::Float64
    dual_objective_value::Float64
    objective_constant::Float64
    solve_time::Float64
end
const OPTIMIZE_NOT_CALLED = -1
Solution() = Solution(OPTIMIZE_NOT_CALLED, Float64[], Float64[], Float64[],
                      Float64[], NaN, NaN, NaN, NaN)

# Used to build the data with allocate-load during `copy_to`.
# When `optimize!` is called, a the data is used to build `ECOSMatrix`
# and the `ModelData` struct is discarded
mutable struct ModelData
    m::Int # Number of rows/constraints
    n::Int # Number of cols/variables
    IA::Vector{Int} # List of conic rows
    JA::Vector{Int} # List of conic cols
    VA::Vector{Float64} # List of conic coefficients
    b::Vector{Float64} # List of conic coefficients
    IG::Vector{Int} # List of equality rows
    JG::Vector{Int} # List of equality cols
    VG::Vector{Float64} # List of equality coefficients
    h::Vector{Float64} # List of equality coefficients
    objective_constant::Float64 # The objective is min c'x + objective_constant
    c::Vector{Float64}
end

# This is tied to ECOS's internal representation
mutable struct ConeData
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    ep::Int # number of primal exponential cone triples
    # The following four field store model information needed to compute `ConstraintPrimal` and `ConstraintDual`
    eqnrows::Dict{Int, Int}   # The number of rows of Zeros
    ineqnrows::Dict{Int, Int} # The number of rows of each vector sets except Zeros
    function ConeData()
        new(0, 0, 0, Int[], 0,
            Dict{Int, UnitRange{Int}}(),
            Dict{Int, UnitRange{Int}}())
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cone::ConeData
    maxsense::Bool
    data::Union{Nothing, ModelData} # only non-Nothing between MOI.copy_to and MOI.optimize!
    sol::Solution
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
        optimizer = new(ConeData(), false, nothing, Solution(), false, Dict{Symbol, Any}())
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawParameter(key), value)
        end
        return optimizer
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "ECOS"

function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    optimizer.options[param.name] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    # TODO: This gives a poor error message if the name of the parameter is invalid.
    return optimizer.options[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

function MOI.is_empty(optimizer::Optimizer)
    !optimizer.maxsense && optimizer.data === nothing
end

function MOI.empty!(optimizer::Optimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    optimizer.sol = Solution()
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.SingleVariable},
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

function MOI.supports_constraint(::Optimizer,
                                 ::Type{MOI.VectorAffineFunction{Float64}},
                                 ::Type{<:Union{MOI.Zeros, MOI.Nonnegatives,
                                                MOI.SecondOrderCone,
                                                MOI.ExponentialCone}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

using SparseArrays

# Computes cone dimensions
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Zeros}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.Zeros)
    ci = cone.f
    cone.f += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Nonnegatives}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.Nonnegatives)
    ci = cone.l
    cone.l += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:MOI.SecondOrderCone}) = cone.l + ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:MOI.ExponentialCone}) = cone.l + cone.q + ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.ExponentialCone)
    ci = 3cone.ep
    cone.ep += 1
    ci
end
constroffset(optimizer::Optimizer, ci::CI) = constroffset(optimizer.cone, ci::CI)
function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
end

# Build constraint matrix
output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, MOI.Zeros}) = 1:optimizer.cone.eqnrows[constroffset(optimizer, ci)]
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.ineqnrows[constroffset(optimizer, ci)]
matrix(data::ModelData, s::MOI.Zeros) = data.b, data.IA, data.JA, data.VA
matrix(data::ModelData, s::Union{MOI.Nonnegatives, MOI.SecondOrderCone, MOI.ExponentialCone}) = data.h, data.IG, data.JG, data.VG
matrix(optimizer::Optimizer, s) = matrix(optimizer.data, s)
# ECOS orders differently than MOI the second and third dimension of the exponential cone
orderval(val, s) = val
function orderval(val, s::Union{MOI.ExponentialCone, Type{MOI.ExponentialCone}})
    val[[1, 3, 2]]
end
orderidx(idx, s) = idx
expmap(i) = (1, 3, 2)[i]
function orderidx(idx, s::MOI.ExponentialCone)
    expmap.(idx)
end
function MOIU.load_constraint(optimizer::Optimizer, ci::MOI.ConstraintIndex, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    I, J, V = findnz(A)
    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    if s isa MOI.Zeros
        optimizer.cone.eqnrows[offset] = length(rows)
    else
        optimizer.cone.ineqnrows[offset] = length(rows)
    end
    i = offset .+ rows
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    b, Is, Js, Vs = matrix(optimizer, s)
    b[i] .= orderval(f.constants, s)
    append!(Is, offset .+ orderidx(I, s))
    append!(Js, J)
    append!(Vs, -V)
end

function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
    optimizer.cone = ConeData()
    VI.(1:nvars)
end

function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.l + cone.q + 3cone.ep
    IA = Int[]
    JA = Int[]
    VA = Float64[]
    b = zeros(cone.f)
    IG = Int[]
    JG = Int[]
    VG = Float64[]
    h = zeros(m)
    c = zeros(nvars)
    optimizer.data = ModelData(m, nvars, IA, JA, VA, b, IG, JG, VG, h, 0., c)
end

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MAX_SENSE
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.Union{MOI.SingleVariable,
                                   MOI.ScalarAffineFunction{Float64}})
end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.SingleVariable)
    MOIU.load(optimizer,
              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
              MOI.ScalarAffineFunction{Float64}(f))
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          optimizer.data.n))
    optimizer.data.objective_constant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
    return nothing
end

function MOI.optimize!(optimizer::Optimizer)
    if optimizer.data === nothing
        # optimize! has already been called and no new model has been copied
        return
    end
    cone = optimizer.cone
    m = optimizer.data.m
    n = optimizer.data.n
    A = ECOS.ECOSMatrix(sparse(optimizer.data.IA, optimizer.data.JA, optimizer.data.VA, cone.f, n))
    b = optimizer.data.b
    G = ECOS.ECOSMatrix(sparse(optimizer.data.IG, optimizer.data.JG, optimizer.data.VG, m, n))
    h = optimizer.data.h
    objective_constant = optimizer.data.objective_constant
    c = optimizer.data.c
    optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded to ECOS
    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:verbose] = false
    end
    ecos_prob_ptr = ECOS.setup(n, m, cone.f, cone.l, length(cone.qa), cone.qa,
                               cone.ep, G, A, c, h, b; options...)
    ret_val = ECOS.solve(ecos_prob_ptr)
    stat = unsafe_load(unsafe_load(ecos_prob_ptr).info)
    solve_time = stat.tsetup + stat.tsolve
    ecos_prob = unsafe_wrap(Array, ecos_prob_ptr, 1)[1]
    primal    = unsafe_wrap(Array, ecos_prob.x, n)[:]
    dual_eq   = unsafe_wrap(Array, ecos_prob.y, cone.f)[:]
    dual_ineq = unsafe_wrap(Array, ecos_prob.z, m)[:]
    slack     = unsafe_wrap(Array, ecos_prob.s, m)[:]
    ECOS.cleanup(ecos_prob_ptr, 0)
    objective_value = (optimizer.maxsense ? -1 : 1) * stat.pcost
    dual_objective_value = (optimizer.maxsense ? -1 : 1) * stat.dcost
    optimizer.sol = Solution(ret_val, primal, dual_eq, dual_ineq, slack, objective_value,
                            dual_objective_value, objective_constant, solve_time)
end

MOI.get(optimizer::Optimizer, ::MOI.SolveTime) = optimizer.sol.solve_time
function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    # Strings from https://github.com/ifa-ethz/ecos/blob/master/include/ecos.h
    flag = optimizer.sol.ret_val
    if flag == OPTIMIZE_NOT_CALLED
        return "Optimize not called"
    elseif flag == ECOS_OPTIMAL
        return "Problem solved to optimality"
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return "Problem solved to inaccurate optimality"
    elseif flag == ECOS_PINF
        return "Found certificate of primal infeasibility"
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of primal infeasibility"
    elseif flag == ECOS_DINF
        return "Found certificate of dual infeasibility"
    elseif flag == ECOS_DINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of dual infeasibility"
    elseif flag == ECOS_MAXIT
        return "Maximum number of iterations reached"
    elseif flag == ECOS_NUMERICS
        return "Search direction unreliable"
    elseif flag == ECOS_OUTCONE
        return "s or z got outside the cone, numerics?"
    elseif flag == ECOS_SIGINT
        return "solver interrupted by a signal/ctrl-c"
    else
        @assert flag == ECOS_FATAL
        return "Unknown problem in solver"
    end
end

# Implements getter for result value and statuses
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    flag = optimizer.sol.ret_val
    if flag == OPTIMIZE_NOT_CALLED
        return MOI.OPTIMIZE_NOT_CALLED
    elseif flag == ECOS.ECOS_OPTIMAL
        return MOI.OPTIMAL
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBLE
    elseif flag == ECOS.ECOS_DINF
        return MOI.DUAL_INFEASIBLE
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.IterationLimit
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_OPTIMAL
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_INFEASIBLE
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_DUAL_INFEASIBLE
    else
        return MOI.OTHER_ERROR
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.objective_value
    if !MOIU.is_ray(MOI.get(optimizer, MOI.PrimalStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end
function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.dual_objective_value
    if !MOIU.is_ray(MOI.get(optimizer, MOI.DualStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_DINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    else
        return MOI.OTHER_RESULT_STATUS
    end
end
# Swapping indices 2 <-> 3 is an involution (it is its own inverse)
const reorderval = orderval
function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(optimizer, attr)
    optimizer.sol.primal[vi.value]
end
MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI}) = MOI.get.(optimizer, Ref(a), vi)
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    MOI.check_result_index_bounds(optimizer, attr)
    rows = constrrows(optimizer, ci)
    return zeros(length(rows))
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, attr)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    return reorderval(optimizer.sol.slack[offset .+ rows], S)
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_DINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    else
        return MOI.OTHER_RESULT_STATUS
    end
end
_dual(optimizer, ci::CI{<:MOI.AbstractFunction, <:MOI.Zeros}) = optimizer.sol.dual_eq
_dual(optimizer, ci::CI) = optimizer.sol.dual_ineq
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, attr)
    offset = constroffset(optimizer, ci)
    rows = constrrows(optimizer, ci)
    return reorderval(_dual(optimizer, ci)[offset .+ rows], S)
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1


text:
using DoubleDouble
using Base.Test

x = sqrt(2.0)
bx = big(x)
sx = Single(x)
dx = Double(x)

y = 0.1
by = big(y)
sy = Single(y)
dy = Double(y)

@test x == sx == dx
@test y == sy == dy

@test zero(Double{Float64}) == Double(0.0, 0.0)
@test one(Double{Float64}) == Double(1.0, 0.0)
@test zero(Double(0.0, 0.0)) == Double(0.0, 0.0)
@test one(Double(0.0, 0.0)) == Double(1.0, 0.0)

dxy = dx*dy
bxy = bx*by
@test sx*sy == dxy
@test x*y == Float64(dxy)
@test dxy == Double(bxy)

@test x+y == Float64(dx+dy)
@test dx+dy == Double(bx+by)

@test x-y == Float64(dx-dy)
@test dx-dy == Double(bx-by)

@test x/y == Float64(dx/dy)
@test dx/dy == Double(bx/by)

@test sqrt(y) == Float64(sqrt(dy))
@test sqrt(dy) == Double(sqrt(by))

#@test rem(dxy,1.0) == Double(rem(bxy,1.0))

## New
@test Double(pi) == Double{Float64}(3.141592653589793, 1.2246467991473532e-16)
@test Double(3.5) == Double{Float64}(3.5, 0.0)
@test Double(3.5) == Double{Float64}(3.5, 0.0)

a = Double(big"3.1")
@test a == Double(3.1, -8.881784197001253e-17)
@test BigFloat(a) == big"3.099999999999999999999999999999995069619342368676216176696466982586064542459781"

@test Single(3) === Single(3.0)
@test Double(3) === Double(3.0, 0.0)
@test Double(big(3)) === Double(3.0, 0.0)

@test convert(Single{Float64}, 1) === Single(1.0)
@test convert(Single{Float32}, 1) === Single(1.0f0)
@test convert(Double{Float64}, 1) === Double(1.0, 0.0)
@test convert(Double{Float32}, 1) === Double(1.0f0, 0.0f0)

@test Single{Float32}(3) === Single{Float32}(3.0f0)
@test Double{Float32}(3) === Double{Float32}(3.0f0, 0.0f0)
@test Single{Float32}(BigFloat(3)) === Single{Float32}(3.0f0)
@test Double{Float32}(BigFloat(3)) === Double{Float32}(3.0f0, 0.0f0)


text:
export make_node

@node typeof(dot) Deterministic [ out, in1, in2 ]

# By default dot-product node uses TinyCorrection() strategy for precision matrix on `in1` and `in2` edges to ensure precision is always invertible
default_meta(::typeof(dot)) = TinyCorrection()


text:
NIScope_Triggering_Dict = Dict()
# Set trigger type
NIScope_Triggering_Dict["edge"]      = NISCOPE_VAL_EDGE_TRIGGER
NIScope_Triggering_Dict["width"]     = NISCOPE_VAL_WIDTH_TRIGGER
NIScope_Triggering_Dict["runt"]      = NISCOPE_VAL_RUNT_TRIGGER
NIScope_Triggering_Dict["glitch"]    = NISCOPE_VAL_GLITCH_TRIGGER
NIScope_Triggering_Dict["tv"]        = NISCOPE_VAL_TV_TRIGGER
NIScope_Triggering_Dict["immediate"] = NISCOPE_VAL_IMMEDIATE_TRIGGER
NIScope_Triggering_Dict["acline"]    = 7
# Set trigger source, niScope channels start from 0, so subtract 1
NIScope_Triggering_Dict["1"]   = "0"
NIScope_Triggering_Dict["2"]   = "1"
NIScope_Triggering_Dict["3"]   = "2"
NIScope_Triggering_Dict["4"]   = "3"
NIScope_Triggering_Dict["5"]   = "4"
NIScope_Triggering_Dict["6"]   = "5"
NIScope_Triggering_Dict["7"]   = "6"
NIScope_Triggering_Dict["8"]   = "7"
NIScope_Triggering_Dict["ext"] = NISCOPE_VAL_EXTERNAL
NIScope_Triggering_Dict["ttl0"]= NISCOPE_VAL_TTL0
# Set trigger slope
NIScope_Triggering_Dict["pos"]    = NISCOPE_VAL_POSITIVE
NIScope_Triggering_Dict["neg"]    = IVISCOPE_VAL_NEGATIVE
NIScope_Triggering_Dict["either"] = NISCOPE_VAL_SLOPE_EITHER
# Trigger coupling
NIScope_Triggering_Dict["ac"]    = NISCOPE_VAL_AC
NIScope_Triggering_Dict["dc"]    = NISCOPE_VAL_DC
NIScope_Triggering_Dict["hfrej"] = NISCOPE_VAL_HF_REJECT
NIScope_Triggering_Dict["lfrej"] = NISCOPE_VAL_LF_REJECT
NIScope_Triggering_Dict["noiserej"] = 5#??
# Ext Trigger impedance
NIScope_Triggering_Dict["50"]    = NISCOPE_VAL_50_OHMS
NIScope_Triggering_Dict["1M"]    = NISCOPE_VAL_1_MEG_OHM

# Set trigger type
function set_trig_type(obj::T;tp::String="edge") where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_TYPE=1250012; # ((1000000 + 250000) + 12)
	tp = NIScope_Triggering_Dict[tp]
	status = niScope_SetAttributeViInt32(obj, "", NISCOPE_ATTR_TRIGGER_TYPE, Int32(tp))
end
# Set trigger source
function set_trig_src(obj::T;src::String="1") where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_SOURCE=1250013; # ((1000000 + 250000) + 13)
	src = NIScope_Triggering_Dict[src]
	status = niScope_SetAttributeViString(obj, "", NISCOPE_ATTR_TRIGGER_SOURCE, src)
end
# Set trigger mode
function set_trig_mode(obj::T;mode::String="normal") where {T<:NIScope}
	@warn "not available for instruments of type ::NIScope; default to SINGLE"
end
# Set trigger slope
function set_trig_slope(obj::T;slp::String="pos") where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_SLOPE=1250018; # ((1000000 + 250000) + 18)
	slp = NIScope_Triggering_Dict[slp]
	status = niScope_SetAttributeViInt32(obj, "", NISCOPE_ATTR_TRIGGER_SLOPE, Int32(slp))
end
# Set trigger level
function set_trig_lev(obj::T;lev=1) where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_LEVEL=1250017; # ((1000000 + 250000) + 17)
	status = niScope_SetAttributeViReal64(obj, "", NISCOPE_ATTR_TRIGGER_LEVEL, Float64(lev))
end
# Set trigger holdoff
function set_trig_holdoff(obj::T;hf=1e-9) where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_HOLDOFF=1250016; # ((1000000 + 250000) + 16)
	status = niScope_SetAttributeViReal64(obj, "", NISCOPE_ATTR_TRIGGER_HOLDOFF, Float64(hf))
end
# Trigger delay time
function set_hoffs_t(obj::T;href=1e-3) where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_DELAY_TIME=1250015; # ((1000000 + 250000) + 15)
	status = niScope_SetAttributeViReal64(obj, "", NISCOPE_ATTR_TRIGGER_DELAY_TIME, Float64(href))
end
# Set trigger coupling
function set_trig_cpl(obj::T;ch::String="1",cpl::String="dc") where {T<:NIScope}
	#NISCOPE_ATTR_TRIGGER_COUPLING=1250014; # ((1000000 + 250000) + 14)
	cpl = NIScope_Triggering_Dict[cpl]
	status = niScope_SetAttributeViInt32(obj, "", NISCOPE_ATTR_TRIGGER_COUPLING, Int32(cpl))
end


text:
module LearningTrees
    using Random
    using DataFrames, CSV, StatsBase

    export
        traintestsplit,
        gini,
        entropy,
        confusionmatrix,
        accuracyscore,
        recallscore,
        precisionscore,
        f1score,

        DecisionTree,
        RandomForest,
        fit!,
        predict

    abstract type Classifier end

    include("utils.jl")
    include("tree.jl")
    include("decisiontree.jl")
    include("randomforest.jl")
end


text:
using Omega
using Test

function testcond()
  x = normal(0.0, 1.0)
  x_ = cond(x, x < 1.0)
  x_samples = rand(x_, 10000, alg = RejectionSample)
  x_samples = filter(s -> s != nothing, x_samples)
  @test maximum(x_samples) < 1.0
  x__ = cond(x_, x > -1.0)
  x__samples = rand(x__, 10000, alg = RejectionSample)
  x__samples = filter(s -> s != nothing, x__samples)
  @test minimum(x__samples) > -1.0
end

testcond()

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
# Basic use:

    VI(logl, [gradlogl], μ, [Σ]; S = 100, iterations = 1, show_every = -1)

Returns approximate Gaussian posterior and log evidence

# Arguments

- `logl` is a function that expresses the joint log-likelihood
- `gradlogl` calculates the gradient for `logl`.
  If not given, it will be obtained via automatic differentiation.
- `μ` is the mean of the initial Gaussian posterior provided as an array.
  If instead, an array of `μ` is passed, they will all be tried out as initial
  candidate solutions and the best one, according to the lower bound, will be selected.
- `Σ` is the covariance of the initial Gaussian posterior.
  If unsecified, its  default values is `0.1*I`.
  If instead, an array of `Σ` is passed (assuming also an array of `μ` is passed), they will all be tried out as initial
  candidate solutions and the best one, according to the lower bound, will be selected.
- `S` is the number of drawn samples that approximate the lower bound integral.
- `Stest` is the number of drawn samples to test the lower bound integral for overfitting.
- `show_every`: report progress every `show_every` number of iterations.

More arguments are explained in the README.md file.

# Example

```julia-repl
# infer posterior of Bayesian linear regression, compare to exact result
julia> using LinearAlgebra, Distributions
julia> D = 4; X = randn(D, 1000); W = randn(D); β = 0.3; α = 1.0;
julia> Y = vec(W'*X); Y += randn(size(Y))/sqrt(β);
julia> Sn = inv(α*I + β*(X*X')) ; mn = β*Sn*X*Y; # exact posterior
julia> posterior, logev = VI( w -> logpdf(MvNormal(vec(w'*X), sqrt(1/β)), Y) + logpdf(MvNormal(zeros(D),sqrt(1/α)), w), randn(D); S = 1_000, iterations = 15);
julia> display([mean(posterior) mn])
julia> display([cov(posterior)  Sn])
julia> display(logev) # display negative log evidence
```

"""
function VI(logl::Function, μ::Array{Float64,1}, Σ = Matrix(0.1*I, length(μ), length(μ)); gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    coreVIfull(logl, [μ], [Σ]; gradlogl = gradlogl, seed = seed, S = S, optimiser=optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=0)

end

function VI(logl::Function, initgaussian::MvNormal; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    VI(logl, mean(initgaussian), cov(initgaussian); gradlogl = gradlogl, seed = seed, S = S, optimiser=optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=0)

end

function VI(logl::Function, μ::Array{Array{Float64,1},1}, Σ = [Matrix(0.1*I, length(μ[1]), length(μ[1])) for _ in 1:length(μ)]; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    @assert(length(μ) == length(Σ))

    coreVIfull(logl, μ, Σ; gradlogl = gradlogl, seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end

#-----------------------------------#
# Call mean field                   #
#-----------------------------------#

function VIdiag(logl::Function, μ::Array{Float64,1}, Σdiag = 0.1*ones(length(μ)); gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    coreVIdiag(logl, [μ], [Σdiag]; gradlogl = gradlogl, seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end

function VIdiag(logl::Function, initgaussian::MvNormal; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    VIdiag(logl, mean(initgaussian), diag(cov(initgaussian)); gradlogl = gradlogl, seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end

function VIdiag(logl::Function, μ::Array{Array{Float64,1},1}, Σdiag = [0.1*ones(length(μ[1])) for _ in 1:length(μ)]; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100,  iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    coreVIdiag(logl, μ, Σdiag; gradlogl = gradlogl, seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end

#-----------------------------------#
# Call VI with spherical covariance #
#-----------------------------------#

function VIfixedcov(logl::Function, μ::Array{Float64,1}, fixedC::Array{Float64,2}; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0, adaptvariance = 1)

    coreVIfixedcov(logl, μ, fixedC, gradlogl = gradlogl, seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations = inititerations, adaptvariance = adaptvariance)

end

#-----------------------------------#
# Call Laplace                      #
#-----------------------------------#

function laplace(logl::Function, x::Array{Float64,1}; gradlogl = x -> ForwardDiff.gradient(logl, x), hesslog = x->ForwardDiff.hessian(logl, x), optimiser=Optim.LBFGS(), iterations=1000, show_every=-1)

    laplace(logl, [x]; gradlogl=gradlogl, hesslog = hesslog, optimiser = optimiser, iterations = iterations, show_every = show_every)

end

function laplace(logl::Function, X::Array{Array{Float64,1},1}; gradlogl = x -> ForwardDiff.gradient(logl, x), hesslog = x->ForwardDiff.hessian(logl, x), optimiser=Optim.LBFGS(), iterations=1000, show_every=-1)

    map(x->coreLaplace(logl, gradlogl, hesslog, x; iterations = iterations, optimiser = optimiser, show_every = show_every), X)

end

#-----------------------------------#
# Call Mixed Variational Inference  #
#-----------------------------------#

function MVI(logl::Function, μ::Array{Float64,1}; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), laplaceiterations=10_000,  seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    MVI(logl, [μ]; gradlogl = gradlogl, seed = seed, S = S, optimiser=optimiser, laplaceiterations=laplaceiterations, iterations=iterations, numerical_verification = numerical_verification, Stest=Stest, show_every=show_every, inititerations=inititerations)

end

function MVI(logl::Function, μ::Array{Array{Float64,1},1}; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), laplaceiterations=10_000,  seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    LAposteriors = laplace(logl, μ; gradlogl = gradlogl, optimiser=optimiser, iterations=laplaceiterations, show_every=show_every)

    coreMVI(logl, gradlogl, LAposteriors; seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end

function MVI(logl::Function, LAposterior::MvNormal; gradlogl = x -> ForwardDiff.gradient(logl, x), optimiser=Optim.LBFGS(), laplaceiterations=10_000, seed = 1, S = 100, iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    coreMVI(logl, gradlogl, [LAposterior]; seed = seed, S = S, optimiser = optimiser, iterations = iterations, numerical_verification = numerical_verification, Stest = Stest, show_every = show_every, inititerations=inititerations)

end


text:
include("$(@__DIR__)/../src/FiniteGroups.jl")
using .FiniteGroups, LinearAlgebra, Test

g = pointgroup("Oh")
rreps = irreps(g, R=true)

rep = oplus(
    rreps[1],   # 1
    rreps[2],   # 2
    rreps[2],   # 3
    rreps[2],   # 4
    rreps[2],   # 5
    rreps[5],   # 6
    rreps[6],   # 7
    rreps[8],   # 8
    rreps[8],   # 9
    rreps[9],   # 10
    rreps[10]   # 11
)
ru = Symmetric(rand(size(rep[1])...)) |> eigvecs
rep = transform_rep(ru', rep, ru)
res = block_decomposition(rep, g)

@testset "Oh Decomp" begin
    @test length(res) == 11
    @test norm(transform_rep(res[1]', rep, res[1]) .- rreps[1]) < 1e-10
    for i = 2:5
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[2]) < 1e-10
    end
    for i = 6:6
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[5]) < 1e-10
    end
    for i = 7:7
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[6]) < 1e-10
    end
    for i = 8:9
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[8]) < 1e-10
    end
    for i = 10:10
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[9]) < 1e-10
    end
    for i = 11:11
        @test norm(transform_rep(res[i]', rep, res[i]) .- rreps[10]) < 1e-10
    end
end


text:
function score(x, y)

end