text:
# Mutated

@inline function mask_unambiguous(x::UInt64)
    return nibble_mask(enumerate_nibbles(x), 0x1111111111111111)
end

@inline function mask_unambiguous(a::UInt64, b::UInt64)
    return mask_unambiguous(a) & mask_unambiguous(b)
end

@inline function mutation_bitcount(x::UInt64, y::UInt64, ::DNAAlphabet{4})
    m = mask_unambiguous(x, y)
    d = (x ⊻ y) & m
    return count_nonzero_nibbles(d), count_1111_nibbles(m)
end

@inline mutation_bitcount(x::UInt64, y::UInt64, ::DNAAlphabet{2}) = BioSequences.mismatch_bitcount(x, y, DNAAlphabet{2}())

# Counting mutations
let
    @info "Compiling bit-parallel mutation counter for LongSequence{<:NucleicAcidAlphabet}"

    counter4 = quote
        c, n = mutation_bitcount(x, y, A())
        count += c
        N += n
    end

    counter2 = :(count += mutation_bitcount(x, y, A()))

    BioSequences.compile_2seq_bitpar(
        :count_mutations_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet{4}),),
        init_code = :(count = 0; N = 0),
        head_code = counter4,
        body_code = counter4,
        tail_code = counter4,
        return_code = :(return count, N)
    ) |> eval

    BioSequences.compile_2seq_bitpar(
        :count_mutations_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet{2}),),
        init_code = :(count = 0),
        head_code = counter2,
        body_code = counter2,
        tail_code = counter2,
        return_code = :(return count)
    ) |> eval
end

count_n_mutations(seqa::LongSequence{<:NucleicAcidAlphabet}, seqb::LongSequence{<:NucleicAcidAlphabet}) = count_n_mutations(promote(seqa, seqb)...)
count_n_mutations(seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet} = count_mutations_bitpar(seqa, seqb)

struct MutationCount
    n_mutations::Int
    n_sites::Int
end

MutationCount() = MutationCount(0, 0)

count_mutations(seqa::LongSequence{<:NucleicAcidAlphabet}, seqb::LongSequence{<:NucleicAcidAlphabet}, args...) = count_mutations(promote(seqa, seqb)..., args...)
count_mutations(seqa::LongSequence{DNAAlphabet{2}}, seqb::LongSequence{DNAAlphabet{2}}) = MutationCount(count_n_mutations(seqa, seqb), min(length(seqa), length(seqb)))
count_mutations(seqa::LongSequence{DNAAlphabet{4}}, seqb::LongSequence{DNAAlphabet{4}}) = MutationCount(count_n_mutations(seqa, seqb)...)

function count_mutations(A::LongSequence{T}, B::LongSequence{T}, win_size::Int) where {T<:DNAAlphabet}
    L = min(length(A), length(B))
    last = win_size
    V = Vector{MutationCount}()
    for i in last:L
        first = last - win_size + 1
        a = A[first:last]
        b = B[first:last]
        push!(V, count_mutations(a, b))
    end
    return V
end

function count_mutations(seqs::Vector{<:LongSequence}, args...)
    N = length(seqs)
    M = [count_mutations(seqs[i], seqs[j], args...) for i in 1:N for j in (i + 1):N]
    return PairwiseListMatrix(M, false, eltype(M)())
end

function count_mutations(seqs::Vector{<:LongSequence}, names::Vector{String}, args...)
    return setlabels(count_mutations(seqs, args...), names)
end

function count_mutations(seqs::Dict{String,L}, args...) where {L<:LongSequence}
    seqsvec = Vector{L}(undef, length(seqs))
    names = Vector{String}(undef, length(seqs))
    i = 1
    for (k, v) in seqs
        seqsvec[i] = v
        names[i] = k
        i += 1
    end
    return count_mutations(seqsvec, names, args...)
end


text:
@testset "NormL1 Tests" begin

    x = Float64.([1; 2; 3; 4; 5])
    xc = complex(Float64.([1; 2; 3; 4; 5]))
    weights_scalar = 2.0
    weights_vector = weights_scalar.*ones(length(x))
    params = Dict{String, Any}()
    tau = 3.14

    @testset "NormL1_primal Tests" begin

        @test NormL1_primal(x, weights_scalar, params) == 30.0
        @test NormL1_primal(x, weights_vector, params) == 30.0

        @test NormL1_primal(xc, weights_scalar, params) == 30.0
        @test NormL1_primal(xc, weights_vector, params) == 30.0

    end

    @testset "NormL1_dual Tests" begin

        @test NormL1_dual(x, weights_scalar, params) == 2.5
        @test NormL1_dual(x, weights_vector, params) == 2.5

        @test NormL1_dual(xc, weights_scalar, params) == 2.5
        @test NormL1_dual(xc, weights_vector, params) == 2.5

    end

    @testset "NormL1_project Tests" begin

        @test isapprox(NormL1_project(x, tau, weights_scalar,
                                        spgOptions())[1], [0.0; 0.0; 0.0; 0.285; 1.285])
        @test isapprox(NormL1_project(x, tau, weights_vector,
                                        spgOptions())[1], [0.0; 0.0; 0.0; 0.285; 1.285])

        @test isapprox(NormL1_project(xc, tau, weights_scalar,
                                        spgOptions())[1], complex.([0.0; 0.0; 0.0; 0.285; 1.285]))
        @test isapprox(NormL1_project(xc, tau, weights_vector,
                                        spgOptions())[1], complex.([0.0; 0.0; 0.0; 0.285; 1.285]))

    end
end

ext: jl


text:
# Progress_bar function.
# To use as a separate file: Should be added as a progress_bar.jl file in the same folder as your code. In the main code, use: include("progress_bar.jl");
# To use in Jupyter notebook, just include the funciton in a cell and run it.

# Input arguments:
# current_step = current position in the loop. Use the loop index.
# total_step = total number of steps in your loop. Note that it is assumed that the loop starts from 1 and ends at total_step.
# Nbar = how many blocks you want to include in the progress bar.

function progress_bar(current_step,total_step,Nblock)

    percentage = round(current_step*100/total_step,1);     # Calculate percentage completion

    print("\r")
    if current_step == 1
        print("Progress = 0.0\% ");
        for i = 1:Nblock
            print("\u2591");                      # Use the unicode character "light shade" to create the empty progress bar.
        end
        print("\r");
    end

    print("Progress = ",percentage, "\% ");         # Print percentage completion

    for i = 1:round(percentage*Nblock/100)        # Map the progress percentage to a 1 to Nblock grid.
        print("\u2593");                          # Use the unicode character "dark shade" to fill the progress bar.
    end

end


text:
Adapt.adapt_storage(to::PyAdaptor{T}, x::CUDA.AnyCuArray) where {T} = CUDA.CuArray(x)
Adapt.adapt_storage(to::PyAdaptor{T}, x::CUDA.StridedCuArray) where {T} = x
Adapt.adapt_storage(to::PyAdaptor{<:CUDA.CuArray}, x::FillArrays.AbstractFill) = CUDA.cu(collect(x))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CFITSIO"
version = v"4.0.0"

