text:
function Base.iszero(x::Float64)
    x ≈ 0.0
end

abstract type AbstractSparseM{Tv, Ti} <: AbstractMatrix{Tv} end

struct SparseCSR{Tv,Ti} <: AbstractSparseM{Tv,Ti}
    m::Ti
    n::Ti
    val::Vector{Tv}
    rowptr::Vector{Ti}
    colind::Vector{Ti}
end

struct SparseCSC{Tv,Ti} <: AbstractSparseM{Tv,Ti}
    m::Ti
    n::Ti
    val::Vector{Tv}
    colptr::Vector{Ti}
    rowind::Vector{Ti}
end

struct SparseCOO{Tv,Ti} <: AbstractSparseM{Tv,Ti}
    m::Ti
    n::Ti
    val::Vector{Tv}
    rowind::Vector{Ti}
    colind::Vector{Ti}
end

SparseCSR(A::Matrix{Tv}) where {Tv} = _tocsr(A, Int)
SparseCSR(A::SparseCSR{Tv,Ti}) where {Tv,Ti} = copy(A)
SparseCSR(A::SparseCSC{Tv,Ti}) where {Tv,Ti} = _tocsr(_tocoo(A))
SparseCSR(A::SparseCOO{Tv,Ti}) where {Tv,Ti} = _tocsr(A)
SparseCSR(A::SparseArrays.SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = _tocsr(_tocoo(_tocsc(A)))

SparseCSC(A::Matrix{Tv}) where {Tv} = _tocsc(A, Int)
SparseCSC(A::SparseCSR{Tv,Ti}) where {Tv,Ti} = _tocsc(_tocoo(A))
SparseCSC(A::SparseCSC{Tv,Ti}) where {Tv,Ti} = copy(A)
SparseCSC(A::SparseCOO{Tv,Ti}) where {Tv,Ti} = _tocsc(A)
SparseCSC(A::SparseArrays.SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = _tocsc(A)

SparseCOO(A::Matrix{Tv}) where {Tv} = _tocoo(A, Int)
SparseCOO(A::SparseCSR{Tv,Ti}) where {Tv,Ti} = _tocoo(A)
SparseCOO(A::SparseCSC{Tv,Ti}) where {Tv,Ti} = _tocoo(A)
SparseCOO(A::SparseCOO{Tv,Ti}) where {Tv,Ti} = copy(A)
SparseCOO(A::SparseArrays.SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = _tocoo(_tocsc(A))

Matrix(A::SparseCSR{Tv,Ti}) where {Tv,Ti} = _dense(A)
Matrix(A::SparseCSC{Tv,Ti}) where {Tv,Ti} = _dense(A)
Matrix(A::SparseCOO{Tv,Ti}) where {Tv,Ti} = _dense(A)

SparseArrays.sparse(A::SparseCSR{Tv,Ti}) where {Tv,Ti} = SparseArrays.sparse(_tocsc(_tocoo(A)))
SparseArrays.sparse(A::SparseCSC{Tv,Ti}) where {Tv,Ti} = SparseArrays.SparseMatrixCSC{Tv,Ti}(A.m, A.n, copy(A.colptr), copy(A.rowind), copy(A.val))
SparseArrays.sparse(A::SparseCOO{Tv,Ti}) where {Tv,Ti} = SparseArrays.sparse(_tocsc(A))

function SparseCOO(m::Ti, n::Ti, elem::AbstractArray{Tuple{Ti,Ti,Tv},1}) where {Tv,Ti}
    rowind = Vector{Ti}()
    colind = Vector{Ti}()
    val = Vector{Tv}()
    prev_index::Tuple{Ti,Ti} = (0,0)
    for (i,j,u) in sort(elem)
        if m < i
            m = i
        end
        if n < j
            n = j
        end
        if prev_index != (i,j)
            push!(rowind, i)
            push!(colind, j)
            push!(val, u)
            prev_index = (i,j)
        else
            val[end] += u
        end
    end
    SparseCOO(m, n, val, rowind, colind)
end

### overload

function Base.show(io::IO, A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(m)*"x"*string(n)*" CSR-SparseMatrix")
    for i = 1:m
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
        end
    end
end

function Base.show(io::IO, A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(A.m)*"x"*string(A.n)*" CSC-SparseMatrix")
    for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
        end
    end
end

function Base.show(io::IO, A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(A.m)*"x"*string(A.n)*" COO-SparseMatrix")
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
    end
end

function Base.show(io::IO, ::MIME"text/plain", A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(m)*"x"*string(n)*" CSR-SparseMatrix")
    for i = 1:m
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(A.m)*"x"*string(A.n)*" CSC-SparseMatrix")
    for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(A.m)*"x"*string(A.n)*" COO-SparseMatrix")
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        println(io, "("*string(i)*","*string(j)*") "*string(A.val[z]))
    end
end

function SparseArrays.nnz(A::AbstractSparseM{Tv,Ti}) where {Tv,Ti}
    return length(A.val)
end

####

function Base.size(A::AbstractSparseM{Tv,Ti}) where {Tv,Ti}
    return (A.m, A.n)
end

function Base.eachindex(A::AbstractSparseM{Tv,Ti}) where {Tv,Ti}
    return eachindex(A.val)
end

function Base.length(A::AbstractSparseM{Tv,Ti}) where {Tv,Ti}
    return length(A.val)
end

function Base.getindex(A::AbstractSparseM{Tv,Ti}, i::Ti) where {Tv,Ti}
    A.val[i]
end

function Base.setindex!(A::AbstractSparseM{Tv,Ti}, value::Tv, i::Ti) where {Tv,Ti}
    A.val[i] = value
end

# function Base.iterate(A::AbstractSparseM{Tv,Ti}, i::Ti = 1) where {Tv,Ti}
#     i == length(A)+1 && return nothing
#     return (A.val[i], i+1)
# end

####

function Base.copy(A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    SparseCSR(A.m, A.n, copy(A.val), A.rowptr, A.colind)
end

function Base.copy(A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    SparseCSC(A.m, A.n, copy(A.val), A.colptr, A.rowind)
end

function Base.copy(A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    SparseCOO(A.m, A.n, copy(A.val), A.rowind, A.colind)
end

####

function Base.similar(A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    SparseCSR(A.m, A.n, similar(A.val), A.rowptr, A.colind)
end

function Base.similar(A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    SparseCSC(A.m, A.n, similar(A.val), A.colptr, A.rowind)
end

function Base.similar(A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    SparseCOO(A.m, A.n, similar(A.val), A.rowind, A.colind)
end

####

function Base.zero(A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    SparseCSR(A.m, A.n, zero(A.val), A.rowptr, A.colind)
end

function Base.zero(A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    SparseCSC(A.m, A.n, zero(A.val), A.colptr, A.rowind)
end

function Base.zero(A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    SparseCOO(A.m, A.n, zero(A.val), A.rowind, A.colind)
end

function Base.zero(A::SparseArrays.SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    SparseArrays.SparseMatrixCSC{Tv,Ti}(A.m, A.n, A.colptr, A.rowval, zero(A.nzval))
end

function _tocsr(A::Matrix{Tv}, ::Type{Ti})::SparseCSR{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    rowptr = Vector{Ti}(undef, m+1)
    colind = Vector{Ti}()
    val = Vector{Tv}()
    rowptr[1] = 1
    for i = 1:m
        for j = 1:n
            if !(iszero(A[i,j]))
                push!(colind, j)
                push!(val, A[i,j])
            end
        end
        rowptr[i+1] = length(val)+1
    end
    SparseCSR(m, n, val, rowptr, colind)
end

function _tocsr(A::SparseCOO{Tv,Ti})::SparseCSR{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    rowptr = Vector{Ti}(undef, m+1)
    colind = Vector{Ti}()
    val = Vector{Tv}()
    p = 1
    i = 1
    rowptr[i] = p
    for x in sort(collect(zip(A.rowind, A.colind, A.val)))
        if i != x[1]
            for u = i+1:x[1]
                rowptr[u] = p
            end
            i = x[1]
        end
        push!(colind, x[2])
        push!(val, x[3])
        p += 1
    end
    for u = i+1:m+1
        rowptr[u] = p
    end
    SparseCSR(m, n, val, rowptr, colind)
end

function _tocsc(A::Matrix{Tv}, ::Type{Ti})::SparseCSC{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    colptr = Vector{Ti}(undef, n+1)
    rowind = Vector{Ti}()
    val = Vector{Tv}()
    colptr[1] = 1
    for j = 1:n
        for i = 1:m
            if !(iszero(A[i,j]))
                push!(rowind, i)
                push!(val, A[i,j])
            end
        end
        colptr[j+1] = length(val)+1
    end
    SparseCSC(m, n, val, colptr, rowind)
end

function _tocsc(A::SparseCOO{Tv,Ti})::SparseCSC{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    colptr = Vector{Ti}(undef, n+1)
    rowind = Vector{Ti}()
    val = Vector{Tv}()
    p = 1
    j = 1
    colptr[j] = p
    for x in sort(collect(zip(A.colind, A.rowind, A.val)))
        if j != x[1]
            for u = j+1:x[1]
                colptr[u] = p
            end
            j = x[1]
        end
        push!(rowind, x[2])
        push!(val, x[3])
        p += 1
    end
    for u = j+1:n+1
        colptr[u] = p
    end
    SparseCSC(m, n, val, colptr, rowind)
end

function _tocsc(A::SparseArrays.SparseMatrixCSC{Tv,Ti})::SparseCSC{Tv,Ti} where {Tv, Ti}
    SparseCSC(A.m, A.n, copy(A.nzval), copy(A.colptr), copy(A.rowval))
end

function _tocoo(A::Matrix{Tv}, ::Type{Ti})::SparseCOO{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    rowind = Vector{Ti}()
    colind = Vector{Ti}()
    val = Vector{Tv}()
    for j = 1:n
        for i = 1:m
            if !(iszero(A[i,j]))
                push!(rowind, i)
                push!(colind, j)
                push!(val, A[i,j])
            end
        end
    end
    SparseCOO(m, n, val, rowind, colind)
end

function _tocoo(A::SparseCSR{Tv,Ti})::SparseCOO{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    rowind = Vector{Ti}()
    colind = Vector{Ti}()
    val = Vector{Tv}()
    for i = 1:m
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            push!(rowind, i)
            push!(colind, j)
            push!(val, A.val[z])
        end
    end
    SparseCOO(m, n, val, rowind, colind)
end

function _tocoo(A::SparseCSC{Tv,Ti})::SparseCOO{Tv,Ti} where {Tv, Ti}
    m, n = size(A)
    rowind = Vector{Ti}()
    colind = Vector{Ti}()
    val = Vector{Tv}()
    for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            push!(rowind, i)
            push!(colind, j)
            push!(val, A.val[z])
        end
    end
    SparseCOO(m, n, val, rowind, colind)
end

function _dense(A::SparseCSR{Tv,Ti})::Matrix{Tv} where {Tv, Ti}
    m, n = size(A)
    M = zeros(m,n)
    for i = 1:m
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            M[i,j] = A.val[z]
        end
    end
    return M
end

function _dense(A::SparseCSC{Tv,Ti})::Matrix{Tv} where {Tv, Ti}
    m, n = size(A)
    M = zeros(m,n)
    for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            M[i,j] = A.val[z]
        end
    end
    return M
end

function _dense(A::SparseCOO{Tv,Ti})::Matrix{Tv} where {Tv, Ti}
    m, n = size(A)
    M = zeros(m,n)
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        M[i,j] = A.val[z]
    end
    return M
end

#################### Adjoint

struct Adjoint{T <: AbstractSparseM}
    parent::T
end

function Base.adjoint(A::AbstractSparseM{Tv,Ti}) where {Tv,Ti}
    Adjoint(A)
end

function Base.adjoint(A::Adjoint{T}) where {T}
    A.parent
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for the Vertex DataFrame.

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

@testset "Edge DataFrame" begin
   g = completegraph(10)

   eit = edges(g)

   float_arr = rand(90)
   int_arr = rand(Int, 90)
   bool_arr = rand(Bool, 90)
   str_arr = map(string, 1:90)
   void_arr = fill(nothing, 90)

   E = g.edata

   # Mutation on emtpy DataFrame
   addedge!(E)
   @test size(E) == (0,0)
   rmedge!(E, 11)
   @test size(E) == (0,0)
   rmedge!(E, 1:90)
   @test size(E) == (0,0)

   ###
   # SETVPROP
   ###
   # Try to create a property only on a subset of properties
   try seteprop!(g, eit[rand(1:90)], 1, :p1); @test false catch @test true end
   try seteprop!(g, eit[rand(1:90, rand(1:90))], 1, :p1); @test false catch @test true end
   try seteprop!(g, :, rand(rand(1:89)), :p1); @test false catch @test true end

   # Create new properties
   seteprop!(g, :, float_arr, :p1)
   seteprop!(g, :, int_arr, :p2)
   seteprop!(g, :, bool_arr, :p3)
   seteprop!(g, :, str_arr, :p4)
   seteprop!(g, :, nothing, :p5)

   # Test property names and types
   @test listeprops(g) == [:p1, :p2, :p3, :p4, :p5]
   @test eltypes(g.edata) == [Float64, Int, Bool, String, Void]

   # Test property values
   @test geteprop(g, :, :p1) == float_arr
   @test geteprop(g, :, :p2) == int_arr
   @test geteprop(g, :, :p3) == bool_arr
   @test geteprop(g, :, :p4) == str_arr
   @test geteprop(g, :, :p5) == void_arr

   # Mutation on DataFrame
   @test addedge!(g, 1=>1) == true
   @test g.indxs[1=>1] == 91
   for eprop in listeprops(g)
      isequal(geteprop(g, 1=>1, eprop), NA)
   end

   rmedge!(g, 1=>1)
   @test size(g.edata) == (90,5)

   # Modify single edge
   e = rand(eit)

   seteprop!(g, e, 0.5, :p1)
   @test geteprop(g, e, :p1) == 0.5

   seteprop!(g, e, "p4", :p4)
   @test geteprop(g, e, :p4) == "p4"

   # Modify a range of edges
   es = rand(eit, 20)
   seteprop!(g, es, 0, :p1)
   @test geteprop(g, es, :p1) == zeros(20)

   es = unique(rand(eit, 50))
   vals = 1 : length(es)
   seteprop!(g, es, vals, :p2)
   @test geteprop(g, es, :p2) == collect(vals)
end

ext: jl


text:
using ChainRulesTestUtils, DifferenceEquations, Distributions, LinearAlgebra, Test, Zygote
using DelimitedFiles
using DiffEqBase
using FiniteDiff: finite_difference_gradient

# Matrices from RBC
A_0_rbc = [-7.824904812740593e-5, 0.0]
A_1_rbc = [0.9568351489231076 6.209371005755285; 3.0153731819288737e-18 0.20000000000000007]
A_2_rbc = cat([-0.00019761505863889124 0.03375055315837927; 0.0 0.0],
              [0.03375055315837913 3.128758481817603; 0.0 0.0]; dims = 3)
B_2_rbc = reshape([0.0; -0.01], 2, 1) # make sure B is a matrix
C_0_rbc = [7.824904812740593e-5, 0.0]
C_1_rbc = [0.09579643002426148 0.6746869652592109; 1.0 0.0]
C_2_rbc = cat([-0.00018554166974717046 0.0025652363153049716; 0.0 0.0],
              [0.002565236315304951 0.3132705036896446; 0.0 0.0]; dims = 3)
D_2_rbc = MvNormal(Diagonal(abs2.([0.1, 0.1])))
u0_2_rbc = zeros(2)

observables_2_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_observables.csv"),
                            ',')' |> collect

# Data and Noise

@testset "basic inference, simulated noise" begin
    prob = QuadraticStateSpaceProblem(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, u0_2_rbc,
                                      (0, size(observables_2_rbc, 2)); C_0 = C_0_rbc, C_1 = C_1_rbc,
                                      C_2 = C_2_rbc, observables_noise = D_2_rbc,
                                      observables = observables_2_rbc)
    @inferred QuadraticStateSpaceProblem(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, u0_2_rbc,
                                         (0, size(observables_2_rbc, 2)); C_0 = C_0_rbc,
                                         C_1 = C_1_rbc, C_2 = C_2_rbc, observables_noise = D_2_rbc,
                                         observables = observables_2_rbc)

    sol = solve(prob)
    @inferred solve(prob)
end

@testset "basic inference, simulated noise, no observations" begin
    T = 20
    prob = QuadraticStateSpaceProblem(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, u0_2_rbc, (0, T);
                                      C_0 = C_0_rbc, C_1 = C_1_rbc, C_2 = C_2_rbc)

    sol = solve(prob)
    @inferred solve(prob)

    # todo: add in regression tests
end


text:
using Smg2s, Test, SparseArrays

@testset "Nilpotent Matrix" begin
    @testset "Nilp 1" begin
        @test NilpMat(Nilp(2,8)) * NilpMat(Nilp(2,8)) * NilpMat(Nilp(2,8)) == spzeros(8, 8)
    end
    @testset "Nilp 2" begin
        vec=[1; 1; 0; 1; 1; 1; 0]
        nilp = Nilp(vec, 8)
        @test NilpMat(nilp) *NilpMat(nilp) * NilpMat(nilp) * NilpMat(nilp) == spzeros(8, 8)
    end
    @testset "Nilp 3" begin
        vec=[1; 1; 0; 1; 1; 1; 0]
        nilp = Nilp(vec, 3, 8)
        @test NilpMat(nilp)*NilpMat(nilp) * NilpMat(nilp) == spzeros(8, 8)
    end
    @testset "Nilp 4" begin
        nilp = Nilp(4, 3, 15)
        @test NilpMat(nilp) * NilpMat(nilp) * NilpMat(nilp) * NilpMat(nilp) * NilpMat(nilp) == spzeros(15, 15)
    end
end


text:
using BenchmarkTools
#using ReconTPCF
using StaticArrays

#floor(Int, 1/200*100)
#pwd(np.ndarray[np.double_t, ndim=2] r, np.ndarray[np.double_t, ndim=1] histo, double rmax, int nbins):
function dist1(r, histo, rmax, nbins)
#Calculate histogram of distances between particles belonging to the same structure (=set of coordinates).
    @inbounds @fastmath for i = 1:(size(r)[1]-1)
        for j = i+1:size(r)[1]
            d=0.
            @simd for k in 1:2
                tmp=r[i,k]-r[j,k]
                d+=tmp*tmp
            end
            d=sqrt(d)
            #d=sqrt(((r[i] - r[j])**2).sum())
            rIndex=ceil(Int64, d/rmax*nbins)
            if rIndex>=nbins
                #print("\n rIndex>=nbins\n")
                continue
            end
            histo[rIndex]+=1
        end
    end
    return histo
end

function dist2(r1, r2, histo, rmax::Int64, nbins::Int64)
#Calculate histogram of distances between particles belonging to different structures(=sets of coordinates).

    @inbounds @fastmath for i in 1:size(r1)[1]
        for j in 1:size(r2)[1]
            d::Float64 = hypot(r1[i,1]-r2[j,1], r1[i,2]-r2[j,2])

            #for k in 1:2
            #    tmp=r1[i,k]-r2[j,k]
            #    d+=tmp*tmp
            #end
            #d=sqrt(d)
            rIndex=ceil(Int64, d/(rmax*nbins))
            if rIndex>=nbins
                #print "\n rIndex>=nbins\n"
                continue
            end
            histo[rIndex]+=1
        end
    end
    return histo
end

function dist3(r1, r2, histo, rmax, nbins)
    @inbounds @fastmath for i in 1:size(r1)[1]
        for j in 1:size(r2)[1]
            d=0.

            tmp1=r1[i,1]-r2[j,1]
            tmp2=r1[i,2]-r2[j,2]
            d=sqrt(tmp1*tmp1 + tmp2*tmp2)
            rIndex=ceil(Int, d/rmax*nbins)
            if rIndex>=nbins
                #print "\n rIndex>=nbins\n"
                continue
            end
            histo[rIndex]+=1
        end
    end
    return histo
end

function dist4(r1, r2, histo, rmax, nbins)
    @inbounds @fastmath for i in 1:size(r1)[1]
        for j in 1:size(r2)[1]
            d=0.

            t::SVector{2,Int} = r1 - r2
            #tmp1=r1[i,1]-r2[j,1]
            #tmp2=r1[i,2]-r2[j,2]
            d=sqrt((t[1]*t[1]) + (t[2]*t[2]))
            rIndex=ceil(Int, d/rmax*nbins)
            if rIndex>=nbins
                #print "\n rIndex>=nbins\n"
                continue
            end
            histo[rIndex]+=1
        end
    end
    return histo
end

function dist5(r1, r2, histo, rmax, nbins, r1len::UInt, r2len::UInt)
    @inbounds @fastmath for i in 1:r1len
        for j in 1:r2len
            d=0.
            for k in 1:2
                tmp=r1[i,k]-r2[j,k]
                d+=tmp*tmp
            end
            d=sqrt(d)
            rIndex=ceil(Int, d/rmax*nbins)
            #if rIndex>=nbins
                #print "\n rIndex>=nbins\n"
            #    continue
            #end
            histo[rIndex]+=1
        end
    end
    return histo
end

function tile_dist(points::Array{Float64,2}, histo::Array{Int64,1}, rmax::Int64, nbins::Int64)
    bsize::Int64 = 100
    N = size(points)[1]
    histo .= 0
    if N <= bsize
        return dist1(points, histo, rmax, nbins)
    else
        histos = [zeros(Int64, length(histo)) for i in 1:Threads.nthreads()]
        #dump(histos[1])
        Threads.@threads for i in 0:(div(N, bsize)-1)
                            dist1(view(points, (i*bsize)+1:((i+1)*bsize), :),
                                 view(histos, Threads.threadid())[1],
                                 rmax,
                                 nbins)
            print("\nirange = ", (i*bsize)+1, ":", ((i+1)*bsize))
            #make thread with self_routine(myblocks)
        end
        @sync for il in 1:bsize:N+1
            for jl in il+bsize:bsize:N-bsize+1
                println("il: ", il, "\t\tir: ", il+bsize-1)
                println("jl: ", jl, "\t\tjr: ", jl+bsize-1)
                Threads.@spawn dist2(view(points, il:il+bsize-1, :),
                        view(points, jl:jl+bsize-1, :),
                        view(histos, Threads.threadid())[1],
                        rmax,
                        nbins)
            end
        end
        for i in 1:Threads.nthreads() histo .+= view(histos, i)[1] end
    end
    return histo
end

using FLoops
points = rand(1100, 2)
points2 = points.+1
hist = zeros(Int, 100)

tile_dist(points, hist, 5, 100)
@trace dist1(points, zeros(Int, 100), 100, 100)
@code_warntype dist2(points, points.+1, zeros(Int, 100), 100, 100)

sum(hist)
dist1(points, hist, 5, 100)
print(f)
view(points, 10:20, :)
points[10:20, :]

function staircase_looper(a, bsize)
    N = length(a)
    for (il, ir) in zip(1:bsize:N-1, bsize:bsize:N-1)
        for (jl, jr) in zip(ir+1:bsize:N-1, ir+bsize:bsize:N-1)
            println("il: ", il, "\t\tir: ", ir)
            println("jl: ", jl, "\t\tjr: ", jr)
        end
    end
end

staircase_looper(zeros(500), 100)

@btime dist3(points, points2, hist, 100, 100)
for i in 1:100
    dist1(points,  hist, 100, 100)
end
@btime dist5(points, points2, hist, 100, 100)

function solvequadratic(a, b, c)
    d = sqrt(b^2 - 4a*c)
    (-b - d) / 2a, (-b + d) / 2a
end

solvequadratic(-0.5, 999.5, -90000.0)

function bench_dist(N)
    points = rand(1:1.0:100, N, 2)
    points2 = points.+1
    hist = zeros(Int, 100)
    return @benchmark dist2($points, $points2, $hist, 100, 100)
end
a = []
range = 1000:1000:10000
for i in range
    #print("for i ops = ", i*i/2 , " ")
    push!(a, bench_dist(i))
end
for (bench, i) in zip(a, range)
    println("for ", abs2(i)/2, " operations, ", minimum(bench).time/(abs2(i)/2))
end

#q = findall(ones(Bool, 10,10))
#s1 = reinterpret(SVector{2,Int}, q)
#t = zeros(Bool, 20,10)
#t[11:end,:] .= true
#z = reinterpret(SVector{2,Int}, findall(t))
#s2 = reinterpret(SVector{2,Int}, q)
#dist4(s1, s2, zeros(Int, 100), 100, 100)
#@inline function euclid(a, b) (a.-b)

@btime ReconTPCF.blas_stat5(q, 100)


text:
"""
Methods to profile the SPRIGHT algorithm and this implementation of it.
"""

using SparseTransforms
using ProgressMeter
using Plots
using Profile

nrange = 4:20
num_runs = 5
methods_list = [:simple, :nso, :nso, :none]
σ = 1e-2
deltas = 0.05:0.05:0.3 |> collect

function method_name(methods_list::Array{Symbol,1})
    if methods_list[3] == :nso
        return "non-sample-optimal"
    elseif methods_list[3] == :mle
        return "maximum likelihood"
    end
end

function bvals(n::Int64, btype::Symbol)
    if btype == :constant
        return 2
    elseif btype == :scaling
        if n < 8
            return 2
        else
            return n ÷ 4
        end
    end
end

function profile_fixed_b(methods_list::Array{Symbol,1}, btype::Symbol)
    brange = map(x -> bvals(x, btype), nrange)
    times = Float64[]
    samples = Float64[]
    @showprogress for (n, b) in zip(nrange, brange)
        t = 0
        s = 0
        for i = 1:num_runs
            # locs = sample(0:2^n-1, 2^b, replace=false)
            # strengths = Float64.(rand(Uniform(0.1, 10), 2^b)) .* (-1) .^ rand(Bool, 2^b)
            # signal = LazySignal(n, locs, strengths, σ)
            signal = get_random_sparse_signal(n, 2^b, σ, 0.1, 10)
            res = @timed spright(signal, methods_list; report=true)
            s += res.value[2]
            t += res.time
        end
        push!(samples, s / (num_runs * 2 ^ n))
        push!(times, t / num_runs)
    end
    return samples, times
end

function check_transform_correctness(target, estimated, σ)
    if length(target.locs) != length(estimated)
        return false
    end
    for (l, s) in zip(target.locs, target.strengths)
        if !haskey(estimated, l) || !isapprox(estimated[l], s, atol=2*σ)
            return false
        end
    end
    return true
end

function profile_delta_sparse(methods_list::Array{Symbol,1}, deltas::Array{Float64,1})
    times_arr = Any[]
    samples_arr = Any[]
    correct_arr = Any[]
    @showprogress for δ in deltas
        times = Float64[]
        samples = Float64[]
        correct = Float64[]
        for n in nrange
            s, t, c = 0, 0, 0
            for i = 1:num_runs
                signal = get_random_delta_sparse_signal(n, σ, δ, 1, 0.5, 10.0)
                res = @timed spright(signal, methods_list; report=true)
                s += res.value[2]
                t += res.time
                c += check_transform_correctness(signal, res.value[1], σ)
            end
            push!(samples, s / (num_runs * 2 ^ n))
            push!(times, t / num_runs)
            push!(correct, c / num_runs)
        end
        push!(times_arr, times)
        push!(samples_arr, samples)
        push!(correct_arr, correct)
    end
    return samples_arr, times_arr, correct_arr
end

function plot_samples(methods_list::Array{Symbol,1}, btype::Symbol)
    brange = map(x -> bvals(x, btype), nrange)
    samples, times = profile(methods_list, btype)
    plot()
    plot!(nrange, (2 .^ brange) .* nrange ./ (2 .^ nrange); label="asymptotic trend")
    plot!(nrange, samples; label="experimental values", xlabel="log2(signal length)", ylabel="sample fraction", title="Sample ratio with $(method_name(methods_list)) decoding and $(string(btype)) sparsity")
end

function plot_delta_profile()
    samples, times = profile_delta_sparse(methods_list, deltas)
    plot()
    plot!(nrange, samples; label=hcat(deltas...), xlabel="log2(signal length)", ylabel="sample fraction", title="Sample ratio with sparsity O(N^delta)")
    brange = map(x -> 2^(-0.7x) * x^2 / 12.5, nrange)
    plot!(nrange, brange; label="asymptotic trend")
end

function plot_delta_times()
    samples, times = profile_delta_sparse(methods_list, deltas)
    plot()
    plot!(nrange, times; label=hcat(deltas...), xlabel="log2(signal length)", ylabel="runtime", title="Runtime with sparsity O(N^delta)")
    brange = map(x -> 2^(0.3x) * x^3, nrange)
    # println(times[1,1])
    # println(times[1][1])
    plot!(nrange, brange * times[1][1] / brange[1] / 5; label="asymptotic trend")
end

function plot_delta_times_m_diff()
    samples, times = profile_delta_sparse(def_methods, [0.3])
    plot()
    plot!(nrange, samples; label="identity subsampling", xlabel="log2(signal length)", ylabel="Sample ratio", title="Sample ratio with sparsity O(N^0.3)")
    def_methods[1] = :random
    samples2, times2 = profile_delta_sparse(def_methods, [0.3])
    plot!(nrange, samples2; label="random subsampling")
    brange = map(x -> 2^(-0.7x) * x^2, nrange)
    plot!(nrange, brange / brange[1]; label="asymptotic trend")
end

function N_logb_N(b)
    return n -> 2^n * n^b
end

function K_logb_N(delta, b)
    return n -> 2^(delta * n) * n^b
end

function plot_method_comparison(plot_val)
    trial_ms = ([:simple, :identity_like, :noiseless, :none],
                [:simple, :identity_like, :mle, :none],
                [:simple, :nso, :nso, :none])
    delta = 0.3
    plot()
    samples_list = []
    times_list = []

    m = trial_ms[1]
    samples, times, corrects = profile_delta_sparse(m, [delta])
    append!(samples_list, samples)
    append!(times_list, times)
    plot_dict = Dict(:samples=>(samples, "Sample ratio"),
                     :times=>(times, "Runtime"),
                     :corrects=>(corrects, "Success rate"))
    y, ylab = plot_dict[plot_val]
    println(y)
    plot!(nrange, y; label=string(m[3]), xlabel="log2(signal length)", ylabel=ylab, title=ylab*" with sparsity O(N^$delta)")

    m = trial_ms[2]
    samples, times, corrects = profile_delta_sparse(m, [delta])
    append!(samples_list, samples)
    append!(times_list, times)
    plot_dict = Dict(:samples=>(samples, "Sample ratio"),
                     :times=>(times, "Runtime"),
                     :corrects=>(corrects, "Success rate"))
    y, ylab = plot_dict[plot_val]
    println(y)
    plot!(nrange, y; label=string(m[3]), xlabel="log2(signal length)", ylabel=ylab, title=ylab*" with sparsity O(N^$delta)")

    m = trial_ms[3]
    samples, times, corrects = profile_delta_sparse(m, [delta])
    append!(samples_list, samples)
    append!(times_list, times)
    plot_dict = Dict(:samples=>(samples, "Sample ratio"),
                     :times=>(times, "Runtime"),
                     :corrects=>(corrects, "Success rate"))
    y, ylab = plot_dict[plot_val]
    println(y)
    plot!(nrange, y; label=string(m[3]), xlabel="log2(signal length)", ylabel=ylab, title=ylab*" with sparsity O(N^$delta)")

    # if plot_val != :corrects
    #     param_lookup = Dict(:samples => (2, delta-1, samples_list[3][10]),
    #                         :times => (3, delta, times_list[3][10]))
    #     b, d, norm = param_lookup[plot_val]
    #     asymp = map(K_logb_N(d, b), nrange)
    #     plot!(nrange, asymp * norm / asymp[10]; label="asymptotic trend")
    # end
    #
    # gui()

    # for m in trial_ms
    #     println(m)
    #     s, t = profile_delta_sparse(m, [0.3])
    #     println(t)
    #     plot!(nrange, t; label=string(m[3]), xlabel="log2(signal length)", ylabel="runtime", title="Runtime with sparsity O(N^0.3)")
    #     # samples_trials.append(s)
    #     # times_trials.append(t)
    # end
end


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofRepFAM
import CytofRepFAM.MCMC
=#

@testset "Test logit / sigmoid" begin
  @test MCMC.logit(.6) ≈ log(.6 / .4)
  @test MCMC.sigmoid(3.) ≈ 1 / (1 + exp(-3.))
end

@testset "Test logsumexp" begin
  x = randn(10)
  @test MCMC.logsumexp(x) ≈ log(sum(exp.(x)))
  @test MCMC.logsumexp(x[1:3]) ≈ MCMC.logsumexp(x[1], x[2], x[3])
end


text:
"""
# Get the categories for an economic data series in the FRED database.

## Arguments

    - `symbol` : String specifying the name (id) of the time series.

## Examples

```jldoctests
julia> get_category("GDPC1")
julia> get_category("FEDFUNDS")
julia> get_category("T10Y2Y")
```
"""
function get_category(symbol::String)

    url = "https://api.stlouisfed.org/fred/series/categories"
    parameters = Dict("api_key" => ENV["API_FRED"],
                      "file_type" => "json",
                      "series_id" => symbol)

    response = HTTP.request("GET", url; query = parameters)
    body = JSON.parse(String(response.body))["categories"][1]["name"]

    return(body)
end


text:
# Helper structs to dispatch known decision variable construction
struct NoSpecifiedConstraint <: MOI.AbstractSet end
Base.copy(set::NoSpecifiedConstraint) = set

struct DecisionSet <: MOI.AbstractScalarSet
    stage::Int
    constraint::Union{MOI.AbstractSet, JuMP.AbstractVectorSet}
    is_recourse::Bool

    function DecisionSet(stage::Integer; constraint::Union{MOI.AbstractSet, JuMP.AbstractVectorSet} = NoSpecifiedConstraint(), is_recourse::Bool = false)
        return new(stage, constraint, is_recourse)
    end
end

struct KnownSet <: MOI.AbstractScalarSet
    stage::Int
end

struct SingleDecisionSet{T} <: MOI.AbstractScalarSet
    stage::Int
    decision::Decision{T}
    constraint::MOI.AbstractSet
    is_recourse::Bool
end
Base.copy(set::SingleDecisionSet) = reuse(set, set.decision)

set_constraint(set::SingleDecisionSet) = set.constraint

struct MultipleDecisionSet{T} <: MOI.AbstractVectorSet
    stage::Int
    decisions::Vector{Decision{T}}
    constraint::MOI.AbstractSet
    is_recourse::Bool
end
MOI.dimension(set::MultipleDecisionSet) = length(set.decisions)
Base.copy(set::MultipleDecisionSet) = reuse(set, set.decisions)

MOIU.variable_function_type(::Type{<:SingleDecisionSet}) = SingleDecision
MOIU.variable_function_type(::Type{<:MultipleDecisionSet}) = VectorOfDecisions

is_decision_type(::Type{SingleDecisionSet}) = true
is_decision_type(::Type{MultipleDecisionSet}) = false

function JuMP.in_set_string(print_mode, set::SingleDecisionSet)
    decision_str(set) = set.is_recourse ? "RecourseDecisions" : "Decisions"
    if state(set.decision) == Taken
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))(value = $(set.decision.value))")
    end
    if state(set.decision) == Known
        return string(JuMP._math_symbol(print_mode, :in), " Known(value = $(set.decision.value))")
    end
    if set.constraint isa NoSpecifiedConstraint
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))")
    else
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function JuMP.in_set_string(print_mode, set::MultipleDecisionSet)
    decision_str(set) = set.is_recourse ? "RecourseDecisions" : "Decisions"
    if all(d -> state(d) == Taken, set.decisions)
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))(values = $([d.value for d in set.decisions]))")
    end
    if all(d -> state(d) == Known, set.decisions)
        return string(JuMP._math_symbol(print_mode, :in), " Known(values = $([d.value for d in set.decisions]))")
    end
    if set.constraint isa NoSpecifiedConstraint
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))")
    else
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function reuse(set::SingleDecisionSet, decision::Decision)
    return SingleDecisionSet(set.stage, decision, copy(set.constraint), set.is_recourse)
end

function reuse(set::MultipleDecisionSet, decisions::Vector{<:Decision})
    return MultipleDecisionSet(set.stage, decisions, copy(set.constraint), set.is_recourse)
end

function decision_set(variable::JuMP.ScalarVariable, set::DecisionSet)
    return SingleDecisionSet(set.stage, Decision(variable.info, Float64), set.constraint, set.is_recourse)
end

function decision_set(variables::Vector{<:JuMP.ScalarVariable}, set::DecisionSet)
    decisions = map(variables) do variable
        Decision(variable.info, Float64)
    end
    if set.constraint isa JuMP.AbstractVectorSet
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   JuMP.moi_set(set.constraint, length(variables)),
                                   set.is_recourse)
    else
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   set.constraint,
                                   set.is_recourse)
    end
end

function decision_set(variable::JuMP.ScalarVariable, set::KnownSet)
    return SingleDecisionSet(set.stage, KnownDecision(variable.info, Float64), NoSpecifiedConstraint(), false)
end

function decision_set(variables, set::KnownSet)
    decisions = map(variables) do variable
        KnownDecision(variable.info, Float64)
    end
    return MultipleDecisionSet(set.stage,
                               decisions,
                               NoSpecifiedConstraint(),
                               false)
end


text:
@test mysignif(SpecificG(40.0, 0.743_056_411E3), TestDigits) ≈ mysignif(SpecificG_Ph(40.0, 2700.0), TestDigits)
@test mysignif(SpecificG(40.0u"MPa", 0.743_056_411E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(40.0u"MPa", 2700.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(60.0, 0.791_137_067E3), TestDigits) ≈ mysignif(SpecificG_Ph(60.0, 2700.0), TestDigits)
@test mysignif(SpecificG(60.0u"MPa", 0.791_137_067E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(60.0u"MPa", 2700.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(60.0, 0.882_756_860E3), TestDigits) ≈ mysignif(SpecificG_Ph(60.0, 3200.0), TestDigits)
@test mysignif(SpecificG(60.0u"MPa", 0.882_756_860E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(60.0u"MPa", 3200.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(40.0, 0.743_056_411E3), TestDigits) ≈ mysignif(SpecificF_Ph(40.0, 2700.0), TestDigits)
@test mysignif(SpecificF(40.0u"MPa", 0.743_056_411E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(40.0u"MPa", 2700.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(60.0, 0.791_137_067E3), TestDigits) ≈ mysignif(SpecificF_Ph(60.0, 2700.0), TestDigits)
@test mysignif(SpecificF(60.0u"MPa", 0.791_137_067E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(60.0u"MPa", 2700.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(60.0, 0.882_756_860E3), TestDigits) ≈ mysignif(SpecificF_Ph(60.0, 3200.0), TestDigits)
@test mysignif(SpecificF(60.0u"MPa", 0.882_756_860E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(60.0u"MPa", 3200.0u"kJ/kg"), TestDigits)

@test SpecificV(40.0, 0.743_056_411E3) ≈ SpecificV_Ph(40.0, 2700.0)
@test SpecificV(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpecificV_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpecificV(60.0, 0.791_137_067E3) ≈ SpecificV_Ph(60.0, 2700.0)
@test SpecificV(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpecificV_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpecificV(60.0, 0.882_756_860E3) ≈ SpecificV_Ph(60.0, 3200.0)
@test SpecificV(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpecificV_Ph(60.0u"MPa", 3200.0u"kJ/kg")

@test SpecificU(40.0, 0.743_056_411E3) ≈ SpecificU_Ph(40.0, 2700.0)
@test SpecificU(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpecificU_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpecificU(60.0, 0.791_137_067E3) ≈ SpecificU_Ph(60.0, 2700.0)
@test SpecificU(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpecificU_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpecificU(60.0, 0.882_756_860E3) ≈ SpecificU_Ph(60.0, 3200.0)
@test SpecificU(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpecificU_Ph(60.0u"MPa", 3200.0u"kJ/kg")

@test SpecificS(40.0, 0.743_056_411E3) ≈ SpecificS_Ph(40.0, 2700.0)
@test SpecificS(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpecificS_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpecificS(60.0, 0.791_137_067E3) ≈ SpecificS_Ph(60.0, 2700.0)
@test SpecificS(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpecificS_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpecificS(60.0, 0.882_756_860E3) ≈ SpecificS_Ph(60.0, 3200.0)
@test SpecificS(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpecificS_Ph(60.0u"MPa", 3200.0u"kJ/kg")

@test SpecificCP(40.0, 0.743_056_411E3) ≈ SpecificCP_Ph(40.0, 2700.0)
@test SpecificCP(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpecificCP_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpecificCP(60.0, 0.791_137_067E3) ≈ SpecificCP_Ph(60.0, 2700.0)
@test SpecificCP(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpecificCP_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpecificCP(60.0, 0.882_756_860E3) ≈ SpecificCP_Ph(60.0, 3200.0)
@test SpecificCP(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpecificCP_Ph(60.0u"MPa", 3200.0u"kJ/kg")

@test SpecificCV(40.0, 0.743_056_411E3) ≈ SpecificCV_Ph(40.0, 2700.0)
@test SpecificCV(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpecificCV_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpecificCV(60.0, 0.791_137_067E3) ≈ SpecificCV_Ph(60.0, 2700.0)
@test SpecificCV(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpecificCV_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpecificCV(60.0, 0.882_756_860E3) ≈ SpecificCV_Ph(60.0, 3200.0)
@test SpecificCV(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpecificCV_Ph(60.0u"MPa", 3200.0u"kJ/kg")

@test SpeedOfSound(40.0, 0.743_056_411E3) ≈ SpeedOfSound_Ph(40.0, 2700.0)
@test SpeedOfSound(40.0u"MPa", 0.743_056_411E3u"K") ≈ SpeedOfSound_Ph(40.0u"MPa", 2700.0u"kJ/kg")

@test SpeedOfSound(60.0, 0.791_137_067E3) ≈ SpeedOfSound_Ph(60.0, 2700.0)
@test SpeedOfSound(60.0u"MPa", 0.791_137_067E3u"K") ≈ SpeedOfSound_Ph(60.0u"MPa", 2700.0u"kJ/kg")

@test SpeedOfSound(60.0, 0.882_756_860E3) ≈ SpeedOfSound_Ph(60.0, 3200.0)
@test SpeedOfSound(60.0u"MPa", 0.882_756_860E3u"K") ≈ SpeedOfSound_Ph(60.0u"MPa", 3200.0u"kJ/kg")


text:
using JuMP

# Inspired from Joey Huchette's test in ConvexHull.jl
function jumpsimplextest(lib::PolyhedraLibrary, n)
    m = Model()
    @variable(m, x[1:n] >= 0)
    @constraint(m, sum(x[i] for i=1:n) == 1)

    poly = polyhedron(m, lib)

    V = zeros(Int, n, n)
    for k in 1:n
        V[k, k] = 1
    end
    generator_fulltest(poly, V)
end

# Inspired from Joey Huchette's test in ConvexHull.jl
function simplexorigtest(lib::PolyhedraLibrary, n)
    m = Model()
    @variable(m, x[1:n] >= 0)
    @constraint(m, sum(x[i] for i=1:n) ≤ 1)

    poly = polyhedron(m, lib)

    V = zeros(Int, n+1, n)
    for k in 1:n
        V[k, k] = 1
    end
    generator_fulltest(poly, V)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LMDB_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LMDB")
JLLWrappers.@generate_main_file("LMDB", UUID("6206cf0b-f360-5984-af49-5437264c140e"))
end  # module LMDB_jll


text: Dict{String, Any}("nest" => Any[Any[Dict{String, Any}("value" => "a", "type" => "string")], Any[Dict{String, Any}("value" => "b", "type" => "string")]])


text:
module Interpolation

export hermite

function hermite(x, y, x0)
    diff_table = Array{eltype(y), 2}(undef, length(x), length(x))
    diff_table[:,1] = y
    for j = 1:length(x)-1
        for i = 1:length(x)-j
            diff_table[i,j+1] = (diff_table[i+1,j]-diff_table[i,j])/(x[i+j]-x[i])
        end
    end

    y0 = zero(eltype(x0))
    s = 1
    for j = 1:length(x)
        y0 += s*diff_table[1,j]
        s *= x0 .- x[j]
    end
    return y0
end

end # module Interpolation

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
function allowed_units(arg::String)
    larg = lowercase(arg)
    return larg == "inch" || larg == "mm" || larg == "cm"
end

function process_units(arg::String)
    larg = lowercase(arg)
    return larg == "inch" ? inch : larg == "mm" ? mm : larg == "cm" ? cm : error("Invalid units, this should not happen.")
end

function allowed_backend(arg::String)
    larg = lowercase(arg)
    return larg == "svg" || larg == "svgjs" || larg == "png" || larg == "pdf" || larg == "ps" || larg == "pgf"
end

function process_backend(arg::String)
    larg = lowercase(arg)
    if larg == "svg"
        return SVG
    end
    if larg == "svgjs"
        return SVGJS
    end
    if larg == "png"
        return PNG
    end
    if larg == "pdf"
        return PDF
    end
    if larg == "ps"
        return PS
    end
    if larg == "pgf"
        return PGF
    end
end

function is_windowed_data(df::DataFrame)
    ns = names(df)
    return (:WindowFirst ∈ ns) && (:WindowLast ∈ ns)
end

function is_distance_data(df::DataFrame)
    ns = names(df)
    return :Value ∈ ns
end

#=
function prepare_dictionary(names::Vector{Symbol})
    d = Dict{Symbol, Vector{Float64}}()
    for name in names
        d[name] = Vector{Float64}()
    end
    return d
end

function clock_collect(names::Vector{Symbol}, dates::Matrix{Float64})
    d = prepare_dictionary(names)
    l = length(names)
    for i = 1:l, j = i+1:l
        push!(d[names[i]], dates[i, j])
        push!(d[names[j]], dates[i, j])
    end
    return d
end
=#

function filter_by_ref(df, ref)
    a = df[:FirstSeq]
    b = df[:SecondSeq]
    idx = [(a[i] == ref) || (b[i] == ref) for i in 1:length(a)]
    return df[idx, :]
end

function sequence_names(df, ref)
    return [df[i, :FirstSeq] != ref ? df[i, :FirstSeq] : df[i, :SecondSeq] for i in 1:nrow(df)]
end

function heatplot_y_order(df, col, bysim)
    level_values = levels(df[:SeqName])
    if bysim
        means = Vector{Float64}(length(level_values))
        vals = Vector{Float64}(df[col])
        m = 1
        for sname in level_values
            idx = df[:SeqName] .== sname
            subvals = vals[idx]
            means[m] = sum(subvals)
            m += 1
        end
        o = sortperm(means)
    else
        o = sortperm(level_values)
    end
    return o
end

function heatplot(df::DataFrame, col::Symbol, legend::String)
    return (plot(df, x = :FirstSeq, y = :SecondSeq, color = col, Geom.rectbin,
                Guide.xlabel("Sequence name"), Guide.ylabel("Sequence name"),
                Guide.colorkey(legend),
                Guide.title("$(legend) between sequences")), df)
end

function heatplot(df::DataFrame, col::Symbol, ref::String, legend::String, sortsim::Bool)
    if ref == "default"
        ref = df[:FirstSeq][1]
    end

    df = filter_by_ref(df, ref)
    df[:RefName] = fill!(Vector{String}(nrow(df)), ref)
    df[:SeqName] = sequence_names(df, ref)
    pool!(df, [:RefName, :SeqName])

    idx = [isnan(i) for i in df[col]]
    df[idx, col] = NA

    complete_df = df[complete_cases(df), :]

    o = heatplot_y_order(complete_df, col, sortsim)

    return (plot(complete_df, x = :WindowFirst, y = :SeqName, color = col, Geom.rectbin,
         Guide.xlabel("Window Start (bp)"), Guide.ylabel("Sequence name"),
         Guide.colorkey(legend), Coord.cartesian(xmin = 0),
         Guide.title("$(legend) between $(ref) and other sequences (sliding window)"),
         Scale.y_discrete(order = o)), df)
end

#using Plots; gr(); sticks(linspace(0.25π,1.5π,5), rand(5), proj=:polar, yerr=.1)

function visualize(args)

    df = readtable(args["inputfile"],
                   separator = ',',
                   header = true)
    pool!(df, [:FirstSeq, :SecondSeq])

    if is_distance_data(df)
        # Data frame contains distances.
        leg = "Genetic distance"
        if is_windowed_data(df)
            # Data is computed across a sliding window.
            p, d = heatplot(df, :Value, args["reference"], leg)
        else
            p, d = heatplot(df, :Value, leg)
        end
    else
        # Data frame contains dates.
        leg = "Divergence time"
        if is_windowed_data(df)
            # Data is computed across a sliding window.
            p, d = heatplot(df, :MidEstimate, args["reference"], leg, args["sortsim"])
        else
            p, d = heatplot(df, :MidEstimate, leg)
        end
    end

    bkend = args["backend"]
    backend = process_backend(bkend)

    unit = process_units(args["units"])
    w = args["width"] * unit
    h = args["height"] * unit

    if args["table"]
        writetable("$(args["outputfile"])_plottable.csv", d)
    end

    draw(backend("$(args["outputfile"]).$bkend", w, h), p)
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHoppingHamiltonianSimple <: AbstractBondHamiltonian{L,1}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BondHoppingHamiltonianSimple{L} <: AbstractBondHamiltonian{L,1}

    # coupling
    coupling :: Float64

end

# export the type
export BondHoppingHamiltonianSimple

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# get bond term of Hamiltonian
function bondterm(
            h :: BondHoppingHamiltonianSimple{L},
            b :: AbstractBond{L,NB}
        ) :: Float64 where {N,L,NB,NS}

    # every bond has the same coupling
    return h.coupling
end

# export the function
export bondterm

function getHoppingHamiltonianSimple(
            unitcell :: U
        ) :: BondHoppingHamiltonianSimple{L} where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B}}

    # just create a new hamiltonian
    return BondHoppingHamiltonianSimple{L}(1.0)
end

export getHoppingHamiltonianSimple


text:
@testset "Checking value of constants from cpxconst.h" begin
    version_split = split(CPLEX.version(), '.')
    version = parse(Int, version_split[1])
    release = parse(Int, version_split[2])
    @test CPLEX.CPX_VERSION_VERSION == version
    @test CPLEX.CPX_VERSION_RELEASE == release
    @test div(CPLEX.CPX_VERSION, 10000) == version * 100 + release
end


text:
const PROPOSALS = [
    Proposals.Uniform(),
    Proposals.RWalk(),
    Proposals.RStagger(),
    Proposals.Slice(),
    Proposals.RSlice()
]

const BOUNDS = [
    Bounds.NoBounds(2),
    Bounds.Ellipsoid(2),
    Bounds.MultiEllipsoid(2)
]

@testset "interface - $(typeof(prop))" for prop in PROPOSALS, bound in BOUNDS
    logl(X) = -sum(x->x^2, X)
    prior(u) = 2u .- 1 # Uniform -1, to 1
    us = rand(rng, 2, 10)
    point, _bound = Bounds.rand_live(rng, bound, us)
    loglstar = logl(prior(point))
    u, v, logL = prop(rng, point, loglstar, _bound, logl, prior)
    # simple bounds checks
    @test all(x -> 0 < x < 1, u)
    @test all(x -> -1 < x < 1, v)

    # check new point actually has better likelihood
    @test logl(v) == logL ≥ loglstar
end

@testset "Uniform" begin
    # printing
    @test sprint(show, Proposals.Uniform()) == "NestedSamplers.Proposals.Uniform"
end

@testset "RWalk" begin
    prop = Proposals.RWalk()
    @test prop.scale == 1
    @test prop.ratio == 0.5
    @test prop.walks == 25

    @test_throws AssertionError Proposals.RWalk(ratio=-0.2)
    @test_throws AssertionError Proposals.RWalk(ratio=1.2)
    @test_throws AssertionError Proposals.RWalk(walks=0)
    @test_throws AssertionError Proposals.RWalk(walks=2, ratio=0.2)
    @test_throws AssertionError Proposals.RWalk(scale=-4)
end

@testset "RStagger" begin
    prop = Proposals.RStagger()
    @test prop.scale == 1
    @test prop.ratio == 0.5
    @test prop.walks == 25

    @test_throws AssertionError Proposals.RStagger(ratio=-0.2)
    @test_throws AssertionError Proposals.RStagger(ratio=1.2)
    @test_throws AssertionError Proposals.RStagger(walks=0)
    @test_throws AssertionError Proposals.RStagger(walks=2, ratio=0.2)
    @test_throws AssertionError Proposals.RStagger(scale=-4)
end

@testset "unitcheck" begin
    @test Proposals.unitcheck(rand(rng, 1000))
    @test !Proposals.unitcheck(randn(rng, 1000))

    # works with tuples, too
    @test Proposals.unitcheck((0.3, 0.6, 0.8))
end

@testset "Slice" begin
    prop = Proposals.Slice()
    @test prop.slices == 5
    @test prop.scale == 1

    @test_throws AssertionError Proposals.Slice(slices=-2)
    @test_throws AssertionError Proposals.Slice(scale=-3)
end

@testset "RSlice" begin
    prop = Proposals.RSlice()
    @test prop.slices == 5
    @test prop.scale == 1

    @test_throws AssertionError Proposals.RSlice(slices=-2)
    @test_throws AssertionError Proposals.RSlice(scale=-3)
end


text:
include("dm.jl")
include("tau_aux.jl")
include("tau.jl")
function simIntDM!(sys::System, pathEnsemble::Ensemble, i::Int64)
    for k in eachindex(pathEnsemble)
        directMethod!(pathEnsemble[k], sys.model.F, sys.model.ν, sys.times[i], sys.model.nr)
    end
end

function simIntTau!(sys::System, pathEnsemble::Ensemble, i::Int64, tauVar::TauVar)
    for k in eachindex(pathEnsemble)
        optimisedTauLeap!(sys.model, pathEnsemble[k], sys.times[i], sys.routine.tauOpt, tauVar)
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
using RvSpectMLBase
using Experimental
using Test

@testset "Experimental.jl" begin
    # Write your tests here.
end


text:
@polly function kernel_syr2k(alpha, beta, C, A, B)
    n,m = size(A)
    for i = 1:n
        for j = 1:i
            C[i,j] += beta
        end
        for k = 1:m, j =1:i
            C[i,j] += A[j,k] * alpha * B[i,k] + B[j,k] * alpha * A[i,k]
        end
    end
end

let
    m = 1000
    n = 1200

    C = zeros(Float32,n,n)
    A = zeros(Float32,n,m)
    B = zeros(Float32,n,m)

    alpha = Float32(1.5)
    beta = Float32(1.2)
    for i = 1:n, j = 1:m
        A[i,j] = ((i*j+1)%n) / n
        B[i,j] = ((i*j+2)%m) / m
    end
    for i = 1:n, j = 1:n
        C[i,j] = ((i*j+3)%n) / m
    end

    SUITE["syr2k"] = @benchmarkable kernel_syr2k($alpha, $beta, C, A, B) setup = (C = copy($C); A = copy($A); B = copy($B))
end


text:
module UCIData

import CSV
using CategoricalArrays
using DataDeps
using DataFrames
using DelimitedFiles
using Dates
import LegacyStrings

const DATA_DIR = joinpath(@__DIR__, "data")

include("convert.jl")

function dataset(datasetname)
  dataset_path = @datadep_str datasetname

  df = CSV.read("$dataset_path/data.csv", DataFrame, header=true,
                typemap=Dict(Date=>String, DateTime=>String))
  for name in names(df)
    if name[1] == 'C'
      transform!(df, name => categorical, renamecols=false)
    end
  end

  df
end

function list_dataset_types()
  filter(x -> isdir(joinpath(DATA_DIR, x)), readdir(DATA_DIR))
end

function list_datasets(datasettype)
  dir = joinpath(DATA_DIR, datasettype)
  datasets = filter(x -> splitext(x)[end] == ".jl", readdir(dir))
  map(x -> splitext(x)[1], datasets)
end

function __init__()
  for datasettype in list_dataset_types()
    for dataset in list_datasets(datasettype)
      include(joinpath(DATA_DIR, datasettype, "$dataset.jl"))
    end
  end
end

end

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
cfg = Circuitscape.init_config()

state = Dict{Int, Symbol}()
for i = 1:10
    state[i] = Symbol("state" * string(i))
end

function step1()
    println()
    println("Step 1: Choose Data Type")
    dt = ["raster", "network"]
    n = request(RadioMenu(vcat(dt)))
    cfg["data_type"] = dt[n]
    step2()
end

function step2()
    println()
    println("Step 2: Choose Modelling Mode")
    is_raster = cfg["data_type"] == "raster"
    if is_raster
        mm = ["PREVIOUS STEP", "pairwise", "advanced", "one-to-all",
                "all-to-one"]
    else
        mm = ["PREVIOUS STEP", "pairwise", "advanced"]
    end
    n = request(RadioMenu(vcat(mm)))
    n == 1 && step1()
    cfg["scenario"] = mm[n]
    step3()
end

function step3()
    println()
    println("Step 3a: Enter path to habitat file")
    opt = ["PREVIOUS STEP", "Enter path manually", "Use filepicker"]
    n = request(RadioMenu(opt))
    n == 1 && step2()
    path = n == 2 ? manualfilepicker() : filepicker()
    cfg["habitat_file"] = path
    println()
    println("Is this a resistance or conductance file?")
    is_res = ["PREVIOUS STEP", "resistances", "conductance"]
    n = request(RadioMenu(is_res))
    n == 1 && step2()
    cfg["habitat_map_is_resistances"] = is_res[n]
    step4()
end

function step4()
    println()
    is_advanced = cfg["scenario"] == "advanced"
    if !is_advanced
        println()
        println("Step 4: Enter path to focal nodes:")
        opt = ["PREVIOUS STEP", "Enter path manually", "Use filepicker"]
        n = request(RadioMenu(opt))
        n == 1 && step3()
        path = n == 2 ? manualfilepicker() : filepicker()
        cfg["point_file"] = path
    else
        println()
        println("Step 4a: Enter path to source file")
        opt = ["PREVIOUS STEP", "Enter path manually", "Use filepicker"]
        n = request(RadioMenu(opt))
        n == 1 && step3()
        path = n == 2 ? manualfilepicker() : filepicker()
        cfg["source_file"] = path
        println()
        println("Step 4b: Enter path to ground file")
        n = ["PREVIOUS STEP", "Enter path manually", "Use filepicker"]
        n == 1 && step3()
        path = n == 2 ? manualfilepicker() : filepicker()
        cfg["ground_file"] = path
    end
    step5()
end

function step5()
    println()
    println("Step 5: Choose solver")
    opt = ["PREVIOUS STEP", "cg+amg", "cholmod"]
    n = request(RadioMenu(opt))
    n == 1 && step4()
    cfg["solver"] = opt[n]
    println()
    step6()
end

function step6()
    println()
    println("Step 6: Choose number of parallel processes")
    opt = collect(1:Sys.CPU_THREADS) |> x -> string.(x)
    n = request(RadioMenu(opt))
    n > 1 && (cfg["parallelize"] = "true")
    cfg["max_parallel"] = string(n)
    step7()
end

function step7()
    println()
    println("Step 7: Choose outputs")
    opt = ["PREVIOUS STEP", "Pick outputs"]
    n = request(RadioMenu(opt))
    n == 1 && step6()
    opt = ["Current maps", "Voltage maps"]
    n = request(MultiSelectMenu(opt))
    1 in n && (cfg["write_cur_maps"] = "true")
    2 in n && (cfg["write_volt_maps"] = "true")
    step8()
end

function step8()
    println()
    println("Step 8: Choose output file name")
    opt = ["PREVIOUS STEP", "Enter output file name"]
    n = request(RadioMenu(opt))
    n == 1 && step7()
    name = readline(stdin)
    step9(name)
end

function step9(name)
    println()
    println("Step 9: Choose output folder")
    opt = ["PREVIOUS STEP", "Enter path manually", "Use folderpicker"]
    n = request(RadioMenu(opt))
    n == 1 && step3()
    path = n == 2 ? manualfolderpicker() : folderpicker()
    cfg["output_file"] = normpath(joinpath(path, name))
    step10(name, path)
end

function step10(name, path)
    println()
    println("Step 10: Would you like to run Circuitscape?")
    opt = ["Yes", "Later"]
    n = request(RadioMenu(opt))
    n == 1 && Circuitscape.compute(cfg)
    if n == 2
        println()
        Circuitscape.write_config(cfg)
        println("$name.ini written to $(normpath(path))")
        println()
        println("Would you to build another problem?")
        l = RadioMenu(["Yes", "No"])
        n = request(l)
        n == 1 && step1()
        n == 2 && return
    end
end

function start()
    logo = raw"""

    `-````````````````````````````````````````````````
    -dd```````````````````````````````````````````````
    `.dh:``````````````````````````````````````.``````
    ```:ddddddddddhhdd+:-..``.-:ddddddhhhhddddhdd.````
    ``````-dddddddd.dddddddddhhhhhddddddddddddddddd```
    ````````````````ddddddddddddddddddddddddddddhd.```
    ```````````.ddddddddddddddddddddddddddddhhhd.`````
    ````````.ddddddddddddd```````````/dddddd/dddd:````
    ````.ddddddddd-----.``````````````````````````````
    ````.dddddd.````.ooo```````/o/```````:oo-`````````
    ``.````````````/oooo/o```oo/oo/o```oooo.:/o```````
    `ooo///o.````:o/`````oooo/`````-oooo.`````ooooooo-
    ``.`````/oo/o/```````oo/````````oo.````````````..`
    `````````oo/````````oo````````.o.`````````````````

    """
    blue    = "\033[34m"
    red     = "\033[31m"
    green   = "\033[32m"
    magenta = "\033[35m"
    normal  = "\033[0m\033[0m"

    welcome_message = raw"""
    Welcome to Circuitscape.
    """
    logo = replace(logo, "d" => "$(red)d$(normal)")
    logo = replace(logo, "h" => "$(red)h$(normal)")
    logo = replace(logo, "o" => "$(magenta)o$(normal)")
    print(logo)
    print(welcome_message)

    step1()
end


text:
g = evolving_graph(Int, String, is_directed = false)
add_edge!(g, 2, 4, "t1")
add_edge!(g, 2, 5, "t1")
add_edge!(g, 2, 3, "t2")
add_edge!(g, 4, 5, "t2")
add_edge!(g, 1, 2, "t3")
add_edge!(g, 1, 4, "t3")
add_edge!(g, 2, 3, "t4")

components = weak_connected_components(g)
comps = weak_connected_components(g, false)

@test length(components) == 1
@test weak_connected(g, 5, 1)
@test !(weak_connected(g, 1, 5))

g2 = evolving_graph(Int, Int)
add_edge!(g2, 1, 2, 1)
add_edge!(g2, 1, 3, 2)
add_edge!(g2, 4, 5, 2)
add_edge!(g2, 2, 3, 3)
add_edge!(g2, 5, 6, 3)

@test length(weak_connected_components(g2)) == 2


text:
using SimpleNeRF: sinusodial_emb, NeRFConfig, NeRFModel, total
using Optimisers
using LinearAlgebra
using Test

bs = 16
freq = 4
coords = randn(Float32, 3, bs)
comb = sinusodial_emb(coords, freq)
@test size(comb) == (2*3*freq, bs)

cfg = NeRFConfig()
model = NeRFModel(;config=cfg)
x = randn(Float32, 3, bs)
d = randn(Float32, 3, bs)
density, rgb = model(x, d)
@test size(density) == (1, bs)
@test size(rgb) == (3, bs)

params = Optimisers.setup(Optimisers.ADAM(), model)
@test sqrt(total(x->sum(x.^2), model)) > 1.0f0


text:
using KernelRidgeRegression
using Gadfly
@everywhere using ParallelDataTransfer

N = 5000
x = rand(1, N) * 4π - 2π
yy = sinc.(x) # vec(sinc.(4 .* x) .+ 0.2 .* sin.(30 .* x))
y = squeeze(yy + 0.1randn(1, N), 1)
xnew = collect(-2.5π:0.01:2.5π)'

sendto(workers(), x = x, y = y)

println("KRR")
BLAS.set_num_threads(2)
mykrr = KernelRidgeRegression.fit(KernelRidgeRegression.KRR, x, y, 1e-3/5000, MLKernels.GaussianKernel(1.0))
ynew = KernelRidgeRegression.predict(mykrr, xnew)
@time mykrr = KernelRidgeRegression.fit(KernelRidgeRegression.KRR, x, y, 1e-3/5000, MLKernels.GaussianKernel(1.0))
@time ynew = KernelRidgeRegression.predict(mykrr, xnew)

println("fastKRR")
myfastkrr = KernelRidgeRegression.fit(KernelRidgeRegression.FastKRR, x, y, 4/5000, 10, MLKernels.GaussianKernel(100.0))
yfastnew = KernelRidgeRegression.predict(myfastkrr, xnew)
@time myfastkrr = KernelRidgeRegression.fit(KernelRidgeRegression.FastKRR, x, y, 4/5000, 10, MLKernels.GaussianKernel(100.0))
@time yfastnew = KernelRidgeRegression.predict(myfastkrr, xnew)

println("ParFastKRR")
BLAS.set_num_threads(1)
myparfastkrr = KernelRidgeRegression.fitPar(
    KernelRidgeRegression.FastKRR,
    N,
    i -> x[:,i], i -> y[i],
    4/5000, 10,
    MLKernels.GaussianKernel(100.0)
)
yparfastnew = KernelRidgeRegression.predict(myfastkrr, xnew)
@time myparfastkrr = KernelRidgeRegression.fitPar(
    KernelRidgeRegression.FastKRR, N,
    i -> x[:,i], i -> y[i],
    4/5000, 10,
    MLKernels.GaussianKernel(100.0)
)
@time yparfastnew = KernelRidgeRegression.predict(myfastkrr, xnew)

draw(SVG("test.svg", 20cm, 20cm), plot(
    layer(x = xnew, y = yparfastnew, Geom.line, Theme(default_color = colorant"red")),
    layer(x = xnew, y = yfastnew,    Geom.line, Theme(default_color = colorant"purple")),
    layer(x = xnew, y = ynew,        Geom.line, Theme(default_color = colorant"green")),
    layer(x = x,    y = y,           Geom.line, Theme(default_color = colorant"blue")),
    Coord.cartesian(ymin = -1.5, ymax = 1.5),
    Guide.manual_color_key(
        "Method",
        ["Data",         "KRR",           "FastKRR",        "ParFastKRR" ],
        [colorant"blue", colorant"green", colorant"purple", colorant"red"]
    )
))


text:
# https://plot.ly/javascript/getting-started

@require Revise begin
    Revise.track(Plots, joinpath(Pkg.dir("Plots"), "src", "backends", "plotly.jl"))
end

const _plotly_attr = merge_with_base_supported([
    :annotations,
    :background_color_legend, :background_color_inside, :background_color_outside,
    :foreground_color_legend, :foreground_color_guide,
    :foreground_color_grid, :foreground_color_axis,
    :foreground_color_text, :foreground_color_border,
    :foreground_color_title,
    :label,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha, :markerstrokestyle,
    :fillrange, :fillcolor, :fillalpha,
    :bins,
    :title, :title_location,
    :titlefontfamily, :titlefontsize, :titlefonthalign, :titlefontvalign,
    :titlefontcolor,
    :legendfontfamily, :legendfontsize, :legendfontcolor,
    :tickfontfamily, :tickfontsize, :tickfontcolor,
    :guidefontfamily, :guidefontsize, :guidefontcolor,
    :window_title,
    :guide, :lims, :ticks, :scale, :flip, :rotation,
    :tickfont, :guidefont, :legendfont,
    :grid, :gridalpha, :gridlinewidth,
    :legend, :colorbar, :colorbar_title,
    :marker_z, :fill_z, :levels,
    :ribbon, :quiver,
    :orientation,
    # :overwrite_figure,
    :polar,
    :normalize, :weights,
    # :contours,
    :aspect_ratio,
    :hover,
    :inset_subplots,
    :bar_width,
    :clims,
    :framestyle,
    :tick_direction,
    :camera,
  ])

const _plotly_seriestype = [
    :path, :scatter, :bar, :pie, :heatmap,
    :contour, :surface, :wireframe, :path3d, :scatter3d, :shape, :scattergl,
]
const _plotly_style = [:auto, :solid, :dash, :dot, :dashdot]
const _plotly_marker = [
    :none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle,
    :cross, :xcross, :pentagon, :hexagon, :octagon, :vline, :hline
]
const _plotly_scale = [:identity, :log10]
is_subplot_supported(::PlotlyBackend) = true
# is_string_supported(::PlotlyBackend) = true
const _plotly_framestyles = [:box, :axes, :zerolines, :grid, :none]
const _plotly_framestyle_defaults = Dict(:semi => :box, :origin => :zerolines)
function _plotly_framestyle(style::Symbol)
    if style in _plotly_framestyles
        return style
    else
        default_style = get(_plotly_framestyle_defaults, style, :axes)
        warn("Framestyle :$style is not supported by Plotly and PlotlyJS. :$default_style was cosen instead.")
        default_style
    end
end

# --------------------------------------------------------------------------------------

function add_backend_string(::PlotlyBackend)
    """
    Pkg.build("Plots")
    """
end

const _plotly_js_path = joinpath(dirname(@__FILE__), "..", "..", "deps", "plotly-latest.min.js")
const _plotly_js_path_remote = "https://cdn.plot.ly/plotly-latest.min.js"

function _initialize_backend(::PlotlyBackend; kw...)
  @eval begin
    import JSON

    _js_code = open(readstring, _plotly_js_path, "r")

    # borrowed from https://github.com/plotly/plotly.py/blob/2594076e29584ede2d09f2aa40a8a195b3f3fc66/plotly/offline/offline.py#L64-L71 c/o @spencerlyon2
    _js_script = """
        <script type='text/javascript'>
            define('plotly', function(require, exports, module) {
                $(_js_code)
            });
            require(['plotly'], function(Plotly) {
                window.Plotly = Plotly;
            });
        </script>
    """

    # if we're in IJulia call setupnotebook to load js and css
    if isijulia()
        display("text/html", _js_script)
    end

    # if isatom()
    #     import Atom
    #     Atom.@msg evaljs(_js_code)
    # end

  end
  # TODO: other initialization
end

# ----------------------------------------------------------------

const _plotly_legend_pos = KW(
    :right => [1., 0.5],
    :left => [0., 0.5],
    :top => [0.5, 1.],
    :bottom => [0.5, 0.],
    :bottomleft => [0., 0.],
    :bottomright => [1., 0.],
    :topright => [1., 1.],
    :topleft => [0., 1.]
    )

plotly_legend_pos(pos::Symbol) = get(_plotly_legend_pos, pos, [1.,1.])
plotly_legend_pos(v::Tuple{S,T}) where {S<:Real, T<:Real} = v

function plotly_font(font::Font, color = font.color)
    KW(
        :family => font.family,
        :size   => round(Int, font.pointsize*1.4),
        :color  => rgba_string(color),
    )
end

function plotly_annotation_dict(x, y, val; xref="paper", yref="paper")
    KW(
        :text => val,
        :xref => xref,
        :x => x,
        :yref => xref,
        :y => y,
        :showarrow => false,
    )
end

function plotly_annotation_dict(x, y, ptxt::PlotText; xref="paper", yref="paper")
    merge(plotly_annotation_dict(x, y, ptxt.str; xref=xref, yref=yref), KW(
        :font => plotly_font(ptxt.font),
        :xanchor => ptxt.font.halign == :hcenter ? :center : ptxt.font.halign,
        :yanchor => ptxt.font.valign == :vcenter ? :middle : ptxt.font.valign,
        :rotation => -ptxt.font.rotation,
    ))
end

# function get_annotation_dict_for_arrow(d::KW, xyprev::Tuple, xy::Tuple, a::Arrow)
#     xdiff = xyprev[1] - xy[1]
#     ydiff = xyprev[2] - xy[2]
#     dist = sqrt(xdiff^2 + ydiff^2)
#     KW(
#         :showarrow => true,
#         :x => xy[1],
#         :y => xy[2],
#         # :ax => xyprev[1] - xy[1],
#         # :ay => xy[2] - xyprev[2],
#         # :ax => 0,
#         # :ay => -40,
#         :ax => 10xdiff / dist,
#         :ay => -10ydiff / dist,
#         :arrowcolor => rgba_string(d[:linecolor]),
#         :xref => "x",
#         :yref => "y",
#         :arrowsize => 10a.headwidth,
#         # :arrowwidth => a.headlength,
#         :arrowwidth => 0.1,
#     )
# end

function plotly_scale(scale::Symbol)
    if scale == :log10
        "log"
    else
        "-"
    end
end

function shrink_by(lo, sz, ratio)
    amt = 0.5 * (1.0 - ratio) * sz
    lo + amt, sz - 2amt
end

function plotly_apply_aspect_ratio(sp::Subplot, plotarea, pcts)
    aspect_ratio = sp[:aspect_ratio]
    if aspect_ratio != :none
        if aspect_ratio == :equal
            aspect_ratio = 1.0
        end
        xmin,xmax = axis_limits(sp[:xaxis])
        ymin,ymax = axis_limits(sp[:yaxis])
        want_ratio = ((xmax-xmin) / (ymax-ymin)) / aspect_ratio
        parea_ratio = width(plotarea) / height(plotarea)
        if want_ratio > parea_ratio
            # need to shrink y
            ratio = parea_ratio / want_ratio
            pcts[2], pcts[4] = shrink_by(pcts[2], pcts[4], ratio)
        elseif want_ratio < parea_ratio
            # need to shrink x
            ratio = want_ratio / parea_ratio
            pcts[1], pcts[3] = shrink_by(pcts[1], pcts[3], ratio)
        end
        pcts
    end
    pcts
end

# this method gets the start/end in percentage of the canvas for this axis direction
function plotly_domain(sp::Subplot, letter)
    figw, figh = sp.plt[:size]
    pcts = bbox_to_pcts(sp.plotarea, figw*px, figh*px)
    pcts = plotly_apply_aspect_ratio(sp, sp.plotarea, pcts)
    i1,i2 = (letter == :x ? (1,3) : (2,4))
    [pcts[i1], pcts[i1]+pcts[i2]]
end

function plotly_axis(axis::Axis, sp::Subplot)
    letter = axis[:letter]
    framestyle = sp[:framestyle]
    ax = KW(
        :visible    => framestyle != :none,
        :title      => axis[:guide],
        :showgrid   => axis[:grid],
        :gridcolor  => rgba_string(plot_color(axis[:foreground_color_grid], axis[:gridalpha])),
        :gridwidth  => axis[:gridlinewidth],
        :zeroline   => framestyle == :zerolines,
        :zerolinecolor => rgba_string(axis[:foreground_color_axis]),
        :showline   => framestyle in (:box, :axes) && axis[:showaxis],
        :linecolor  => rgba_string(plot_color(axis[:foreground_color_axis])),
        :ticks      => axis[:tick_direction] == :out ? "outside" : "inside",
        :mirror     => framestyle == :box,
        :showticklabels => axis[:showaxis],
    )

    if letter in (:x,:y)
        ax[:domain] = plotly_domain(sp, letter)
        ax[:anchor] = "$(letter==:x ? :y : :x)$(plotly_subplot_index(sp))"
    end

    ax[:tickangle] = -axis[:rotation]
    lims = axis_limits(axis)
    ax[:range] = map(scalefunc(axis[:scale]), lims)

    if !(axis[:ticks] in (nothing, :none))
        ax[:titlefont] = plotly_font(guidefont(axis))
        ax[:type] = plotly_scale(axis[:scale])
        ax[:tickfont] = plotly_font(tickfont(axis))
        ax[:tickcolor] = framestyle in (:zerolines, :grid) || !axis[:showaxis] ? rgba_string(invisible()) : rgb_string(axis[:foreground_color_axis])
        ax[:linecolor] = rgba_string(axis[:foreground_color_axis])

        # flip
        if axis[:flip]
            ax[:autorange] = "reversed"
        end

        # ticks
        ticks = get_ticks(axis)
        if ticks != :auto
            ttype = ticksType(ticks)
            if ttype == :ticks
                ax[:tickmode] = "array"
                ax[:tickvals] = ticks
            elseif ttype == :ticks_and_labels
                ax[:tickmode] = "array"
                ax[:tickvals], ax[:ticktext] = ticks
            end
        end
    else
        ax[:showticklabels] = false
        ax[:showgrid] = false
    end

    ax
end

function plotly_polaraxis(axis::Axis)
    ax = KW(
        :visible => axis[:showaxis],
        :showline => axis[:grid],
    )

    if axis[:letter] == :x
        ax[:range] = rad2deg.(axis_limits(axis))
    else
        ax[:range] = axis_limits(axis)
        ax[:orientation] = -90
    end

    ax
end

function plotly_layout(plt::Plot)
    d_out = KW()

    w, h = plt[:size]
    d_out[:width], d_out[:height] = w, h
    d_out[:paper_bgcolor] = rgba_string(plt[:background_color_outside])
    d_out[:margin] = KW(:l=>0, :b=>20, :r=>0, :t=>20)

    d_out[:annotations] = KW[]

    for sp in plt.subplots
        spidx = plotly_subplot_index(sp)

        # add an annotation for the title... positioned horizontally relative to plotarea,
        # but vertically just below the top of the subplot bounding box
        if sp[:title] != ""
            bb = plotarea(sp)
            tpos = sp[:title_location]
            xmm = if tpos == :left
                left(bb)
            elseif tpos == :right
                right(bb)
            else
                0.5 * (left(bb) + right(bb))
            end
            titlex, titley = xy_mm_to_pcts(xmm, top(bbox(sp)), w*px, h*px)
            title_font = font(titlefont(sp), :top)
            push!(d_out[:annotations], plotly_annotation_dict(titlex, titley, text(sp[:title], title_font)))
        end

        d_out[:plot_bgcolor] = rgba_string(sp[:background_color_inside])

        # set to supported framestyle
        sp[:framestyle] = _plotly_framestyle(sp[:framestyle])

        # if any(is3d, seriesargs)
        if is3d(sp)
            azim = sp[:camera][1] - 90 #convert azimuthal to match GR behaviour
            theta = 90 - sp[:camera][2] #spherical coordinate angle from z axis
            d_out[:scene] = KW(
                Symbol("xaxis$spidx") => plotly_axis(sp[:xaxis], sp),
                Symbol("yaxis$spidx") => plotly_axis(sp[:yaxis], sp),
                Symbol("zaxis$spidx") => plotly_axis(sp[:zaxis], sp),

                #2.6 multiplier set camera eye such that whole plot can be seen
                :camera => KW(
                    :eye => KW(
                        :x => cosd(azim)*sind(theta)*2.6,
                        :y => sind(azim)*sind(theta)*2.6,
                        :z => cosd(theta)*2.6,
                    ),
                ),
            )
        elseif ispolar(sp)
            d_out[Symbol("angularaxis$spidx")] = plotly_polaraxis(sp[:xaxis])
            d_out[Symbol("radialaxis$spidx")] = plotly_polaraxis(sp[:yaxis])
        else
            d_out[Symbol("xaxis$spidx")] = plotly_axis(sp[:xaxis], sp)
            d_out[Symbol("yaxis$spidx")] = plotly_axis(sp[:yaxis], sp)
        end

        # legend
        d_out[:showlegend] = sp[:legend] != :none
        xpos,ypos = plotly_legend_pos(sp[:legend])
        if sp[:legend] != :none
            d_out[:legend] = KW(
                :bgcolor  => rgba_string(sp[:background_color_legend]),
                :bordercolor => rgba_string(sp[:foreground_color_legend]),
                :font     => plotly_font(legendfont(sp)),
                :x => xpos,
                :y => ypos
            )
        end

        # annotations
        append!(d_out[:annotations], KW[plotly_annotation_dict(ann...; xref = "x$spidx", yref = "y$spidx") for ann in sp[:annotations]])

        # series_annotations
        for series in series_list(sp)
            anns = series[:series_annotations]
            for (xi,yi,str,fnt) in EachAnn(anns, series[:x], series[:y])
                push!(d_out[:annotations], plotly_annotation_dict(
                    xi,
                    yi,
                    PlotText(str,fnt); xref = "x$spidx", yref = "y$spidx")
                )
            end
        end

        # # arrows
        # for sargs in seriesargs
        #     a = sargs[:arrow]
        #     if sargs[:seriestype] in (:path, :line) && typeof(a) <: Arrow
        #         add_arrows(sargs[:x], sargs[:y]) do xyprev, xy
        #             push!(d_out[:annotations], get_annotation_dict_for_arrow(sargs, xyprev, xy, a))
        #         end
        #     end
        # end

        if ispolar(sp)
            d_out[:direction] = "counterclockwise"
        end

        d_out
    end

    # turn off hover if nothing's using it
    if all(series -> series.d[:hover] in (false,:none), plt.series_list)
        d_out[:hovermode] = "none"
    end

    d_out
end

function plotly_layout_json(plt::Plot)
    JSON.json(plotly_layout(plt))
end

function plotly_colorscale(grad::ColorGradient, α)
    [[grad.values[i], rgb_string(grad.colors[i])] for i in 1:length(grad.colors)]
end
plotly_colorscale(c, α) = plotly_colorscale(cgrad(alpha=α), α)
# plotly_colorscale(c, alpha = nothing) = plotly_colorscale(cgrad(), alpha)

const _plotly_markers = KW(
    :rect       => "square",
    :xcross     => "x",
    :x          => "x",
    :utriangle  => "triangle-up",
    :dtriangle  => "triangle-down",
    :star5      => "star-triangle-up",
    :vline      => "line-ns",
    :hline      => "line-ew",
)

function plotly_subplot_index(sp::Subplot)
    spidx = sp[:subplot_index]
    spidx == 1 ? "" : spidx
end

# the Shape contructor will automatically close the shape. since we need it closed,
# we split by NaNs and then construct/destruct the shapes to get the closed coords
function plotly_close_shapes(x, y)
    xs, ys = nansplit(x), nansplit(y)
    for i=1:length(xs)
        shape = Shape(xs[i], ys[i])
        xs[i], ys[i] = coords(shape)
    end
    nanvcat(xs), nanvcat(ys)
end

plotly_data(v) = collect(v)
plotly_data(surf::Surface) = surf.surf
plotly_data(v::AbstractArray{R}) where {R<:Rational} = float(v)

plotly_surface_data(series::Series, a::AbstractVector) = a
plotly_surface_data(series::Series, a::AbstractMatrix) = transpose_z(series, a, false)
plotly_surface_data(series::Series, a::Surface) = plotly_surface_data(series, a.surf)

#ensures that a gradient is called if a single color is supplied where a gradient is needed (e.g. if a series recipe defines marker_z)
as_gradient(grad::ColorGradient, α) = grad
as_gradient(grad, α) = cgrad(alpha = α)

# get a dictionary representing the series params (d is the Plots-dict, d_out is the Plotly-dict)
function plotly_series(plt::Plot, series::Series)
    st = series[:seriestype]
    if st == :shape
        return plotly_series_shapes(plt, series)
    end

    sp = series[:subplot]
    d_out = KW()

    # these are the axes that the series should be mapped to
    spidx = plotly_subplot_index(sp)
    d_out[:xaxis] = "x$spidx"
    d_out[:yaxis] = "y$spidx"
    d_out[:showlegend] = should_add_to_legend(series)

    x, y = plotly_data(series[:x]), plotly_data(series[:y])
    d_out[:name] = series[:label]

    isscatter = st in (:scatter, :scatter3d, :scattergl)
    hasmarker = isscatter || series[:markershape] != :none
    hasline = st in (:path, :path3d)
    hasfillrange = st in (:path, :scatter, :scattergl) &&
        (isa(series[:fillrange], AbstractVector) || isa(series[:fillrange], Tuple))

    # for surface types, set the data
    if st in (:heatmap, :contour, :surface, :wireframe)
        for letter in [:x,:y,:z]
            d_out[letter] = plotly_surface_data(series, series[letter])
        end
    end

    d_out[:colorbar] = KW(:title => sp[:colorbar_title])

    clims = sp[:clims]
    if is_2tuple(clims)
        d_out[:zmin], d_out[:zmax] = clims
    end

    # set the "type"
    if st in (:path, :scatter, :scattergl)
        d_out[:type] = st==:scattergl ? "scattergl" : "scatter"
        d_out[:mode] = if hasmarker
            hasline ? "lines+markers" : "markers"
        else
            hasline ? "lines" : "none"
        end
        if series[:fillrange] == true || series[:fillrange] == 0 || isa(series[:fillrange], Tuple)
            d_out[:fill] = "tozeroy"
            d_out[:fillcolor] = rgba_string(series[:fillcolor])
        elseif isa(series[:fillrange], AbstractVector)
            d_out[:fill] = "tonexty"
            d_out[:fillcolor] = rgba_string(series[:fillcolor])
        elseif !(series[:fillrange] in (false, nothing))
            warn("fillrange ignored... plotly only supports filling to zero and to a vector of values. fillrange: $(series[:fillrange])")
        end
        d_out[:x], d_out[:y] = x, y

    elseif st == :bar
        d_out[:type] = "bar"
        d_out[:x], d_out[:y], d_out[:orientation] = if isvertical(series)
            x, y, "v"
        else
            y, x, "h"
        end
        d_out[:marker] = KW(:color => rgba_string(series[:fillcolor]))

    elseif st == :heatmap
        d_out[:type] = "heatmap"
        # d_out[:x], d_out[:y], d_out[:z] = series[:x], series[:y], transpose_z(series, series[:z].surf, false)
        d_out[:colorscale] = plotly_colorscale(series[:fillcolor], series[:fillalpha])
        d_out[:showscale] = hascolorbar(sp)

    elseif st == :contour
        d_out[:type] = "contour"
        # d_out[:x], d_out[:y], d_out[:z] = series[:x], series[:y], transpose_z(series, series[:z].surf, false)
        # d_out[:showscale] = series[:colorbar] != :none
        d_out[:ncontours] = series[:levels]
        d_out[:contours] = KW(:coloring => series[:fillrange] != nothing ? "fill" : "lines")
        d_out[:colorscale] = plotly_colorscale(series[:linecolor], series[:linealpha])
        d_out[:showscale] = hascolorbar(sp)

    elseif st in (:surface, :wireframe)
        d_out[:type] = "surface"
        # d_out[:x], d_out[:y], d_out[:z] = series[:x], series[:y], transpose_z(series, series[:z].surf, false)
        if st == :wireframe
            d_out[:hidesurface] = true
            wirelines = KW(
                :show => true,
                :color => rgba_string(series[:linecolor]),
                :highlightwidth => series[:linewidth],
            )
            d_out[:contours] = KW(:x => wirelines, :y => wirelines, :z => wirelines)
            d_out[:showscale] = false
        else
            d_out[:colorscale] = plotly_colorscale(series[:fillcolor], series[:fillalpha])
            d_out[:opacity] = series[:fillalpha]
            if series[:fill_z] != nothing
                d_out[:surfacecolor] = plotly_surface_data(series, series[:fill_z])
            end
            d_out[:showscale] = hascolorbar(sp)
        end

    elseif st == :pie
        d_out[:type] = "pie"
        d_out[:labels] = pie_labels(sp, series)
        d_out[:values] = y
        d_out[:hoverinfo] = "label+percent+name"

    elseif st in (:path3d, :scatter3d)
        d_out[:type] = "scatter3d"
        d_out[:mode] = if hasmarker
            hasline ? "lines+markers" : "markers"
        else
            hasline ? "lines" : "none"
        end
        d_out[:x], d_out[:y] = x, y
        d_out[:z] = plotly_data(series[:z])

    else
        warn("Plotly: seriestype $st isn't supported.")
        return KW()
    end

    # add "marker"
    if hasmarker
        d_out[:marker] = KW(
            :symbol => get(_plotly_markers, series[:markershape], string(series[:markershape])),
            # :opacity => series[:markeralpha],
            :size => 2 * series[:markersize],
            # :color => rgba_string(series[:markercolor]),
            :line => KW(
                :color => rgba_string(series[:markerstrokecolor]),
                :width => series[:markerstrokewidth],
            ),
        )

        # gotta hack this (for now?) since plotly can't handle rgba values inside the gradient
        if series[:marker_z] == nothing
            d_out[:marker][:color] = rgba_string(series[:markercolor])
        else
            # grad = ColorGradient(series[:markercolor], alpha=series[:markeralpha])
            # grad = as_gradient(series[:markercolor], series[:markeralpha])
            cmin, cmax = get_clims(sp)
            # zrange = zmax == zmin ? 1 : zmax - zmin # if all marker_z values are the same, plot all markers same color (avoids division by zero in next line)
            d_out[:marker][:color] = [clamp(zi, cmin, cmax) for zi in series[:marker_z]]
            d_out[:marker][:cmin] = cmin
            d_out[:marker][:cmax] = cmax
            d_out[:marker][:colorscale] = plotly_colorscale(series[:markercolor], series[:markeralpha])
            d_out[:marker][:showscale] = hascolorbar(sp)
        end
    end

    # add "line"
    if hasline
        d_out[:line] = KW(
            :color => rgba_string(series[:linecolor]),
            :width => series[:linewidth],
            :shape => if st == :steppre
                "vh"
            elseif st == :steppost
                "hv"
            else
                "linear"
            end,
            :dash => string(series[:linestyle]),
            # :dash => "solid",
        )
    end

    plotly_polar!(d_out, series)
    plotly_hover!(d_out, series[:hover])

    if hasfillrange
        # if hasfillrange is true, return two dictionaries (one for original
        # series, one for series being filled to) instead of one
        d_out_fillrange = deepcopy(d_out)
        d_out_fillrange[:showlegend] = false
        if isa(series[:fillrange], AbstractVector)
            d_out_fillrange[:y] = series[:fillrange]
            delete!(d_out_fillrange, :fill)
            delete!(d_out_fillrange, :fillcolor)
        else
            # if fillrange is a tuple with upper and lower limit, d_out_fillrange
            # is the series that will do the filling
            d_out_fillrange[:x], d_out_fillrange[:y] =
                concatenate_fillrange(series[:x], series[:fillrange])
            d_out_fillrange[:line][:width] = 0
            delete!(d_out, :fill)
            delete!(d_out, :fillcolor)
        end

        return [d_out_fillrange, d_out]
    else
        return [d_out]
    end
end

function plotly_series_shapes(plt::Plot, series::Series)
    d_outs = []

    # TODO: create a d_out for each polygon
    # x, y = series[:x], series[:y]

    # these are the axes that the series should be mapped to
    spidx = plotly_subplot_index(series[:subplot])
    base_d = KW()
    base_d[:xaxis] = "x$spidx"
    base_d[:yaxis] = "y$spidx"
    base_d[:name] = series[:label]
    # base_d[:legendgroup] = series[:label]

    x, y = plotly_data(series[:x]), plotly_data(series[:y])
    for (i,rng) in enumerate(iter_segments(x,y))
        length(rng) < 2 && continue

        # to draw polygons, we actually draw lines with fill
        d_out = merge(base_d, KW(
            :type => "scatter",
            :mode => "lines",
            :x => vcat(x[rng], x[rng[1]]),
            :y => vcat(y[rng], y[rng[1]]),
            :fill => "tozeroy",
            :fillcolor => rgba_string(_cycle(series[:fillcolor], i)),
        ))
        if series[:markerstrokewidth] > 0
            d_out[:line] = KW(
                :color => rgba_string(_cycle(series[:linecolor], i)),
                :width => series[:linewidth],
                :dash => string(series[:linestyle]),
            )
        end
        d_out[:showlegend] = i==1 ? should_add_to_legend(series) : false
        plotly_polar!(d_out, series)
        plotly_hover!(d_out, _cycle(series[:hover], i))
        push!(d_outs, d_out)
    end
    d_outs
end

function plotly_polar!(d_out::KW, series::Series)
    # convert polar plots x/y to theta/radius
    if ispolar(series[:subplot])
        theta, r = filter_radial_data(pop!(d_out, :x), pop!(d_out, :y), axis_limits(series[:subplot][:yaxis]))
        d_out[:t] = rad2deg.(theta)
        d_out[:r] = r
    end
end

function plotly_hover!(d_out::KW, hover)
    # hover text
    if hover in (:none, false)
        d_out[:hoverinfo] = "none"
    elseif hover != nothing
        d_out[:hoverinfo] = "text"
        d_out[:text] = hover
    end
end

# get a list of dictionaries, each representing the series params
function plotly_series_json(plt::Plot)
    slist = []
    for series in plt.series_list
        append!(slist, plotly_series(plt, series))
    end
    JSON.json(slist)
    # JSON.json(map(series -> plotly_series(plt, series), plt.series_list))
end

# ----------------------------------------------------------------

const _use_remote = Ref(false)

function html_head(plt::Plot{PlotlyBackend})
    jsfilename = _use_remote[] ? _plotly_js_path_remote : ("file://" * _plotly_js_path)
    # "<script src=\"$(joinpath(dirname(@__FILE__),"..","..","deps","plotly-latest.min.js"))\"></script>"
    "<script src=\"$jsfilename\"></script>"
end

function html_body(plt::Plot{PlotlyBackend}, style = nothing)
    if style == nothing
        w, h = plt[:size]
        style = "width:$(w)px;height:$(h)px;"
    end
    uuid = Base.Random.uuid4()
    html = """
        <div id=\"$(uuid)\" style=\"$(style)\"></div>
        <script>
        PLOT = document.getElementById('$(uuid)');
        Plotly.plot(PLOT, $(plotly_series_json(plt)), $(plotly_layout_json(plt)));
        </script>
    """
    html
end

function js_body(plt::Plot{PlotlyBackend}, uuid)
    js = """
          PLOT = document.getElementById('$(uuid)');
          Plotly.plot(PLOT, $(plotly_series_json(plt)), $(plotly_layout_json(plt)));
    """
end

# ----------------------------------------------------------------

function _show(io::IO, ::MIME"image/png", plt::Plot{PlotlyBackend})
    # show_png_from_html(io, plt)
    error("png output from the plotly backend is not supported.  Please use plotlyjs instead.")
end

function _show(io::IO, ::MIME"image/svg+xml", plt::Plot{PlotlyBackend})
    write(io, html_head(plt) * html_body(plt))
end

function _display(plt::Plot{PlotlyBackend})
    standalone_html_window(plt)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5d783968-f503-11ea-2a7f-c178016c9aee
using Pkg, DrWatson

# ╔═╡ 5da06014-f503-11ea-3247-9d56efd51315
using StatisticalRethinking

# ╔═╡ 1d31f6f8-f503-11ea-3be1-471a142ee605
md"## Clip-03-17-19t.jl"

# ╔═╡ 5d9ff58c-f503-11ea-1e45-85013f83f4b2
@quickactivate "SR2TuringPluto"

# ╔═╡ 5dad7252-f503-11ea-0fcd-2d550e087c4b
md"### snippet 3.11"

# ╔═╡ 5dae18ec-f503-11ea-3d5a-07f5924669ae
begin
	p_grid = range(0, step=0.00001, stop=1)
	prior = ones(length(p_grid))
	likelihood = [pdf(Binomial(3, p), 3) for p in p_grid]
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
	samples = sample(p_grid, Weights(posterior), length(p_grid));
end

# ╔═╡ 5db96b02-f503-11ea-38fe-41d84e68ba6e
md"### snippet 3.14-16"

# ╔═╡ 5dba1d40-f503-11ea-22f2-4d79ce1517fe
begin
	p1 = density(samples;
	  xlab="Proportion water (p)", ylab="Density", lab="density",
	  title="Mean, median & mode", leg=:topleft);
	vline!(p1, [mode(samples)], lab="mode")
	vline!(p1, [median(samples)], lab="median")
	vline!(p1, [mean(samples)], lab="mean")
end

# ╔═╡ 5dc56b70-f503-11ea-0c94-53ba51a38824
md"### snippet 3.17-19"

# ╔═╡ b2bfb106-f503-11ea-2d25-85732a809965
loss = map(p -> sum(posterior .* abs.(p .- p_grid)), p_grid);

# ╔═╡ 5dc62824-f503-11ea-0605-3ddfb6928156
begin
	m = findmin(loss)
	p2 = plot(loss;
	  xlab="Decision (p)", ylab="Expected proportional loss",
	  title="Loss value", lab="Loss function");
	scatter!([m[2]], [m[1]], lab="p_grid[$(m[2])]=$(round(m[1], digits=3))")
end;

# ╔═╡ 5dd1f8e8-f503-11ea-0bfd-759d3c65ed80
plot(p1, p2, layout=(1, 2))

# ╔═╡ 5ddca9c8-f503-11ea-2527-9d51a7cbf3a8
md"## End of clip-03-17-19t.jl"

# ╔═╡ Cell order:
# ╟─1d31f6f8-f503-11ea-3be1-471a142ee605
# ╠═5d783968-f503-11ea-2a7f-c178016c9aee
# ╠═5d9ff58c-f503-11ea-1e45-85013f83f4b2
# ╠═5da06014-f503-11ea-3247-9d56efd51315
# ╟─5dad7252-f503-11ea-0fcd-2d550e087c4b
# ╠═5dae18ec-f503-11ea-3d5a-07f5924669ae
# ╠═5db96b02-f503-11ea-38fe-41d84e68ba6e
# ╠═5dba1d40-f503-11ea-22f2-4d79ce1517fe
# ╟─5dc56b70-f503-11ea-0c94-53ba51a38824
# ╠═b2bfb106-f503-11ea-2d25-85732a809965
# ╠═5dc62824-f503-11ea-0605-3ddfb6928156
# ╠═5dd1f8e8-f503-11ea-0bfd-759d3c65ed80
# ╟─5ddca9c8-f503-11ea-2527-9d51a7cbf3a8


text:
using Test
using AcceleratedArrays
using SplitApplyCombine

@test isempty(detect_ambiguities(Base, AcceleratedArrays))

@testset "AcceleratedArrays" begin
    include("Interval.jl")
    include("MaybeVector.jl")
    include("SingleVector.jl")

    include("HashIndex.jl")
    include("UniqueHashIndex.jl")
    include("SortIndex.jl")
    include("UniqueSortIndex.jl")
end


text:
Base.fill(v; each::Integer) = collect(Iterators.flatten([fill(x, each) for x in v]))

mutable struct IndexedCombo
  keys::AbstractArray{Symbol,1}
  domain::AbstractArray
  index::AbstractArray{Integer,1}
  size::AbstractArray{Integer,1}
end

IndexedCombo() = IndexedCombo(Symbol[],[], Int64[], [0,0])
IndexedCombo()

function IndexedCombo(; kwargs...)
    if isempty(kwargs)
        return IndexedCombo(Symbol[],[], Int64[], [0,0])
    else
        keys = []; domain = []
        for (kw, val) in kwargs;
            push!(keys, kw)
            push!(domain, val)
        end
    end

    IndexedCombo(keys, domain, collect(1:*(length.(domain)...)),
      [*(length.(domain)...), length(keys)])
end

x = IndexedCombo(a=1:5,b=1:3,c=1:2)

function IndexedCombo(x::Union{AbstractArray{Pair{Symbol, Any}}, Array{Pair{Symbol,UnitRange{Int64}},1}})
    if length(x) == 0
        return IndexedCombo(Symbol[],[], Int64[], [0,0])
    else
        keys = []; domain = []
        for (kw, val) in x;
            push!(keys, kw)
            push!(domain, val)
        end
    end
    IndexedCombo(keys, domain, collect(1:*(length.(domain)...)),
      [*(length.(domain)...), length(keys)])
end
x = IndexedCombo([:a=>1:5,:b=>1:3,:c=>1:2])

function expand(x::IndexedCombo; kwargs...)
  if isempty(kwargs)
    return x
  elseif x.size[2]==0
    return IndexedCombo([kwargs...])
  else
    keys = []; domain = []
    for (kw, val) in kwargs;
        push!(keys, kw)
        push!(domain, val)
    end
  end

  foreach(y -> (y ∈ x.keys) && throw("key :$y already defined!") , keys)

  expander = *(length.(domain)...)
  outindex = (fill(x.index, each = expander) .- 1) .* expander +
    collect(Iterators.flatten(fill(1:expander,x.size[1])))
  IndexedCombo([x.keys..., keys...], [x.domain..., domain...], outindex,
    [*(length.(domain)..., x.size[1]), length(keys) + x.size[2]])

end
x = expand(IndexedCombo(), a=1:3,b=1:2)
@time x=expand(expand(IndexedCombo(), a=1:3,b=1:2), c=1:3, d=1:2, e=1:6, f=1:6)

function expand(x::IndexedCombo, y::Union{AbstractArray{Pair{Symbol, Any}}, Array{Pair{Symbol,UnitRange{Int64}},1}})

    (length(x.keys) == 0) && return IndexedCombo(y)
    (length(y) == 0) && return indexset

    keys = []; domain = []
    for (kw, val) in y;
        push!(keys, kw)
        push!(domain, val)
    end

    foreach(y -> (y ∈ x.keys) && throw("key :$y already defined!") , keys)

    expander = *(length.(domain)...)
    outindex = (fill(x.index, each = expander) .- 1) .* expander +
      collect(Iterators.flatten(fill(1:expander,x.size[1])))
    IndexedCombo([x.keys..., keys...], [x.domain..., domain...], outindex,
      [*(length.(domain)..., x.size[1]), length(keys) + x.size[2]])
end

x = expand(expand(IndexedCombo(), [:a=>1:3,:b=>1:2]), [:c=>1:4, :d=>1:5])

Base.keys(x::IndexedCombo)                = x.keys
domain(x::IndexedCombo)                   = x.domain
Base.size(x::IndexedCombo)                = x.size
Base.size(x::IndexedCombo,i::Integer)     = x.size[i]
solutions(x::IndexedCombo)                = x.index |> length

function Base.getindex(x::IndexedCombo, row::Integer, col::Integer)
    (col==0) && (row==0) && return :Valid
    (row==0) && return keys(x)[col]
    (col==0) && return row ∈ x.index
    # Divisor is calculating based on how many values remains how many times to repeat the current value
    divisor = (col+1 > size(x)[2] ? 1 : x.domain[(col+1):size(x)[2]] .|> length |> prod)
    # indexvalue finds the index to select from the column of interest
    indexvalue = (Integer(ceil(row / divisor)) - 1) % length(domain(x)[col]) + 1
    domain(x)[col][indexvalue]
end

[x[i,j] for i in 1:x.size[1], j in 0:x.size[2]]

Base.collect(x::IndexedCombo; bool::Bool=true, varnames::Bool=true) =
  [x[i,j] for i in !varnames:size(x)[1], j in !bool:size(x)[2]]

function mysplit(x, grouplength::Integer; filler::Array=[0,-1])
  (2*grouplength >= length(x)) && return x
  top = x[1 .+ 0:grouplength]
  bottom = x[reverse(length(x) .- (0:(grouplength-1)))]
  [top..., filler[[1,2,1]]..., bottom...]
end

mysplit([1,2,3,5,6,7,8], 2)
mysplit(1:15, 4)
mysplit(1:15, 5)

function Base.println(x::IndexedCombo;
    n::Integer=10,
    footer::Bool = true,
    header::Bool = true,
    everything::Bool=false)

  isodd(n) && (n=round(n÷2)*2)

  println()
  header && (join(x[0,:, bool=true], "\t") |> println)
  indexset = everything ? (1:size(x,1)) : x.index
  for i in mysplit(indexset, n÷2)
      (i > 0) && (join(x[i,:, bool=true], "\t") |> println)
      (i == 0) && (repeat( "..\t", size(x,2)) |> println)
      (i .< 0) && ("$(length(indexset) - n) rows not printed" |> println)
  end
  footer && (join(x[0,:, bool=true], "\t") |> println)
end

printall(x::IndexedCombo; I...) = println(x, everything = true)

x.index = x.index[x.index .% 9 .== 0]
println(x)
println(x, n=7)
println(x, n=15)
println(x, n=7, everything=true)

x[x.index .% 18 .== 0]

Base.getindex(x::IndexedCombo, row::Int, ::Colon; bool::Bool=false) =
  [ x[row,i] for i = !bool:size(x,2) ]

Base.getindex(x::IndexedCombo, ::Colon, col::Union{Int64,Symbol}; bool::Bool=false) =
  [ x[i,col] for i = 1:size(x)[1]]

function Base.getindex(x::IndexedCombo, row::Integer, col::Symbol)
  keymatch = findall(y -> y == col , x.keys)
  length(keymatch)==0 && throw("Symbol :$col not found")
  x[row, keymatch...]
end

Base.getindex(x::IndexedCombo, row::Integer, col::String) = x[row, Symbol(col)]
Base.getindex(x::IndexedCombo, ::Colon, col::String) = x[:, Symbol(col)]

Base.getindex(x::IndexedCombo, ::Colon, col::Union{Int64,Symbol}) =
  [ x[i,col] for i = 1:size(x)[1] ]

Base.getindex(x::IndexedCombo, ::Colon, ::Colon; bool=false, varnames=false) =
  collect(x,bool=bool,varnames=varnames)

Base.getindex(x::IndexedCombo, ::Colon) =  x.index
Base.getindex(x::IndexedCombo, y::Union{Int64,UnitRange}) = [x.index ∈ i for i ∈ y]

Base.getindex(x::IndexedCombo, y::BitArray{1}) =
  x.index[[ind ∈ (1:size(x,1))[y] for ind in x.index]]

x.index = x[(1:size(x,1)) .% 2 .== 0]
x.index = x[(1:size(x,1)) .% 3 .== 0]
x.index = x[(1:size(x,1)) .% 4 .== 0]
x.index = x[(1:size(x,1)) .% 5 .== 0]

x |> println

  [x[i,j] for i in (1:size(x)[1])[x[:] .& y], j in 1:size(x)[2]]

  [x[i,j] for i in (1:size(x)[1])[x[:] .& y], j in 1:size(x)[2]]

Base.getindex(x::IndexedCombo, ::Colon, ::Colon, ::Colon)   =
  [x[i,j] for i in (1:size(x)[1])[x[:]], j in 1:size(x)[2]]

Base.getindex(x::IndexedCombo, ::Colon, ::Colon, y::Union{Int64,Symbol,String}) =
  [x[i,y] for i in (1:size(x)[1])[x[:]]]

# Set index!
Base.setindex!(x::IndexedCombo, y::Union{Int64,UnitRange}) =  x.logical[y]

Base.setindex!(x::IndexedCombo, y::Bool, z::Integer)   = x.logical[z] = y
Base.setindex!(x::IndexedCombo, y::Bool, z::Union{UnitRange, AbstractArray}) =
  x.logical[z] .= y
Base.setindex!(x::IndexedCombo, y::Union{Array{Bool},Array{Bool,1},BitArray{1}}, ::Colon) = x.logical[:] .= y

Base.setindex!(x::IndexedCombo, y::Union{Array{Bool},Array{Bool,1}}, z::Union{UnitRange, AbstractArray}) =
  x.logical[z] = y

Base.fill(v; each::Integer) = collect(Iterators.flatten([fill(x, each) for x in v]))

function Base.range(x::IndexedCombo)
  p = Dict()
  for i in 1:size(x)[2]; p[x.keys[i]] = sort(unique(x[:,:,i])); end
  p
end

showfeasible(x::IndexedCombo) = x[:,:,:]

###################### Testing

x = IndexedCombo(a=1:3,b=1:2,c=2:4)

x[2] = false
x[3:4] = false
x[17:18] = fill(false,2)

x[:] === x.logical
x[:,:]
collect(x)
x[:,:,:]
x[:,:,:a]

x[:,:a] == x[:,"a"]

x[2,:a]
x[2,:]
x[2]
range(x)

# Testing Expand Function
x = IndexedCombo(a=1:2,b=1:2)
x[1:2] = false
collect(x)
range(x)

x2 = expand(x,c=1:2)
collect(x2)
range(x2)

# Throw an error
#x3 = expand(x, a=1:2)

# Expanding on an empty set generates a new set
expand(IndexedCombo(), a=1:2, b=1) |> collect