# Collection of sources required to build CFITSIO
sources = [
    ArchiveSource("http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c/cfitsio-$(version).tar.gz",
                  "b2a8efba0b9f86d3e1bd619f662a476ec18112b4f27cc441cc680a4e3777425e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cfitsio*
atomic_patch -p1 ../patches/configure_in.patch
atomic_patch -p1 ../patches/Makefile_in.patch
autoreconf
if [[ "${target}" == *-mingw* ]]; then
    # This is ridiculous: when CURL is enabled, CFITSIO defines a macro,
    # `TBYTE`, that has the same name as a mingw macro.  Let's rename all
    # `TBYTE` to `_TBYTE`.
    sed -i 's/\<TBYTE\>/_TBYTE/g' $(grep -lr '\<TBYTE\>')
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-reentrant
make -j${nproc} shared
make install
# Delete the static library
rm ${prefix}/lib/libcfitsio.a
# On Windows platforms, we need to move our .dll files to bin
if [[ "${target}" == *-mingw* ]]; then
    mkdir -p ${libdir}
    mv ${prefix}/lib/*.dll ${libdir}/.
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcfitsio", :libcfitsio)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LibCURL_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
using TruncatedDistributions
using Test
using Distributions
using HCubature
using PDMats

include("exampleDists.jl")

function test_distributions(dist_gen)
       for dg in dist_gen
              d, r, properties = dg()

              d_basic = BasicBoxTruncatedMvNormal(d.μ, d.Σ, r.a, r.b)
              @test length(d_basic) == properties["length"]
              @test tp(d_basic) ≈ properties["tp"]
              @test mean(d_basic) ≈ properties["mean"]

              d_recursive = RecursiveMomentsBoxTruncatedMvNormal(d.μ, d.Σ, r.a, r.b)
              @test length(d_recursive) == properties["length"]
              @test tp(d_recursive) ≈ properties["tp"]
              @test mean(d_recursive) ≈ properties["mean"]
       end
end

test_distributions(distribution_generators2)


text:
"""
NOT IMPLEMENTED YET
"""
module Floquet

end # module


text:
@everywhere using SeisDownload, Dates
using Dates, JSON

#==================================================#
# Input Parameters
MAX_MEM_PER_CPU = 2.0 # [GB] maximum allocated memory for one cpu
DownloadType = "Earthquake" # Choise of "Noise" or "Earthquake"

IsLocationBox = false

method  = "FDSN" # Method to download data.

datasource = "GEONET" # currently, only one src can be specified.

IsResponseRemove = false #whether instrumental response is removed or not
pre_filt    = (0.001, 0.002, 10.0, 20.0) #prefilter tuple used obspy remove_response: taper between f1 and f2, f3 and f4 with obspy

fodir       = "./dataset"
foname      = "NZ_GEONETquake" # data is saved at ./dataset/$foname.jld2
#==================================================#

# read the dictionary of the phase picks from Calum's json file.
dict2 = Dict()
f = open("at_least_five_stations_repicked_with_origins.json")
dict2 = JSON.parse(f)
close(f)

# start reading the catalog
Ne=length(dict2["events"])

for i=1:Ne
    # number of phase picked for that event
    Npicks = length(dict2["events"][i]["picks"])
    println("index of quake ", i, "Number of picks  ",Npicks)

    #event information
    publicid    = nothing#dict2["events"][i]["origins"][1]["publicID"]
    catalog     = "Chamberlain19"#dict2["events"][i]["origins"][1]["catalog"]
    latitude    =dict2["events"][i]["origins"][1]["latitude"]
    longitude   =dict2["events"][i]["origins"][1]["longitude"]
    depth       =dict2["events"][i]["origins"][1]["depth"]/1E3
    type=nothing
    magnitude=nothing
    isempty(dict2["events"][i]["magnitudes"])  ? magnitude : magnitude=dict2["events"][i]["magnitudes"][1]

    time       =dict2["events"][i]["origins"][1]["time"]
    origin      = Dict("publicid" => publicid, "catalog" => catalog, "time" => time,
    "latitude"=>latitude, "longitude"=>longitude, "depth"=>depth, "magnitude"=>magnitude, "type"=>type)

    # create a dictionary for each quake with channel information and time.
    ptime=Vector{String}(undef,Npicks)
    net=Vector{String}(undef,Npicks)
    sta=Vector{String}(undef,Npicks)
    cha=Vector{String}(undef,Npicks)
    src=Vector{String}(undef,Npicks)

    pickphase_dict=[]
    # loop through all phase picked
    df = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.ss")
    for j=1:Npicks
        ptime[j]=dict2["events"][i]["picks"][j]["time"][1:end-5]  # start time of pick
        net[j] = dict2["events"][i]["picks"][j]["waveform_id"]["network_code"] # network code
        sta[j] = dict2["events"][i]["picks"][j]["waveform_id"]["station_code"] # station code
        cha[j] = dict2["events"][i]["picks"][j]["waveform_id"]["channel_code"] # channel code
        src[j] = datasource # data center

        # prints for sanity checkes
        println(ptime[j]," ",net[j]," ",sta[j]," ",cha[j])
#         # make datetime object of ptime:
        starttime=DateTime(ptime[j],df)
        starttime = starttime - Dates.Second(60)
        endtime  = starttime + Dates.Second(30*60)
        println(starttime,endtime)
        pdict_temp = Dict("starttime" => starttime, "endtime" => endtime,
         "net" => net[j], "sta" => sta[j], "loc" => loc[j], "cha" => cha[j], "src" => src[j])

        push!(pickphase_dict, pdict_temp)
    end

    event_temp = Dict("origin" => origin, "pickphase" => pickphase_dict)
    push!(event, event_temp)
end

mkpath(fodir)
fopath=joinpath(fodir, foname*".jld2")

#if lat-log box or not
IsLocationBox ? reg=locationbox : reg=[]
IsResponseRemove ? pre_filt = pre_filt : pre_filt = []

InputDictionary = Dict([
      "DownloadType"    => DownloadType,
      "method"          => method,
      "event"           => event,
      "IsLocationBox"   => IsLocationBox,
      "reg"             => reg,
      "IsResponseRemove"=> IsResponseRemove,
      "pre_filt"        => pre_filt,
      "fopath"          => fopath
    ])

# mass request with input Dictionary
seisdownload(InputDictionary, MAX_MEM_PER_CPU=float(MAX_MEM_PER_CPU))


text:
function test_cuvec(backend::Backend, T)
  println("-- Testing CuVec Utilities{$T}")
  dims = (5,6,7,8)
  len = prod(dims)
  eps = 1e-5

  X = rand(T, dims)
  Y = rand(T, dims)
  X_blob = make_blob(backend, X)
  Y_blob = make_blob(backend, Y)

  println("    > mul!")
  Vec.mul!(X, Y)
  CuVec.mul!(backend, T, X_blob.ptr.p, Y_blob.ptr.p, len)
  X2 = similar(X)
  copy!(X2, X_blob)
  @test all(abs.(X-X2) .< eps)

  println("    > pow!")
  Vec.pow!(X, 2)
  CuVec.pow!(backend, T, X_blob.ptr.p, 2, len)
  copy!(X2, X_blob)
  @test all(abs.(X-X2) .< eps)

  Vec.pow!(X, convert(T, 0.75))
  CuVec.pow!(backend, T, X_blob.ptr.p, convert(T, 0.75), len)
  copy!(X2, X_blob)
  @test all(abs.(X-X2) .< eps)

  println("    > log!")
  X = max(X, convert(T,1e-20))
  copy!(X_blob, X)
  CuVec.log!(backend, X_blob)
  X2 = to_array(X_blob)
  @test all(abs.(X2 - log(X)) .< eps)
end

function test_cuvec(backend::Backend)
  test_cuvec(backend, Float32)
  test_cuvec(backend, Float64)
end

test_cuvec(backend_gpu)


text:
using Clang.Generators
using libsodium_jll

cd(@__DIR__)

include_dir = joinpath(libsodium_jll.artifact_dir, "include") |> normpath

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

headers = detect_headers(include_dir, args)

ctx = create_context(headers, args, options)

build!(ctx)


text:
module CvGeom

# package code goes here

export hom, ihom, homNorm, homNormRobust, isHomEq, norm1!, fitConic
#export Point2, InterestPoint, SRF

function hom(A)
  sz = size(A)
  o = ones(eltype(A), tuple(1, sz[2:end]...))
  return cat(1, A, o)
end

function ihom{T<:Number}(A::Vector{T})
  return A[1:end-1] / A[end]
end

function ihom{T<:Number}(A::Array{T, 2})
  return A[1:end-1,:] ./ A[end, :]
end

function homNorm{T<:Number}(a::Vector{T})
  return a / a[end]
end

function homNormRobust{T<:Number}(a::Vector{T})
  bestMag = abs(a[1])
  bestIdx = 1
  N = size(a)[1]
  for idx=2:N
    mag = abs(a[idx])
    if mag > bestMag
      bestMag = mag
      bestIdx = idx
    end
  end
  return a / a[bestIdx]
end

function norm1!{T<:FloatingPoint}(a::Vector{T})
  a = a / norm(a)
end

function isHomEq(a::Vector, b::Vector)
  if size(a) != size(b)
    return false
  end

  return homNormRobust(a) == homNormRobust(b)

end

function isPoint2OnLine{T<:Number}(p::Vector{T})
  @assert size(p) == 2
  return false
end

function fitConic{T<:Number}(points::Array{T,2})
  d,n = size(points)
  @assert d==2
  @assert n>=5

  A = Array(Float64, (n, 6))
  for i=1:n
    x,y = points[:,i]
    A[i,:] = [x*x x*y y*y x y 1]
  end

  _,d,V = svd(A, thin=false)

  return V[:,end], A, d
end

type Point2{T<:Real}
  x::T
  y::T

  Point2() = new()
  Point2(v) = new(v, v)
  Point2(x, y) = new(x, y)
end

# flesh out parametric constructors
Point2{T<:Real}(v::T) = Point2{T}(v)
Point2{T<:Real}(x::T, y::T) = Point2{T}(x,y)
# some additional constructors for convenience
Point2(x::Real, y::Real) = Point2(promote(x,y)...)
# conversion between Point2 types
import Base.convert
convert{T<:Real}(::Type{Point2{T}}, pt::Point2) = Point2{T}(convert(T, pt.x), convert(T, pt.y))

type InterestPoint{T<:Real}
  pt::Point2{T}
  scale::T
  orientation::T

  InterestPoint() = new(Point2{T}(), -1, -1)
  InterestPoint(pt, scale, orientation) = new(pt, scale, orientation)
end

# flesh out parametric constructors
InterestPoint{T<:Real}(pt::Point2{T}, s::T, o::T) = InterestPoint{T}(pt, s, o)
# some additional constructors for convenience
InterestPoint{T<:Real}(pt::Point2{T}) = (ip = InterestPoint{T}(); ip.pt = pt; return ip)
InterestPoint{T<:Real}(x::T, y::T, s::T, o::T) = InterestPoint(Point2(x,y), s, o)
InterestPoint(x::Real, y::Real, s::Real, o::Real) = InterestPoint(promote(x,y,s,o)...)
InterestPoint{T<:Real}(::Type{T}, x::Real, y::Real, s::Real, o::Real) = InterestPoint(convert((T,T,T,T), (x,y,s,o)))

#########################################
## Match (with and without image index)
#########################################
type Match_0{T<:Number}
    img1_idx::Unsigned
    img2_idx::Unsigned
    img1_pt_idx::Unsigned
    img2_pt_idx::Unsigned
    score::T
end
Match_0{T<:Number}(img1_idx, img2_idx, img1_pt_idx, img2_pt_idx, score::T) =
    Match_0{T}(img1_idx, img2_idx, img1_pt_idx, img2_pt_idx, score)

type Match_1{T<:Number}
    idx1::Uint
    idx2::Uint
    score::T
end
Match_1{T<:Number}(idx1::Uint, idx2::Uint, score::T) =
  Match_1{T}(idx1, idx2, score)

########################################
## countBits
########################################
export countBits
function countBits_naive{T<:Unsigned}(a::T)
    sum = 0
    while(bool(a))
        sum += a&1
        a >>>= 1
    end
    return sum
end

LUT8  = Uint8[countBits_naive(a) for a in 0x00:0xFF]
LUT16 = Uint8[countBits_naive(a) for a in 0x0000:0xFFFF]
@assert LUT8 == LUT16[1:256]

function countBits(a::Uint8)
    @inbounds return int(LUT8[a+1])
end

function countBits(a::Uint16)
    @inbounds return int(LUT16[a+1])
end

function countBits{T<:Unsigned}(a::T)
    sum = 0
    while(bool(a))
        @inbounds sum += LUT16[uint16(a) + 1]
        a >>>= 16
    end
    return sum
end

function countBits!(ret, a)
  assert(size(ret) == size(a))
  for (i,v) in enumerate(a)
    @inbounds ret[i] = countBits(v)
  end
end

function countBits!{S, T<:Unsigned}(ret::Vector{S}, a::Vector{T})
  assert(size(ret) == size(a))
  @simd for i=1:length(a)
    @inbounds ret[i] = countBits(a[i])
  end
end

function countBits{T<:Unsigned}(a::AbstractArray{T})
    ret = similar(a, Uint8)
    countBits!(ret, a)
    return ret
end

#########################################
## Hamming distance
#########################################
export distance_hamming, distance2_hamming, distance2_hamming!
distance_hamming{T<:Unsigned}(a::T, b::T) = uint8(countBits(a$b))
function distance_hamming{T<:Unsigned}(p1::AbstractVector{T}, p2::AbstractVector{T})
  if length(p1) != length(p2)
    error("vectors must be same length")
  end
  d = 0
  for i=1:length(p1)
    @inbounds d += countBits(p1[i] $ p2[i])
  end
  return uint8(d)
end
function distance2_hamming!{T<:Unsigned}(p1::Matrix{T}, p2::Matrix{T}, dist2::Matrix{Uint8})
  assert(size(p1,1) == size(p2,1))
  d, m = size(p1)
  n = size(p2,2)
  assert(size(dist2) == (m,n))

  for j=1:n
    @inbounds p2j = p2[:,j]
    for i=1:m
      s = 0
      @simd for k=1:d
#         @inbounds s += countBits(p1[d,i] $ p2[d,j])
        @inbounds s += LUT16[(p1[d,i] $ p2[d,j])+1]
      end
      @inbounds dist2[i,j] = s
#       dist2[i,j] = distance_hamming(p1[:,i], p2[:,j])
#       dist2[i,j] = sum(countBits(p1[:,i] $ p2j))
    end
  end
end

function distance2_hamming!(p1::Matrix{Uint8}, p2::Matrix{Uint8}, dist2::Matrix{Uint8})
  p116 = reinterpret(Uint16, p1, (size(p1,1)>>1, size(p1,2)))
  p216 = reinterpret(Uint16, p2, (size(p2,1)>>1, size(p2,2)))

  distance2_hamming!(p116, p216, dist2)
end

function distance2_hamming{T<:Unsigned}(p1::Matrix{T}, p2::Matrix{T})
    d1, n1 = size(p1)
    d2, n2 = size(p2)
    assert(d1 == d2)
    ret = Array(Uint8, d1, n1, n2)
    for j=1:n2
        for i=1:n1
            for k=1:d1
                ret[k,i,j] = countBits_naive(p1[k,i] $ p2[k,j])
            end
        end
    end
    return sum(ret,1)
end

######################################
## nearest neighbors
######################################
export knnSearch, knnSearch!

function knnSearch(tDescs, qDescs, k, dist_fn)
    assert(size(tDescs,1) == size(qDescs,1))
    assert(k >= 1)
    nT = size(tDescs,2)
    nQ = size(qDescs,2)
    assert(nT >= k)
    ret = Array(Match_1{Uint8}, k, nQ)

    for i=1:nQ
        # initialize first match with first train descriptor
        ret[1,i] = Match_1(1, i, dist_fn(tDescs[:,1], qDescs[:,i]))
        # initialize next k-1 matches with next k-1 train descriptors (with sorting)
        for j = 2:k
            mji = Match_1(j, i, dist_fn(tDescs[:,j], qDescs[:,i]))
            l = 1
            while l < j
                if mji.score < ret[l,i].score
                    ret[l+1:j,i] = ret[l:j-1,i]
                    break;
                end
                l += 1
            end
            ret[l,i] = mji
        end
        # go through the rest of the train descriptors
        for j=k+1:nT
            mji = Match_1(j, i, dist_fn(tDescs[:,j], qDescs[:,i]))
            l = 1
            for l = 1:k
                if mji.score < ret[l,i].score
                    ret[l+1:k,i] = ret[l:k-1,i]
                    ret[l,i] = mji
                    break;
                end
            end
        end
    end

    return ret
end

function knnSearch!{T, U<:Number}(tDescs::Matrix{T}, qDescs::Matrix{T}, k, dist_fn::Function, tIdx::Matrix{Int}, dist::Matrix{U})
  assert(size(tDescs,1) == size(qDescs,1))
  d = size(tDescs,1)
  assert(k >= 1)
  nT = size(tDescs,2)
  nQ = size(qDescs,2)
  assert(nT >= k)
  assert(size(tIdx) == size(dist) == (k, nQ))
  assert(typeof(dist_fn(tDescs[:,1], qDescs[:,1])) == eltype(dist))

  for i=1:nQ
    # initialize first match with first train descriptor
    tIdx[i] = 1
    q = sub(qDescs, 1:d, i)
    dist[i] = dist_fn(tDescs[:,1], q)
  end

end

end # module


text:
#################################
# uniform mixture from template #
#################################

import Gen: geometric
import LinearAlgebra
import NearestNeighbors
import Memoize: @memoize
import StatsFuns: logsumexp

struct UniformMixtureFromTemplate <: Distribution{Matrix{Float64}} end

"""
    Y = uniform_mixture_from_template(X, p_outlier, radius, bounds)

Distribution over variable-sized point clouds representing i.i.d. samples from
a mixture of uniform distributions over:

* balls centered at each point of the template point cloud `X` (inliers), and
* a big box with bounds `bounds` (outliers).

That is, this model is an infinite mixture (for `n`) of finite mixtures (for
each of the `n` points in the point cloud).

The distribution on `n` is geometric, with expected value proportional to the
number of points in the template point cloud `X`.

The inner mixture component corresponding to the box `bounds` has weight
`p_outlier`, and for each `i`, the inner mixture component corresponding to the
sphere centered at `X[:,i]` with radius `radius` has weight `(1 - p_outlier) /
m`, where `size(X) = (3, m)`.

This distribution is "collapsed," in the sense that the `logpdf` marginalizes
out the inlier/outlierness.

The returned matrix `Y` has dimensions `3 × n`.
"""
const uniform_mixture_from_template = UniformMixtureFromTemplate()

"""
Parameter of the geometric distribution for number of points in the random
point cloud `Y`, given the number of points in the template point cloud `X`
"""
function _UniformMixtureFromTemplate_p_numpoints_Y(numpoints_X::Integer)
    return min(1 - eps(), 3 / numpoints_X)
end

# Reference implementation of the sampler -- to clarify the semantics of the
# distribution, even if it's not used in code.
function Gen.random(::UniformMixtureFromTemplate, X::Matrix{Float64},
                    p_outlier::Float64, radius::Float64,
                    bounds::Tuple)
    @assert size(X, 1) == 3
    m = size(X, 2)
    n = Gen.random(geometric, _UniformMixtureFromTemplate_p_numpoints_Y(m))
    mins = [bounds[1], bounds[3], bounds[5]]
    maxs = [bounds[2], bounds[5], bounds[6]]

    outlier = (rand(n) .< p_outlier)
    return reduce(hcat, (if outlier[i]
                             mins + rand(3) .* (maxs - mins)
                         else
                             v = LinearAlgebra.normalize(randn(3))
                             X[:, rand(1:m)] + radius * rand() * v
                         end for i in 1:n);
                  init=zeros(3, 0))
end

@memoize function get_tree_from_cloud(Y)
    NearestNeighbors.KDTree(Y)
end

function Gen.logpdf(
        ::UniformMixtureFromTemplate, Y::Matrix{Float64}, X::Matrix{Float64},
        p_outlier::Float64, radius::Float64, bounds::Tuple)
    size(X, 1) == 3 || error("X must have size 3 × something, each column represents a point.  Got size(X) = $(size(X))")
    size(Y, 1) == 3 || error("Y must have size 3 × something, each column represents a point.  Got size(Y) = $(size(Y))")
    all((Y .≥ [bounds[1], bounds[3], bounds[5]]) .&
        (Y .≤ [bounds[2], bounds[4], bounds[6]])
       ) || @warn("Some points in Y were out of bounds")

    tree = get_tree_from_cloud(Y)

    m = size(X, 2)
    n = size(Y, 2)
    logp_numpoints = Gen.logpdf(geometric, n,
                                _UniformMixtureFromTemplate_p_numpoints_Y(m))
    r = radius

    # all_idxs is an array of arrays, where the array all_idxs[i] contains all indices j such that
    # the points X[:,i] and Y[:,j] are within a distance of r units of each other.
    all_idxs = NearestNeighbors.inrange(tree, X, r)

    # For each point in Y, count how many points in X it is within r units of.
    num_close_pts = zeros(n)
    for idxs in all_idxs
        num_close_pts[idxs] .= num_close_pts[idxs] .+ 1
    end

    # Compute volume of inlier sphere, outlier bounds box, and the difference.
    Volume_inlier = 4.0 /3 * pi * r^3
    Volume_outlier = (bounds[2] - bounds[1]) * (bounds[4] - bounds[3]) * (bounds[6] - bounds[5])

    # Consider each column `Y[:,j]` to represent a single sample.
    # Then `logp_y_and_inlier` is an array whose `j`th entry is
    #     log( p( the jth sample is an inlier && the jth sample equals Y[:,j] ) )
    logp_y_and_inlier = log.( ((1-p_outlier) / Volume_inlier / m) .* num_close_pts )
    logp_y_and_outlier = log(p_outlier / Volume_outlier)

    # Log probability of each point (i.e., each column of `Y`), with the
    # inlier/outlierness marginalized out.
    @assert ndims(logp_y_and_inlier) == 1
    logp_y = logsumexp.(logp_y_and_inlier,
                        logp_y_and_outlier)

    # Return product of likelihoods across points.
    return logp_numpoints + sum(logp_y)
end

struct UniformMixtureFromTemplateMultiCloud <: Distribution{Matrix{Float64}} end
const uniform_mixture_from_template_multi_cloud = UniformMixtureFromTemplateMultiCloud()

function Gen.logpdf(
        ::UniformMixtureFromTemplateMultiCloud, Y::Matrix{Float64}, X::Array{Matrix{Float64}},
        p_outlier::Float64, radius::Float64, bounds::Tuple)
   log_pdfs = [logpdf(uniform_mixture_from_template, Y, X[i], p_outlier, radius, bounds) for i in 1:length(X)]
   sum(log_pdfs)
end

export uniform_mixture_from_template, uniform_mixture_from_template_multi_cloud


text:
#### solver attributes
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsDuals) = true
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsAddConstraintAfterSolve) = true
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsAddVariableAfterSolve) = true
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsDeleteConstraint) = true
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsDeleteVariable) = true
MathOptInterface.getattribute(m::Union{MosekSolver,MosekModel},::MathOptInterface.SupportsConicThroughQuadratic) = false # though actually the solver does

#### objective
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.ObjectiveValue) = getprimalobj(m.task,m.solutions[attr.resultindex].whichsol)
MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.ObjectiveValue) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ObjectiveValue) = attr.resultindex > 0 && attr.resultindex <= length(m.solutions)

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ObjectiveBound) = getintinf(m.task,MSK_IINF_MIO_NUM_RELAX) > 0
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.ObjectiveBound) = getdouinf(m.task,MSK_DINF_MIO_OBJ_BOUND)

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.RelativeGap) = getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP) > -1.0
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.RelativeGap) = getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP)

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.SolveTime) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.SolveTime) = true
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.SolveTime) = getdouinf(m.task,MSK_DINF_OPTIMIZER_TIME)

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.ObjectiveSense) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ObjectiveSense)  = true
MathOptInterface.cansetattribute(m::MosekSolver,attr::MathOptInterface.ObjectiveSense) = true
MathOptInterface.cansetattribute(m::MosekModel,attr::MathOptInterface.ObjectiveSense)  = true

# NOTE: The MOSEK interface currently only supports Min and Max
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.ObjectiveSense)
    sense = getobjsense(m.task)
    if sense == MSK_OBJECTIVE_SENSE_MINIMIZE
        MathOptInterface.MinSense
    else
        MathOptInterface.MaxSense
    end
end

function MathOptInterface.setattribute!(m::MosekModel,attr::MathOptInterface.ObjectiveSense, sense::MathOptInterface.OptimizationSense)
    if sense == MathOptInterface.MinSense
        setobjsense(m.task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    elseif sense == MathOptInterface.MaxSense
        setobjsense(m.task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
    else
        error("Sense '$sense' is not supported")
    end
end

#### Solver/Solution information

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.SimplexIterations) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.SimplexIterations) = true
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.SimplexIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_SIMPLEX_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_SIM_PRIMAL_ITER) + getintinf(m.task,MSK_IINF_SIM_DUAL_ITER))
    end
end

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.BarrierIterations) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.BarrierIterations) = true
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.BarrierIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_INTPNT_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_INTPNT_ITER))
    end
end

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.NodeCount) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.NodeCount) = true
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.NodeCount)
        Int(getintinf(m.task,MSK_IINF_MIO_NUM_BRANCH))
end

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.RawSolver) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.RawSolver) = true
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.RawSolver) = m.task

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.ResultCount) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ResultCount) = true
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.ResultCount) = length(m.solutions)

#### Problem information

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.NumberOfVariables) = true
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.NumberOfVariables) = true
MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.NumberOfVariables) = m.publicnumvar

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.NumberOfConstraints) = true
MathOptInterface.cangetattribute{F,D}(m::MosekModel,attr::MathOptInterface.NumberOfConstraints{F,D}) = true
MathOptInterface.getattribute{F,D}(m::MosekModel,attr::MathOptInterface.NumberOfConstraints{F,D}) = length(select(m.constrmap,F,D))

#MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.ListOfVariableReferences) = false
#MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.ListOfConstraints) = false

MathOptInterface.cangetattribute{F,D}(m::MosekSolver,attr::MathOptInterface.ListOfConstraintReferences{F,D}) = true
MathOptInterface.getattribute{F,D}(m::MosekSolver,attr::MathOptInterface.ListOfConstraintReferences{F,D}) = keys(select(m.constrmap,F,D))

#### Warm start values

MathOptInterface.cansetattribute(m::MosekSolver,attr::MathOptInterface.VariablePrimalStart) = true
MathOptInterface.cansetattribute(m::MosekModel,attr::MathOptInterface.VariablePrimalStart) = true
function MathOptInterface.setattribute!(m::MosekModel,attr::MathOptInterface.VariablePrimalStart, v :: MathOptInterface.VariableReference, val::Float64)
    subj = getindexes(m.x_block,ref2id(v))

    xx = Float64[val]
    for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
        putxxslice(m.task,sol,Int(subj),Int(subj+1),xx)
    end
end

function MathOptInterface.setattribute!(m::MosekModel,attr::MathOptInterface.VariablePrimalStart, vs::Vector{MathOptInterface.VariableReference}, vals::Vector{Float64})
    subj = Array{Int}(length(vs))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(vs[i]),subj,i)
    end

    for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
        if solutiondef(m.task,sol)
            xx = getxx(m.task,sol)
            xx[subj] = vals
            putxx(m.task,sol,xx)
        else
            xx = zeros(Float64,getnumvar(m.task))
            xx[subj] = vals
            putxx(m.task,sol,xx)
        end
    end
end

MathOptInterface.cansetattribute(m::MosekSolver,attr::MathOptInterface.ConstraintPrimalStart) = false # not sure what exactly this would be...
MathOptInterface.cansetattribute(m::MosekModel,attr::MathOptInterface.ConstraintPrimalStart) = false

MathOptInterface.cansetattribute(m::MosekSolver,attr::MathOptInterface.ConstraintDualStart) = false # for now
MathOptInterface.cansetattribute(m::MosekModel,attr::MathOptInterface.ConstraintDualStart) = false

# function MathOptInterface.setattribute!(m::MosekModel,attr::MathOptInterface.ConstraintDualStart, vs::Vector{MathOptInterface.ConstraintReference}, vals::Vector{Float64})
#     subj = Array{Int}(length(vs))
#     for i in 1:length(subj)
#         getindexes(m.x_block,ref2id(vs[i]),subj,i)
#     end

#     for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
#         if solutiondef(m.task,sol)
#             xx = getxx(m.task,sol)
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         else
#             xx = zeros(Float64,getnumvar(m.task))
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         end
#     end
# end

#### Variable solution values

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.VariablePrimal) = attr.N > 0 && attr.N <= length(m.solutions)
MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.VariablePrimal) = true

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.VariablePrimal,vs::Vector{MathOptInterface.VariableReference}) = MathOptInterface.cangetattribute(m,attr)
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.VariablePrimal,v::MathOptInterface.VariableReference) = MathOptInterface.cangetattribute(m,attr)

function MathOptInterface.getattribute!(output::Vector{Float64},m::MosekModel,attr::MathOptInterface.VariablePrimal, vs::Vector{MathOptInterface.VariableReference})
    subj = Array{Int}(length(vs))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(vs[i]),subj,i)
    end

    output[1:length(output)] = m.solutions[attr.N].xx[subj]
end

function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.VariablePrimal, vs::Vector{MathOptInterface.VariableReference})
    output = Vector{Float64}(length(vs))
    MathOptInterface.getattribute!(output,m,attr,vs)
    output
end

function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.VariablePrimal, vref::MathOptInterface.VariableReference)
    subj = getindexes(m.x_block,ref2id(vref))[1]
    m.solutions[attr.N].xx[subj]
end

MathOptInterface.cangetattribute(m::MosekSolver,attr::MathOptInterface.VariableBasisStatus) = false # is a solution selector missing?
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.VariableBasisStatus)  = false

#### Constraint solution values
function cangetattribute(m::MosekModel,attr::MathOptInterface.ConstraintPrimal)
    num = 0
    if solutiondef(m.task,MSK_SOL_ITG) num += 1 end
    if solutiondef(m.task,MSK_SOL_BAS) num += 1 end
    if solutiondef(m.task,MSK_SOL_ITR) num += 1 end
    attr.N > 0 && attr.N <= num
end

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ConstraintPrimal, c) = cangetattribute(m,attr)
MathOptInterface.cangetattribute{F,D}(m::MosekModel,attr::MathOptInterface.ConstraintPrimal, cref::MathOptInterface.ConstraintReference{F,D}) = cangetattribute(m,attr)
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ConstraintPrimal, cref::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.LessThan{Float64}}) = cangetattribute(m,attr)

function MathOptInterface.getattribute{D}(
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintPrimal,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.SingleVariable,D})

    conid = ref2id(cref)
    idxs  = getindexes(m.xc_block,conid)
    subj  = m.xc_idxs[idxs[1]]

    m.solutions[attr.N].xx[subj]
end

# Semidefinite domain for a variable
function MathOptInterface.getattribute!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintPrimal,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,MathOptInterface.PositiveSemidefiniteConeTriangle})

    cid = ref2id(cref)
    assert(cid < 0)

    # It is in fact a real constraint and cid is the id of an ordinary constraint
    barvaridx = - m.c_slack[-cid]
    output[1:length(output)] = getbarxj(m.task,m.solutions[attr.N].whichsol,barvaridx)
end

# Any other domain for variable vector
function MathOptInterface.getattribute!{D}(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintPrimal,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,D})

    xcid = ref2id(cref)
    assert(xcid > 0)

    mask = m.xc_bounds[xcid]
    idxs = getindexes(m.xc_block,xcid)
    subj = m.xc_idxs[idxs]

    output[1:length(output)] = m.solutions[attr.N].xx[subj]
end

function MathOptInterface.getattribute{D}(m     ::MosekModel,
                                          attr  ::MathOptInterface.ConstraintPrimal,
                                          cref  ::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D})
    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)[1]
    m.solutions[attr.N].xc[subi]
end

function MathOptInterface.getattribute!{D}(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintPrimal,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D})

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    if     m.c_block_slack[cid] == 0 # no slack
        output[1:length(output)] = m.solutions[attr.N].xc[subi] + m.c_constant[subi]
    elseif m.c_block_slack[cid] >  0 # qcone slack
        xid = m.c_block_slack[cid]
        xsubj = getindexes(m.x_block, xid)
        output[1:length(output)] = m.solutions[attr.N].xx[xsubj]  + m.c_constant[subi]
    elseif m.c_block_slack[cid]  # psd slack
        xid = - m.c_block_slack[cid]
        output[1:length(output)] = getbarxj(m.task,m.solutions[attr.N].whichsol,Int32(xid)) + m.c_constant[subi]
    end
end

MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ConstraintDual) = attr.N > 0 && attr.N <= length(m.solutions)
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.ConstraintDual, crefs::Vector{MathOptInterface.ConstraintReference}) = MathOptInterface.cangetattribute(m,attr)
MathOptInterface.cangetattribute(m::MosekModel,
                                 attr::MathOptInterface.ConstraintDual,
                                 cref::MathOptInterface.ConstraintReference{F,D}) where {F,D} = MathOptInterface.cangetattribute(m,attr)

function MathOptInterface.getattribute(
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintDual,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.SingleVariable,D}) where { D <: MathOptInterface.AbstractSet }

    xcid = ref2id(cref)
    idxs = getindexes(m.xc_block,xcid) # variable ids

    assert(blocksize(m.xc_block,xcid) > 0)

    subj  = getindexes(m.x_block, m.xc_idxs[idxs][1])[1]
    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        if     m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
            m.solutions[attr.N].slx[subj] - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_lower) != 0
            m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_upper) != 0
            - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_cone) != 0
            m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    else
        if     m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
            m.solutions[attr.N].sux[subj] - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_lower) != 0
            - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_upper) != 0
            m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_cone) != 0
            - m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    end
end

# Semidefinite domain for a variable
function MathOptInterface.getattribute!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintDual,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,MathOptInterface.PositiveSemidefiniteConeTriangle})

    whichsol = getsolcode(m,attr.N)
    cid = ref2id(cref)
    assert(cid < 0)

    # It is in fact a real constraint and cid is the id of an ordinary constraint
    barvaridx = - m.c_slack[-cid]
    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        output[1:length(output)] = getbarsj(m.task,whichsol,barvaridx)
    else
        output[1:length(output)] = - getbarsj(m.task,whichsol,barvaridx)
    end
end

# Any other domain for variable vector
function MathOptInterface.getattribute!{D}(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintDual,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,D})

    xcid = ref2id(cref)
    assert(xcid > 0)

    mask = m.xc_bounds[xcid]
    idxs = getindexes(m.xc_block,xcid)
    subj = m.xc_idxs[idxs]

    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        if     m.xc_bounds[conid] & boundflag_lower != 0 && m.xc_bounds[conid] & boundflag_upper != 0
            output[1:length(output)] = m.solutions[attr.N].slx[subj] - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[conid] & boundflag_lower) != 0
            output[1:length(output)] = m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[conid] & boundflag_upper) != 0
            output[1:length(output)] = - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[conid] & boundflag_cone) != 0
            output[1:length(output)] = m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    else
        if     m.xc_bounds[conid] & boundflag_lower != 0 && m.xc_bounds[conid] & boundflag_upper != 0
            output[1:length(output)] = m.solutions[attr.N].sux[subj] - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[conid] & boundflag_lower) != 0
            output[1:length(output)] = - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[conid] & boundflag_upper) != 0
            output[1:length(output)] = m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[conid] & boundflag_cone) != 0
            output[1:length(output)] = - m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    end
end

function MathOptInterface.getattribute{D}(m     ::MosekModel,
                                          attr  ::MathOptInterface.ConstraintDual,
                                          cref  ::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D})
    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)[1]

    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        m.solutions[attr.N].y[subi]
    else
        - m.solutions[attr.N].y[subi]
    end
end

function MathOptInterface.getattribute(
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintDual,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D}) where { D <: MathOptInterface.AbstractSet }

    n = blocksize(m.c_block,ref2id(cref))
    res = Vector{Float64}(n)
    MathOptInterface.getattribute!(res,m,attr,cref)
    res
end
function MathOptInterface.getattribute!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MathOptInterface.ConstraintDual,
    cref  ::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D}) where { D <: MathOptInterface.AbstractSet }

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        if     m.c_block_slack[cid] == 0 # no slack
            output[1:length(output)] = m.solutions[attr.N].y[subi]
        elseif m.c_block_slack[cid] >  0 # qcone slack
            xid = m.c_block_slack[cid]
            subj = getindexes(m.x_block, xid)
            output[1:length(output)] = m.solutions[attr.N].snx[xsubj]
        elseif m.c_block_slack[cid]  # psd slack
            xid = - m.c_block_slack[cid]
            output[1:length(output)] = getbarsj(m.task,m.solutions[attr.N].whichsol,Int32(xid))
        end
    else
        if     m.c_block_slack[cid] == 0 # no slack
            output[1:length(output)] = - m.solutions[attr.N].y[subi]
        elseif m.c_block_slack[cid] >  0 # qcone slack
            xid = m.c_block_slack[cid]
            subj = getindexes(m.x_block, xid)
            output[1:length(output)] = - m.solutions[attr.N].snx[subj]
        elseif m.c_block_slack[cid]  # psd slack
            xid = - m.c_block_slack[cid]
            output[1:length(output)] = - getbarsj(m.task,m.solutions[attr.N].whichsol,Int32(xid))
        end
    end
end

solsize{F <: MathOptInterface.AbstractScalarFunction,D}(m::MosekModel, cref :: MathOptInterface.ConstraintReference{F,D}) = 1
function solsize{D}(m::MosekModel, cref :: MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,D})
    cid = ref2id(cref)
    if cid < 0
        blocksize(m.c_block,-cid)
    else
        blocksize(m.xc_block,cid)
    end
end

function solsize{D}(m::MosekModel, cref :: MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction,D})
    cid = ref2id(cref)
    blocksize(m.c_block,cid)
end

function MathOptInterface.getattribute{F <: MathOptInterface.AbstractScalarFunction,D}(m::MosekModel,attr::MathOptInterface.ConstraintPrimal, cref :: MathOptInterface.ConstraintReference{F,D})
    cid = ref2id(cref)
    output = Vector{Float64}(solsize(m,cref))
    MathOptInterface.getattribute!(output,m,attr,cref)
    output
end

#### Status codes
MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.TerminationStatus) = true
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.TerminationStatus)
    if     m.trm == MSK_RES_OK
        MathOptInterface.Success
    elseif m.trm == MSK_RES_TRM_MAX_ITERATIONS
        MathOptInterface.IterationLimit
    elseif m.trm == MSK_RES_TRM_MAX_TIME
        MathOptInterface.TimeLimit
    elseif m.trm == MSK_RES_TRM_OBJECTIVE_RANGE
        MathOptInterface.ObjectiveLimit
    elseif m.trm == MSK_RES_TRM_MIO_NEAR_REL_GAP
        MathOptInterface.AlmostSuccess
    elseif m.trm == MSK_RES_TRM_MIO_NEAR_ABS_GAP
        MathOptInterface.AlmostSuccess
    elseif m.trm == MSK_RES_TRM_MIO_NUM_RELAXS
        MathOptInterface.OtherLimit
    elseif m.trm == MSK_RES_TRM_MIO_NUM_BRANCHES
        MathOptInterface.NodeLimit
    elseif m.trm == MSK_RES_TRM_NUM_MAX_NUM_INT_SOLUTIONS
        MathOptInterface.SolutionLimit
    elseif m.trm == MSK_RES_TRM_STALL
        MathOptInterface.SlowProgress
    elseif m.trm == MSK_RES_TRM_USER_CALLBACK
        MathOptInterface.Interrupted
    elseif m.trm == MSK_RES_TRM_MAX_NUM_SETBACKS
        MathOptInterface.OtherLimit
    elseif m.trm == MSK_RES_TRM_NUMERICAL_PROBLEM
        MathOptInterface.SlowProgress
    elseif m.trm == MSK_RES_TRM_INTERNAL
        MathOptInterface.OtherError
    elseif m.trm == MSK_RES_TRM_INTERNAL_STOP
        MathOptInterface.OtherError
    else
        MathOptInterface.OtherError
    end
end

function MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.PrimalStatus)
    if attr.N < 0 || attr.N > length(m.solutions)
        false
    else
        solsta = m.solutions[attr.N].solsta
        if     solsta == MSK_SOL_STA_UNKNOWN true
        elseif solsta == MSK_SOL_STA_OPTIMAL true
        elseif solsta == MSK_SOL_STA_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL true
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER false
        elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER true
        elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL true
        elseif solsta == MSK_SOL_STA_NEAR_INTEGER_OPTIMAL true
        else false
        end
    end
end
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.PrimalStatus)
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL
        MathOptInterface.NearlyFeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_FEAS
        MathOptInterface.NearlyFeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_FEAS
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_AND_DUAL_FEAS
        MathOptInterface.NearFeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MathOptInterface.InfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
        MathOptInterface.NearlyInfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MathOptInterface.ReductionCertificate
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_INTEGER_OPTIMAL
        MathOptInterface.NearlyFeasiblePoint
    else
        MathOptInterface.UnknownResultStatus
    end
end

function MathOptInterface.cangetattribute(m::MosekModel,attr::MathOptInterface.DualStatus)
    if attr.N < 0 || attr.N > length(m.solutions)
        false
    else
        solsta = m.solutions[attr.N].solsta
        if     solsta == MSK_SOL_STA_UNKNOWN true
        elseif solsta == MSK_SOL_STA_OPTIMAL true
        elseif solsta == MSK_SOL_STA_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL true
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER true
        elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER false
        elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL false
        elseif solsta == MSK_SOL_STA_NEAR_INTEGER_OPTIMAL false
        else false
        end
    end
end
function MathOptInterface.getattribute(m::MosekModel,attr::MathOptInterface.DualStatus)
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL
        MathOptInterface.NearlyFeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_FEAS
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_FEAS
        MathOptInterface.NearlyFeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_AND_DUAL_FEAS
        MathOptInterface.NearFeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MathOptInterface.InfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
        MathOptInterface.NearlyInfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MathOptInterface.ReductionCertificate
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MathOptInterface.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MathOptInterface.FeasiblePoint
    elseif solsta == MSK_SOL_STA_NEAR_INTEGER_OPTIMAL
        MathOptInterface.NearlyFeasiblePoint
    else
        MathOptInterface.UnknownResultStatus
    end
end


text:
if Sys.islinux()
    exit() # using the jll package now
end

const PROJECT_DIR = (@__DIR__) |> dirname
const TORCH_LIB_DIR = joinpath(PROJECT_DIR, "csrc/libtorch/lib")
const TORCH_LIB_BUILD_DIR = joinpath(PROJECT_DIR, "deps/lib")
const JULIA_THC_GENERATOR = joinpath(PROJECT_DIR, "src/thc/thc-generator.jl")

function build_locally()
    LIBTORCH_URL = if Sys.islinux()
        "https://download.pytorch.org/libtorch/cpu/libtorch-cxx11-abi-shared-with-deps-1.7.1%2Bcpu.zip"
    elseif Sys.isapple()
        "https://download.pytorch.org/libtorch/cpu/libtorch-macos-1.7.1.zip"
        # elseif Sys.iswindows()
        #     "https://download.pytorch.org/libtorch/cpu/libtorch-win-shared-with-deps-1.4.0.zip"
    else
        error("Your OS $(Sys.MACHINE) is not supported.")
    end

    if !isdir(TORCH_LIB_DIR)
        zipfile = download(LIBTORCH_URL)
        cd(joinpath(PROJECT_DIR, "csrc")) do
            run(`unzip $(zipfile)`)
        end
        isdir(TORCH_LIB_DIR) || error("Failed to get libtorch.")
    end

    isdir(TORCH_LIB_BUILD_DIR) || mkdir(TORCH_LIB_BUILD_DIR)

    cd(TORCH_LIB_BUILD_DIR) do
        cmd_cmake = `cmake -DCMAKE_PREFIX_PATH=$(joinpath(PROJECT_DIR, "csrc/libtorch")) ../../csrc`
        run(cmd_cmake)
        run(`make torch_capi`)
    end
end

build_locally()

if !isempty(get(ENV, "THARRAYS_DEV", ""))
    JULIA_EXE = joinpath(Sys.BINDIR, "julia")
    run(`$(JULIA_EXE) $(JULIA_THC_GENERATOR)`)
end


text:
using ADI: find_angles, compute_pa_thresh

@testset "helpers" begin
    angles = 10:10:90
    fwhm = 4.7
    r = 15
    idx = 6

    full_inds = axes(angles, 1)
    @test find_angles(angles, idx, 3) == setdiff(full_inds, idx)
    @test find_angles(angles, idx, 3, limit=4) == [4, 5, 7, 8]

    @test compute_pa_thresh(angles, r, fwhm, 1.5) ≈ 26.449102384427057
    @test compute_pa_thresh(angles, r, fwhm, 5) ≈ 36
end

cube, angles = BetaPictoris[:cube, :pa]
cube .-= minimum(cube) # rescale so min is > 0

@testset "framewise - $(nameof(typeof(alg)))" for alg in [PCA(10), NMF(3), Classic()]
    fr_alg = Framewise(alg)
    S = reconstruct(fr_alg, cube; angles=angles, fwhm=4.7, r=15)
    @test size(S) == size(cube)

    @test fr_alg(cube, angles; r=15, fwhm=4.7) ≈ collapse!(cube .- S, angles)
end


text:
#InspectDR: Edit plot layout within Gtk plot layer
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const LAYOUTHANDLE_ALLOWANCE=3 #Distance (px): Allowance for marker hit test
const EDITINGLAYOUT_LINESTYLE = LineStyle(:dash, 3.0, COLOR_BLACK)
const EDITINGLAYOUT_TEXTBOX_AREAATTR = AreaAttributes(
	line=LineAttributes(:solid, 2, COLOR_BLACK),
	fillcolor=ARGB32(1, 1, 1, 0.8)
)
const LAYOUTELEM_LIST = [
	:valloc_top, :valloc_mid, :valloc_bottom,
	:halloc_left, :halloc_right,
	:halloc_legend,
]

#==Types
===============================================================================#
#User input states
mutable struct ISWaitingToEditLayout <: InputState
end
mutable struct ISEditingLayout <: InputState #Mouse-move
	elem::Symbol
	v0::Float64 #Initial value (before editing)
	p0::Point2D #Start
	pmouse::Point2D #move
end
function ISEditingLayout(elem::Symbol, lyt, p0::Point2D)
	v0 = lyt[elem]
	ISEditingLayout(elem, v0, p0, p0)
end

#Info used to draw control elements for layout elements
struct LayoutControlInfo
	elem::Symbol
	ctrlline::Symbol #Where is the control line?: xmin, xmax, ymin or ymax
	bb::BoundingBox
	visible::Bool
end
LayoutControlInfo(elem, ctrlline, bb) = LayoutControlInfo(elem, ctrlline, bb, true)

#==Helper functions
===============================================================================#
isvert(lci::LayoutControlInfo) = lci.ctrlline == :ymin || lci.ctrlline == :ymax

function _linepoints(lci::LayoutControlInfo)
	bb = lci.bb
	xmin, xmax = lci.bb.xmin, lci.bb.xmax
	ymin, ymax = lci.bb.ymin, lci.bb.ymax

	if :xmin == lci.ctrlline
		xmax = xmin
	elseif :xmax == lci.ctrlline
		xmin = xmax
	elseif :ymin == lci.ctrlline
		ymax = ymin
	elseif :ymax == lci.ctrlline
		ymin = ymax
	else
#		xmax = xmin
#		ymax = ymin
	end

	return (Point2D(xmin, ymin), Point2D(xmax, ymax))
end

function _computenewval(s::ISEditingLayout, lci::LayoutControlInfo)
	v = s.v0
	Δ = s.pmouse - s.p0
	if :xmin == lci.ctrlline
		v -= Δ.x
	elseif :xmax == lci.ctrlline
		v += Δ.x
	elseif :ymin == lci.ctrlline
		v -= Δ.y
	elseif :ymax == lci.ctrlline
		v += Δ.y
	end
	return round(v)
end

#==Layout control elements
===============================================================================#
function LayoutControlInfo(::DS{:valloc_top}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:valloc_top]
	return LayoutControlInfo(:valloc_top, :ymax,
		BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymin+v),
	)
end
function LayoutControlInfo(::DS{:valloc_mid}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:valloc_mid]
	visible = false
	return LayoutControlInfo(:valloc_mid, :ymax,
		BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymin+v),
		visible
	)
end
function LayoutControlInfo(::DS{:valloc_bottom}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:valloc_bottom]
	return LayoutControlInfo(:valloc_bottom, :ymin,
		BoundingBox(bb.xmin, bb.xmax, bb.ymax-v, bb.ymax),
	)
end
function LayoutControlInfo(::DS{:halloc_left}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:halloc_left]
	return LayoutControlInfo(:halloc_left, :xmax,
		BoundingBox(bb.xmin, bb.xmin+v, bb.ymin, bb.ymax),
	)
end
function LayoutControlInfo(::DS{:halloc_right}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:halloc_right]
	xmax = bb.xmax# - lyt[:halloc_legend]
	return LayoutControlInfo(:halloc_right, :xmin,
		BoundingBox(xmax-v, xmax, bb.ymin, bb.ymax),
		!lyt[:enable_legend]
	)
end
function LayoutControlInfo(::DS{:halloc_legend}, rplot::RPlot2D, lyt)
	bb = rplot.bb;	v = lyt[:halloc_legend]
	return LayoutControlInfo(:halloc_legend, :xmin,
		BoundingBox(bb.xmax-v, bb.xmax, bb.ymin, bb.ymax),
		lyt[:enable_legend]
	)
end

#==Hit test
===============================================================================#
function hittest(lci::LayoutControlInfo, x::Float64, y::Float64)
	if !lci.visible; return false; end
	v = getfield(lci.bb, lci.ctrlline)
	if isvert(lci)
		return abs(y-v) <= LAYOUTHANDLE_ALLOWANCE
	else
		return abs(x-v) <= LAYOUTHANDLE_ALLOWANCE
	end
end

#==Draw functions
===============================================================================#
function editlyt_drawelem(ctx::CairoContext, rplot::RPlot2D, lyt, elem::Symbol)
	lci = LayoutControlInfo(DS(elem), rplot, lyt)
	if !lci.visible; return; end
	p1, p2 = _linepoints(lci)

	Cairo.save(ctx) #-----
	setlinestyle(ctx, EDITINGLAYOUT_LINESTYLE)
	drawline(ctx, p1, p2)
	Cairo.stroke(ctx)
	Cairo.restore(ctx) #-----
	return
end

function editlyt_drawall(ctx::CairoContext, rplot::RPlot2D, lyt)
	for elem in LAYOUTELEM_LIST
		editlyt_drawelem(ctx, rplot, lyt, elem)
	end
end

function editlyt_draweleminfo(ctx::CairoContext, lyt, elem::Symbol)
	v = lyt[elem]
	msg = ":$elem = $v"

	Cairo.save(ctx) #-----
	font = lyt[:font_annotation]
	setfont(ctx, font)
	w, h = textextents_wh(ctx, msg)
	Δ = h * 0.2
	bbann = BoundingBox(0, w+Δ, 0, h+Δ)
	ctr = Point2D(bbann.xmax/2, bbann.ymax/2)
	drawrectangle(ctx, bbann, EDITINGLAYOUT_TEXTBOX_AREAATTR)
	render(ctx, msg, ctr, font, align=ALIGN_HCENTER|ALIGN_VCENTER)
	Cairo.stroke(ctx)
	Cairo.restore(ctx) #-----
end

function drawoverlay(s::ISWaitingToEditLayout, ctx::CairoContext, rplot::RPlot2D, lyt)
	editlyt_drawall(ctx, rplot, lyt)
end
function drawoverlay(s::ISEditingLayout, ctx::CairoContext, rplot::RPlot2D, lyt)
	editlyt_drawelem(ctx, rplot, lyt, s.elem)
	editlyt_draweleminfo(ctx, lyt, s.elem)
end

#==State-dependent functions
===============================================================================#
function hittest(s::ISWaitingToEditLayout, pwidget::PlotWidget, x::Float64, y::Float64)
	lyt = pwidget.src.layout
	for elem in LAYOUTELEM_LIST
		lci = LayoutControlInfo(DS(elem), pwidget.rplot, lyt)
		if hittest(lci, x, y)
			return lci
		end
	end
	return nothing
end

function editlyt_hover(s::ISWaitingToEditLayout, pwidget::PlotWidget, x::Float64, y::Float64)
	lci = hittest(s, pwidget, x, y)
	if nothing == lci
		gdk_window_set_cursor(pwidget.canvas, CURSOR_DEFAULT)
		return
	end

	csr = isvert(lci) ? CURSOR_ROWRESIZE : CURSOR_COLRESIZE
	gdk_window_set_cursor(pwidget.canvas, csr)
	return lci
end

function editlyt_mousedown(s::ISWaitingToEditLayout, pwidget::PlotWidget, x::Float64, y::Float64)
	lci = editlyt_hover(s, pwidget, x, y) #Get right cursor, etc
	if nothing == lci
		return
	end
	pwidget.state = ISEditingLayout(lci.elem, pwidget.src.layout, Point2D(x,y))
	Gtk.draw(pwidget.canvas)
end

#Exit ISWaitingToEditLayout editing state
function editlyt_exit(s::ISWaitingToEditLayout, pwidget::PlotWidget)
	setstate_normal(pwidget)
	refresh(pwidget, refreshdata=true) #Proper refresh when exiting
end

#Complete a layout mouse-move operation:
function editlyt_mousemove(s::ISEditingLayout, pwidget::PlotWidget, x::Float64, y::Float64)
	s.pmouse = Point2D(x, y)
	lci = LayoutControlInfo(DS(s.elem), pwidget.rplot, pwidget.src.layout)
	pwidget.src.layout[s.elem] = _computenewval(s, lci)
	Gtk.draw(pwidget.canvas)
	return
end

function editlyt_mouseup(s::ISEditingLayout, pwidget::PlotWidget, x::Float64, y::Float64)
	#Already modified value.  Just revert to wait state:
	pwidget.state = ISWaitingToEditLayout()
	refresh(pwidget, refreshdata=true)
end

function editlyt_cancel(s::ISEditingLayout, pwidget::PlotWidget)
	pwidget.state = ISWaitingToEditLayout() #before resetting layout value
	pwidget.src.layout[s.elem] = s.v0
	Gtk.draw(pwidget.canvas)
end

#==Event handlers for state: ISWaitingToEditLayout
===============================================================================#
function handleevent_keypress(s::ISWaitingToEditLayout, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval
		editlyt_exit(s, pwidget)
	end
end
function handleevent_mousemove(s::ISWaitingToEditLayout, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	editlyt_hover(s, pwidget, event.x, event.y)
end
function handleevent_mousepress(s::ISWaitingToEditLayout, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	set_focus(pwidget) #In case not in focus

	if 1==event.button
		editlyt_mousedown(s, pwidget, event.x, event.y)
	end
end

#==Event handlers for state: ISEditingLayout (change properties)
===============================================================================#
function handleevent_keypress(s::ISEditingLayout, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval
		editlyt_cancel(s, pwidget)
	end
end

function handleevent_mouserelease(s::ISEditingLayout, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	if 1==event.button
		editlyt_mouseup(s, pwidget, event.x, event.y)
	end
end
function handleevent_mousemove(s::ISEditingLayout, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	#handleevent_plothover(pwidget, event)
	editlyt_mousemove(s, pwidget, event.x, event.y)
end
function handleevent_mousepress(s::ISEditingLayout, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	set_focus(pwidget) #In case not in focus
end

#==Switch to layout-editing state
===============================================================================#
function setstate_editlayout(pwidget::PlotWidget)
	pwidget.state = ISWaitingToEditLayout()
end

#Last line


text:
using Documenter, Andes

makedocs(modules = [Andes], sitename = "Andes.jl")

deploydocs(repo = "github.com/cuihantao/Andes.jl.git")


text:
# A method of rand to sample from a triangular distribution
# https://en.wikipedia.org/wiki/Triangular_distribution

import Distributions.TriangularDist
import Distributions.params
import Distributions.rand
#=
function rand(rng::AbstractRNG, d::TriangularDist)
    (a, b, c) = params(d)
    b_m_a = b - a
    u = rand(rng)
    if b_m_a * u < (c - a)
        return a + sqrt(u * b_m_a * (c - a))
    else
        return b - sqrt((1 - u) * b_m_a * (b - c))
    end
end
=#

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
import UUIDs

"""Recursively validates that the object and fields have UUIDs."""
function validate_uuids(obj::T) where {T}
    if !(obj isa Component)
        return true
    end

    result = true
    if !(IS.get_uuid(obj) isa Base.UUID)
        result = false
        @error "object does not have a UUID" obj
    end

    for fieldname in fieldnames(T)
        val = getfield(obj, fieldname)
        if !validate_uuids(val)
            result = false
        end
    end

    return result
end

function validate_uuids(obj::T) where {T <: AbstractArray}
    result = true
    for elem in obj
        if !validate_uuids(elem)
            result = false
        end
    end

    return result
end

function validate_uuids(obj::T) where {T <: AbstractDict}
    result = true
    for elem in values(obj)
        if !validate_uuids(elem)
            result = false
        end
    end
    return result
end

@testset "Test internal values" begin
    sys_rts = create_rts_system()
    @test validate_uuids(sys_rts)
end


text:
function _nancov(x::AbstractVector, y::AbstractVector, corrected::Bool, μᵪ::Number, μᵧ::Number)
    # Calculate covariance
    σᵪᵧ = ∅ = zero(promote_type(typeof(μᵪ), typeof(μᵧ), Int))
    n = 0
    @turbo for i ∈ indices((x,y))
            δᵪ = x[i] - μᵪ
            δᵧ = y[i] - μᵧ
            δ² = δᵪ * δᵧ
            notnan = δ²==δ²
            n += notnan
            σᵪᵧ += ifelse(notnan, δ², ∅)
    end
    σᵪᵧ = σᵪᵧ / (n-corrected)
    return σᵪᵧ
end

"""
```julia
nancov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
```
Compute the covariance between the vectors `x` and `y`.
As `Statistics.cov`, but ignoring `NaN`s.

If `corrected` is `true` as is the default, _Bessel's correction_ will be applied,
such that the sum is scaled by `n-1` rather than `n`, where `n = length(x)`.
"""
function nancov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
    # Check lengths
    nᵪ = length(x)
    nᵧ = length(y)
    @assert nᵪ == nᵧ

    μᵪ = _nanmean(x,:)
    μᵧ = _nanmean(y,:)
    σᵪᵧ = _nancov(x, y, corrected, μᵪ, μᵧ)
    return σᵪᵧ
end

"""
```julia
nancov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
```
Compute the covariance matrix of the matrix `X`, along dimension `dims`.
As `Statistics.cov`, but ignoring `NaN`s.

If `corrected` is `true` as is the default, _Bessel's correction_ will be applied,
such that the sum is scaled by `n-1` rather than `n`, where `n = length(x)`.
"""
function nancov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Σ = similar(X, Tₒ, (m, m))
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means for each column
        μ = ntuple(m) do d
            nanmean(view(X,:,d))
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _nancov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    elseif dims == 2
        # Precalculate means for each row
        μ = ntuple(m) do d
            nanmean(view(X,d,:))
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _nancov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    else
        throw("Dimension not in range")
    end
    return Σ
end
export nancov

"""
```julia
nancor(x::AbstractVector, y::AbstractVector)
```
Compute the (Pearson's product-moment) correlation between the vectors `x` and `y`.
As `Statistics.cor`, but ignoring `NaN`s.

Equivalent to `nancov(x,y) / (nanstd(x) * nanstd(y))`.
"""
function nancor(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
    # Check lengths
    nᵪ = length(x)
    nᵧ = length(y)
    @assert nᵪ == nᵧ

    μᵪ = nanmean(x)
    μᵧ = nanmean(y)
    σᵪ = nanstd(x, mean=μᵪ, corrected=corrected)
    σᵧ = nanstd(y, mean=μᵧ, corrected=corrected)
    σᵪᵧ = _nancov(x, y, corrected, μᵪ, μᵧ)
    ρᵪᵧ = σᵪᵧ / (σᵪ * σᵧ)

    return ρᵪᵧ
end

"""
```julia
nancor(X::AbstractMatrix; dims::Int=1)
```
Compute the (Pearson's product-moment) correlation matrix of the matrix `X`,
along dimension `dims`. As `Statistics.cor`, but ignoring `NaN`s.
"""
function nancor(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Ρ = similar(X, Tₒ, (m, m))
    # Diagonal must be unity
    @inbounds for i = 1:m
        Ρ[i,i] = one(Tₒ)
    end
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            nanmean(view(X,:,d))
        end
        σ = ntuple(m) do d
            nanstd(view(X,:,d), mean=μ[d], corrected=corrected)
        end
        # Fill off-diagonals symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _nancov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    elseif dims == 2
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            nanmean(view(X,d,:))
        end
        σ = ntuple(m) do d
            nanstd(view(X,d,:), mean=μ[d], corrected=corrected)
        end
        @inbounds for i = 1:m
            for j = 1:i-1
                σᵢⱼ = _nancov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    else
        throw("Dimension not in range")
    end
    return Ρ
end
export nancor


text:
# this file contains code for the simulated annealing TSP heuristic

"""
.. simulated_annealing(distmat; ...) ..

Use a simulated annealing strategy to return a closed tour. The temperature
decays exponentially from init_temp to final_temp.

Optional arguments:

- steps: number of steps to take; defaults to 50n^2 where n is number of cities

- num_starts: number of times to run the simulated annealing algorithm, each time
	starting with a random path. Defaults to 1.

- init_temp: initial temperature which controls initial chance of accepting an
	inferior tour.

- final_temp: final temperature which controls final chance of accepting an
	inferior tour; lower values roughly correspond to a longer period of 2-opt.

- init_path: path to start the annealing from. A Nullable{Vector{Int}}. An empty
	Nullable corresponds to picking a random path; if the Nullable contains a value
	then this path will be used. Defaults to a random path.
"""
function simulated_annealing{T <: Real}(distmat::Matrix{T}; steps = 50*length(distmat),
										num_starts = 1,
										init_temp = exp(8), final_temp = exp(-6.5),
										init_path::Nullable{Vector{Int}} = Nullable{Vector{Int}}())

	# check inputs
	check_square(distmat, "Must pass a square distance matrix to simulated_annealing.")

	# cooling rate: we multiply by a constant mult each step
	cool_rate = (final_temp / init_temp)^(1 / (steps - 1))

	# do SA with a single starting location
	function sahelper()
		temp = init_temp / cool_rate # divide by cool_rate so when we first multiply we get init_temp
		n = size(distmat, 1)
		path = isnull(init_path) ? randpath(n) : copy(get(init_path))
		cost_cur = pathcost(distmat, path)

		for i in 1:steps
			temp *= cool_rate

			# take a step
			# keep first and last cities fixed
			first, last = rand(2:n), rand(2:n)
			if first > last
				first, last = last, first
			end
			cost_other = pathcost_rev(distmat, path, first, last)
			@fastmath accept = cost_other < cost_cur ? true : rand() < exp((cost_cur - cost_other) / temp)
			# should we accept?
			if accept
				reverse!(path, first, last)
				cost_cur = cost_other
			end
		end

		return path, cost_cur
	end

	path, cost = sahelper()
	for _ in 2:num_starts
		otherpath, othercost = sahelper()
		if othercost < cost
			cost = othercost
			path = otherpath
		end
	end

	return path, cost
end

function randpath(n)
	path = 1:n |> collect |> shuffle
	push!(path, path[1]) # loop
	return path
end


text:
function mse_loss(θ, y, predict)
    ŷ = predict(θ)
    loss = sum(abs2, y .- ŷ), ŷ
end

function polar_loss(θ, y, w, predict)
    ŷ = predict(θ)
    N = size(ŷ, 2)
    loss = sum(w[1] * [normed_ld(y[:,i], ŷ[:,i]) for i in 1:N] +
               w[2] * [cosine_distance(y[:,i], ŷ[:,i]) for i in 1:N])
    return loss, ŷ
end

normed_ld(a,b) = abs(norm(a)-norm(b))/(norm(a)+norm(b))

cosine_similarity(a, b) = dot(a,b) / (norm(a) * norm(b))
cosine_distance(a, b) = (1 - cosine_similarity(a, b))/2

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
using MCMCDiagnostics
using Test

using Random
using Statistics: mean, var

# consistent testing
Random.seed!(UInt32[0x3d50884f, 0xd6560f94, 0x6c04ab37, 0xb1c52878])

@testset "IID samples" begin
    chains = [randn(1000) for _ in 1:10]
    @test 9900 ≤ effective_sample_size(vcat(chains...)) ≤ 10100
    @test 1 ≤ potential_scale_reduction(chains...) ≤ 1.001
end

@testset "IID dispersed ragged PSRF" begin
    μs = -10:1:10
    chains = [randn(1000) * rand(1:10) .+ μ for μ in μs]
    B = var(μs)
    W = mean(var.(chains))
    expected_R̂ = √(1 + B / W)
    R̂ = potential_scale_reduction(chains...)
    @test 0.99 ≤ (R̂ / expected_R̂) ≤ 1.01
end

@testset "AR(1)" begin
    "Simulate an AR(1) process (coefficient ρ, N draws, std σ)."
    function simulate_ar1(ρ, N; σ = 1.0)
        x = Vector{Float64}(undef, N)
        z = σ * randn() / √(1 - ρ^2)
        for i in 1:N
            z = ρ*z + randn()*σ
            x[i] = z
        end
        x
    end

    "Theoretical ESS."
    ar1_ess_factor(ρ) = 1/(1 + 2*ρ/(1-ρ))

    """
    Estimate of effective sample size for an AR(1) process with
    coefficient ρ, divided by the theoretical value.  Should be
    invariant to σ, and around 1.
    """
    function rel(ρ, N)
        effective_sample_size(simulate_ar1(ρ, N, σ=randn()^2+0.5))/(ar1_ess_factor(ρ)*N)
    end

    for ρ in range(0.1; length = 10, stop = 0.4)
        @test 0.95 ≤ rel(ρ,100000) ≤ 1.05
    end

    for ρ in range(0.5; length = 10, stop = 0.9)
        @test 0.95 ≤ rel(ρ,1000000) ≤ 1.05
    end

    # testing cap, since for ρ < 0, ESS ≥ N
    for ρ in range(-0.5; length = 10, stop = -0.9)
        effective_sample_size(simulate_ar1(ρ, 1000, σ=randn()^2+0.5)) == 1000
    end
end


text:
function print_summary(result::TSPDroneResult)
    report = generate_summary(result)
    println(report)
end

fmt(n) = n #round(n, digits=5)

function generate_summary(result::TSPDroneResult)
    combined_nodes = intersect(result.truck_route, result.drone_route)
    operations = [(combined_nodes[i], combined_nodes[i+1]) for i in 1:length(combined_nodes)-1]

    op_costs = Float64[]
    io = IOBuffer()

    for i in 1:length(operations)
        op = operations[i]
        truck_length, truck_sub_route =
            operation_length(op[1], op[2], result.truck_route, result.Ct)
        drone_length, drone_sub_route =
            operation_length(op[1], op[2], result.drone_route, result.Cd)

        op_length = max(truck_length, drone_length)
        push!(op_costs, op_length)

        println(io, "Operation #$(i):")
        println(io, "  - Truck        = $(fmt(truck_length)) : $(truck_sub_route)")
        println(io, "  - Drone        = $(fmt(drone_length)) : $(drone_sub_route)")
        println(io, "  - Length       = $(fmt(op_length))")
        if drone_length > result.flying_range
            write(io, "  - flying_range : ")
            printstyled(io, "FAIL\n", color=:red)
        end

    end

    @assert sum(op_costs) == result.total_cost

    write(io, "Total Cost = $(fmt(sum(op_costs)))")

    report = String(take!(io))

    return report
end

function operation_length(origin, destination, route, cost_mtx)
    o_idx = findfirst(x -> x == origin, route)
    d_idx = findfirst(x -> x == destination, route)
    sub_route = route[o_idx:d_idx]

    length = 0.0
    for i in o_idx:d_idx-1
        length += cost_mtx[route[i], route[i+1]]
    end

    return length, sub_route
end


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################

using Coverage, Base.Test, Compat
using Compat.String

# test our filename matching. These aren't exported functions but it's probably
# a good idea to have explicit tests for them, as they're used to match files
# that get deleted

@test Coverage.iscovfile("test.jl.cov")
@test Coverage.iscovfile("test.jl.2934.cov")
@test Coverage.iscovfile("/home/somebody/test.jl.2934.cov")
@test !Coverage.iscovfile("test.ji.2934.cov")
@test !Coverage.iscovfile("test.jl.2934.cove")
@test !Coverage.iscovfile("test.jicov")
@test !Coverage.iscovfile("test.c.cov")
@test Coverage.iscovfile("test.jl.cov", "test.jl")
@test !Coverage.iscovfile("test.jl.cov", "other.jl")
@test Coverage.iscovfile("test.jl.8392.cov", "test.jl")
@test Coverage.iscovfile("/somedir/test.jl.8392.cov", "/somedir/test.jl")
@test !Coverage.iscovfile("/otherdir/test.jl.cov", "/somedir/test.jl")

cd(Pkg.dir("Coverage")) do
    datadir = joinpath("test", "data")
    # Process a saved set of coverage data...
    r = process_file(joinpath(datadir,"Coverage.jl"))
    # ... and memory data
    analyze_malloc(datadir)
    lcov = IOBuffer()
    # we only have a single file, but we want to test on the Vector of file results
    LCOV.write(lcov, FileCoverage[r])
    open(joinpath(datadir, "expected.info")) do f
        @test takebuf_string(lcov) == readstring(f)
    end

    # Test a file from scratch
    srcname = joinpath("test", "data","testparser.jl")
    covname = srcname*".cov"
    # clean out any previous coverage files. Don't use clean_folder because we
    # need to preserve the pre-baked coverage file Coverage.jl.cov
    clean_file(srcname)
    cmdstr = "include(\"$srcname\"); using Base.Test; @test f2(2) == 4"
    run(`julia --code-coverage=user -e $cmdstr`)
    r = process_file(srcname, datadir)
    # The next one is the correct one, but julia & JuliaParser don't insert a line number after the 1-line @doc -> test
    # See https://github.com/JuliaLang/julia/issues/9663 (when this is fixed, can uncomment the next line on julia 0.4)
    target = Union{Int64,Void}[nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing]
    #target = Union{Int64,Void}[nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing]
    @test r.coverage[1:length(target)] == target

    covtarget = (sum(x->x != nothing && x > 0, target), sum(x->x != nothing, target))
    @test get_summary(r) == covtarget
    @test get_summary(process_folder(datadir)) != covtarget

    #json_data = Codecov.build_json_data(Codecov.process_folder("data"))
    #@test typeof(json_data["coverage"]["data/Coverage.jl"]) == Array{Union{Int64,Void},1}
    open("fakefile",true,true,true,false,false)
    @test isempty(Coverage.process_cov("fakefile",datadir))
    rm("fakefile")
end

######################
# codecovio.jl tests #
######################

"""
extracts the api URL from stdout in a codecov.io submit call
very helpful for testing codecovio.jl
"""
function extract_codecov_url(fun)
    originalSTDOUT = STDOUT

    (outRead, outWrite) = redirect_stdout()

    fun()

    close(outWrite)

    data = String(readavailable(outRead))

    close(outRead)
    redirect_stdout(originalSTDOUT)

    lines = split(data, "\n")

    url = "None"
    get_next = false
    for line in lines
        if get_next
            url = line
            get_next = false
        end
        if contains(line, "Codecov.io API URL")
            get_next = true
        end
    end

    #println("url: $(url)")
    @assert url != "None" "unable to find codecov api url in stdout, check for changes in codecovio.jl"
    return url
end

# empty file coverage for testing
fcs = FileCoverage[]

# setup base system ENV vars for testing
withenv(
    "CODECOV_URL" => nothing,
    "CODECOV_TOKEN" => nothing,
    "TRAVIS_BRANCH" => nothing,
    "TRAVIS_COMMIT" => nothing,
    "TRAVIS_PULL_REQUEST" => nothing,
    "TRAVIS_JOB_ID" => nothing,
    "TRAVIS_REPO_SLUG" => nothing,
    "TRAVIS_JOB_NUMBER" => nothing
    ) do

    # test local submission process

    # default values
    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
    @test contains(codecov_url, "codecov.io")
    @test contains(codecov_url, "commit")
    @test contains(codecov_url, "branch")
    @test !contains(codecov_url, "service")

    # default values in depreciated call
    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs) )
    @test contains(codecov_url, "codecov.io")
    @test contains(codecov_url, "commit")
    @test contains(codecov_url, "branch")
    @test !contains(codecov_url, "service")

    # env var url override
    withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
        @test contains(codecov_url, "enterprise-codecov-1.com")
        @test contains(codecov_url, "commit")
        @test contains(codecov_url, "branch")
        @test !contains(codecov_url, "service")

        # function argument url override
        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
        @test contains(codecov_url, "enterprise-codecov-2.com")
        @test contains(codecov_url, "commit")
        @test contains(codecov_url, "branch")
        @test !contains(codecov_url, "service")

        # env var token
        withenv( "CODECOV_TOKEN" => "token_name_1" ) do

            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
            @test contains(codecov_url, "enterprise-codecov-1.com")
            @test contains(codecov_url, "token=token_name_1")
            @test !contains(codecov_url, "service")

            # function argument token url override
            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true, token="token_name_2") )
            @test contains(codecov_url, "enterprise-codecov-1.com")
            @test contains(codecov_url, "token=token_name_2")
            @test !contains(codecov_url, "service")
        end
    end

    # test travis-ci submission process

    #setup travis env
    withenv(
        "TRAVIS_BRANCH" => "t_branch",
        "TRAVIS_COMMIT" => "t_commit",
        "TRAVIS_PULL_REQUEST" => "t_pr",
        "TRAVIS_JOB_ID" => "t_job_id",
        "TRAVIS_REPO_SLUG" => "t_slug",
        "TRAVIS_JOB_NUMBER" => "t_job_num"
        ) do

        # default values
        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
        @test contains(codecov_url, "codecov.io")
        @test contains(codecov_url, "service=travis-org")
        @test contains(codecov_url, "branch=t_branch")
        @test contains(codecov_url, "commit=t_commit")
        @test contains(codecov_url, "pull_request=t_pr")
        @test contains(codecov_url, "job=t_job_id")
        @test contains(codecov_url, "slug=t_slug")
        @test contains(codecov_url, "build=t_job_num")

        # env var url override
        withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
            @test contains(codecov_url, "enterprise-codecov-1.com")
            @test contains(codecov_url, "service=travis-org")
            @test contains(codecov_url, "branch=t_branch")
            @test contains(codecov_url, "commit=t_commit")
            @test contains(codecov_url, "pull_request=t_pr")
            @test contains(codecov_url, "job=t_job_id")
            @test contains(codecov_url, "slug=t_slug")
            @test contains(codecov_url, "build=t_job_num")

            # function argument url override
            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
            @test contains(codecov_url, "enterprise-codecov-2.com")
            @test contains(codecov_url, "service=travis-org")
            @test contains(codecov_url, "branch=t_branch")
            @test contains(codecov_url, "commit=t_commit")
            @test contains(codecov_url, "pull_request=t_pr")
            @test contains(codecov_url, "job=t_job_id")
            @test contains(codecov_url, "slug=t_slug")
            @test contains(codecov_url, "build=t_job_num")

            # env var token
            withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                @test contains(codecov_url, "enterprise-codecov-1.com")
                @test contains(codecov_url, "token=token_name_1")
                @test contains(codecov_url, "branch=t_branch")
                @test contains(codecov_url, "commit=t_commit")
                @test contains(codecov_url, "pull_request=t_pr")
                @test contains(codecov_url, "job=t_job_id")
                @test contains(codecov_url, "slug=t_slug")
                @test contains(codecov_url, "build=t_job_num")

                # function argument token url override
                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, token="token_name_2") )
                @test contains(codecov_url, "enterprise-codecov-1.com")
                @test contains(codecov_url, "token=token_name_2")
                @test contains(codecov_url, "branch=t_branch")
                @test contains(codecov_url, "commit=t_commit")
                @test contains(codecov_url, "pull_request=t_pr")
                @test contains(codecov_url, "job=t_job_id")
                @test contains(codecov_url, "slug=t_slug")
                @test contains(codecov_url, "build=t_job_num")
            end
        end
    end

end


text:
using Altro
using TrajectoryOptimization
using StaticArrays
using BenchmarkTools
using RobotZoo
using RobotDynamics
using StaticArrays, LinearAlgebra
using ForwardDiff
using Test
const TO = TrajectoryOptimization
const RD = RobotDynamics

## Test Functions
function Πsoc(x)
    v = x[1:end-1]
    s = x[end]
    a = norm(v)
    if a <= -s
        return zero(x)
    elseif a <= s
        return x
    elseif a >= abs(s)
        x̄ = append!(v, a)
        return 0.5*(1 + s/a) * x̄
    end
    throw(ErrorException("Invalid second-order cone"))
end
function in_soc(x)
    v = x[1:end-1]
    s = x[end]
    a = norm(v)
	a <= s
end

function auglag(f,g, x,z,μ)
    # ceq = h(x)
    L0 = f(x) #+ y'ceq + 0.5*μ*ceq'ceq

    cones = g(x)
    p = length(cones)
    for i = 1:p
        cone = cones[i]
		proj = Πsoc(z[i] - μ*cones[i])
		pen = proj'proj - z[i]'z[i]
        L0 += 1 / (2μ) * pen
    end
    return L0
end

function DI_objective(Q,R,Qf,xf,dt)
    g = Float64[]
    h = Float64[]
    c = 0.0
    for k = 1:N-1
        append!(g, -Q*xf)
        append!(g, zeros(m))
        append!(h, Q.diag)
        append!(h, R.diag)
        c += 0.5*xf'Q*xf
	end
	g .*= dt
	h .*= dt
	c *= dt
    append!(g, -Qf*xf)
    # append!(g, zeros(m))
    append!(h, Qf.diag)
    # append!(h, R.diag)
    c += 0.5*xf'Qf*xf
    H = Diagonal(h)
    return H,g,c
end

function DI_dynamics(x0, N; dt=0.1, gravity::Bool=false)
    n = length(x0)
    m = n ÷ 2
    D = m
    NN = N*(n+m)

    #  Continuous dynamics
    Ac = zeros(n,n)
    Bc = zeros(n,m)
    for i = 1:D
        Ac[i,D+i] = 1
        Bc[D+i,i] = 1
    end

    # Euler integration
    Ad = I + Ac*dt
    Bd = Bc*dt

    iA = 1:n
    jA = iA
    jB = n .+ (1:m)
    jA2 = jA .+ (n+m)

    A = zeros(N*n, NN)
    b = zeros(N*n)
    A[iA,jA] .= I(n)
    b[iA] .= -x0
    iA = iA .+ n
    for k = 1:N-1
        A[iA,jA] .= Ad
        A[iA,jB] .= Bd
        A[iA,jA2] .= -I(n)
        if gravity
            b[iA[D]] = -9.81
        end

        # Advance inds
        iA = iA .+ n
        jA = jA .+ (n+m)
        jB = jB .+ (n+m)
        jA2 = jA2 .+ (n+m)
    end
    return A,b
end

function DI_cones(::Val{D},N) where D
    uinds = SVector{D}(1:D) .+ 2D
    qinds = [uinds .+ (k-1)*3D for k = 1:N-1]
end
DI_cones(D,N) = DI_cones(Val(D),N)  # WARNING: type unstable

## Setup
D,N = 2,11
n,m = 2D,D
dt = 0.1
tf = (N-1)*dt
x0 = @SVector fill(0., n)
xf = [(@SVector fill(1.,D)); (@SVector fill(0, D))]
Q = Diagonal(@SVector fill(1.0, n))
R = Diagonal(@SVector fill(1e-1, m))
Qf = (N-1)*Q * 100
qinds = DI_cones(D,N)
u_bnd = 6.0

## Define batch functions
H_obj,g_obj,c_obj = DI_objective(Q,R,Qf,xf,dt)
Adyn,bdyn = DI_dynamics(x0,N, dt=dt)

di_obj(x) = 0.5*x'H_obj*x + g_obj'x + c_obj
di_dyn(x) = Adyn*x + bdyn
di_soc(x) = [push(x[qi],u_bnd) for qi in qinds]

## Define using ALTRO
model = RobotZoo.DoubleIntegrator(D)

obj = LQRObjective(Q, R, Qf, xf, N)
cons = ConstraintList(n,m,N)
cone = TO.SecondOrderCone()
# cone = TO.Inequality()
u_bnd = 6.0
connorm = NormConstraint(n, m, u_bnd, cone, :control)
add_constraint!(cons, connorm, 1:N-1)
prob = Problem(model, obj, xf, tf, x0=x0, constraints=cons, integration=RD.Euler)
solver = Altro.AugmentedLagrangianSolver(prob)
conSet = TO.get_constraints(solver)
alobj = TO.get_objective(solver)

## Initialization
inds = LinearIndices(zeros(3D,N))
xinds = [SVector{n}(z[1:n]) for z in eachcol(inds)]
uinds = [SVector{m}(z[n+1:end]) for z in eachcol(inds)][1:end-1]
NN = 3D*N-D
P = (N-1)*m
x0 = rand(NN)
z = rand.(length.(di_soc(x0)))
μ = 2.4

X0 = [x0[xi] for xi in xinds]
U0 = [x0[ui] for ui in uinds]
Z0 = Traj(X0,U0,fill(dt,N))
initial_trajectory!(solver, Z0)

## Test objective and soc constraints
cost(obj, Z0) ≈ di_obj(x0)
conSet.convals[1].vals ≈ di_soc(x0)

# copy multipliers to ALconSet
cval = conSet.convals[1]
for i = 1:N-1
	cval.λ[i] .= z[i]
	cval.μ[i] .= μ
end

# test augmented Lagrangian value
LA(x) = auglag(di_obj, di_soc, x, z, μ)
cost(solver) ≈ LA(x0)
!(cost(solver) ≈ di_obj(x0))   # make sure the AL cost isn't the same as the normal cost

λbar = z .- μ .* di_soc(x0)
statuses = [TO.cone_status(TO.SecondOrderCone(), λ) for λ in λbar]
@test :below ∈ statuses

# E = TO.QuadraticObjective(n,m,N)
E = TO.CostExpansion(n,m,N)
TO.cost_expansion!(E, alobj, Z0)
grad = vcat([[e.q; e.r] for e in E]...)[1:end-m]
@test grad ≈ ForwardDiff.gradient(LA, x0)

hess_blocks = vcat([[e.Q, e.R] for e in E]...)
hess = cat(hess_blocks..., dims=(1,2))[1:end-m, 1:end-m]
@test hess ≈ ForwardDiff.hessian(LA, x0)

@test all(in_soc.(di_soc(x0)))
@test !all(in_soc.(z .- μ*di_soc(x0)))

# Make it go outside the cone
z = rand.(length.(di_soc(x0)))
z .*= 100
λbar = z .- μ .* di_soc(x0)
statuses = [TO.cone_status(TO.SecondOrderCone(), λ) for λ in λbar]
@test :outside ∈ statuses
@test :in ∈ statuses
for i = 1:N-1
	cval.λ[i] .= z[i]
	cval.μ[i] .= μ
end
LA(x) = auglag(di_obj, di_soc, x, z, μ)

# E = TO.QuadraticObjective(n,m,N)
E = TO.CostExpansion(n,m,N)
TO.cost_expansion!(E, alobj, Z0)
grad = vcat([[e.q; e.r] for e in E]...)[1:end-m]
@test grad ≈ ForwardDiff.gradient(LA, x0)

hess_blocks = vcat([[e.Q, e.R] for e in E]...)
hess = cat(hess_blocks..., dims=(1,2))[1:end-m, 1:end-m]
@test hess ≈ ForwardDiff.hessian(LA, x0)

## Solve it
solver = Altro.AugmentedLagrangianSolver(prob, verbose=0, show_summary=false,
	projected_newton=true, penalty_initial=100.0, penalty_scaling=50,
	cost_tolerance_intermediate=1e-1)
initial_controls!(solver, [rand(D) for k = 1:N])
solve!(solver)
norm.(controls(solver))
@test all(abs.(norm.(controls(solver))[[1,2,N-2,N-1]] .- u_bnd) .< 1e-6)

## Rocket landing problem
D,N = 3,101
n,m = 2D,D
dt = 0.1
tf = (N-1)*dt
x0 = SA_F64[10,10,100, 0,0,-10]
xf = @SVector zeros(n)
Q = Diagonal(@SVector fill(1.0, n))
R = Diagonal(@SVector fill(1e-1, m))
Qf = (N-1)*Q * 100

##
model = RobotZoo.DoubleIntegrator(3, gravity=SA[0,0,-9.81])
obj = LQRObjective(Q,R,Qf,xf,N)
cons = ConstraintList(n,m,N)
u_bnd = 100.0
connorm = NormConstraint(n, m, u_bnd, TO.SecondOrderCone(), :control)
add_constraint!(cons, GoalConstraint(xf), N)
add_constraint!(cons, connorm, 1:N-1)

prob = Problem(model, obj, xf, tf, x0=x0, constraints=cons)
solver = ALTROSolver(prob, projected_newton=false, show_summary=false)
solve!(solver)
@test iterations(solver) == 12

@test abs(maximum(norm.(controls(solver))) - u_bnd) < 1e-6
@test norm(states(solver)[end] - xf) < 1e-6

# prob = Problem(model, obj, xf, tf, x0=x0, constraints=cons)
# solver = ALTROSolver(prob, show_summary=true, verbose=2, projected_newton=false)
# @test benchmark_solve!(solver).allocs == 0
# benchmark_solve!(solver)