text:
__precompile__()
module SimpleFileServer
immutable File
    name::ASCIIString
    hash::ASCIIString
    size::Int64
end
abstract FileCommand

# TODO: is this really safe?
function safePath(f::Function,dirname::ASCIIString, path::ASCIIString)
    if contains(path, "..")
        throw(ErrorException("Path not allowed to contain '..'"))
    elseif contains(path, "~")
        throw(ErrorException("Path not allowed to contain '~'"))
    elseif startswith(path, "/") || endswith(path, "/")
        throw(ErrorException("Path not allowed to start or end with '/'"))

    else
        joinpath(dirname, path)|> f
    end
end

immutable Download <: FileCommand
    name::ASCIIString
    n_bytes_offset::Nullable{Tuple{Integer,Integer}}
end
immutable Delete <: FileCommand
    name::ASCIIString
end
immutable Upload <: FileCommand
    f::File
end

module Server
using SHA
using SimpleFileServer: Upload, Download, Delete, File, safePath, FileCommand
using Logging
@Logging.configure(level=DEBUG)
function HandleFileCommand(conn, args::Download, base::AbstractString)
    safePath(base, args.name) do path
        try
            if isfile(path)
                open(path, "r") do f
                    s =filesize(path)
                    n_bytes, offset = isnull(args.n_bytes_offset) ? (s, 0) : args.n_bytes_offset.value
                    serialize(conn, n_bytes)
                    @debug s, offset
                    seek(f, offset)
                    @debug "seeked"
                    buffer = Vector{UInt8}(Mmap.PAGESIZE)
                    while n_bytes != 0
                        b = readbytes!(f, buffer, min(Mmap.PAGESIZE, n_bytes))
                        write(conn, buffer[1:b])
                        @debug "remaining $n_bytes"
                        n_bytes -= b
                    end
                    @debug "sent"
                    flush(conn)
                end
            else
                serialize(conn, -1)
            end
        catch err
            #TODO Don't use this -- what's a better solution thought?
            serialize(conn, -1)
        finally
            flush(conn)
        end
    end
end
function HandleFileCommand(conn, args::Delete, base::AbstractString)
    try
        safePath(base, args.name) do path
            @debug path
            rm(path)
        end
        serialize(conn, Nullable())
    catch err
        serialize(conn, Nullable(err))
    end
end
function HandleFileCommand(conn, up::Upload, base::AbstractString)
    args = up.f
    s = args.size
    safePath(base, args.name) do path
        if length(args.hash) != 64
            serialize(conn,Nullable(ErrorException("Malformed Hash")))
        elseif isfile(path)
            serialize(conn,Nullable(ErrorException("File Already Exists")))
        end
        open(path, "w+") do f
            buffer = Vector{UInt8}(Mmap.PAGESIZE)
            while s != 0
                nbytes = readbytes!(conn, buffer, min(s, Mmap.PAGESIZE))
                write(f, buffer[1:nbytes])
                s -= nbytes
            end
            flush(f)
        end
        if SHA.sha256(open(path)) != args.hash
            rm(path)
            serialize(conn,Nullable(ErrorException("File did not Hash Properly")))
        end
        serialize(conn,Nullable{ErrorException}())
    end
end

function start(port::Int64, base::AbstractString, setup::Function)
    @info "Starting SimpleFileServer on port $port in path $base"
    base = base |>  expanduser |> realpath # Need to get rid of links
    server = listen(port)
    setup()
    while  true
        conn = accept(server)
        @debug "SimpleFileServer on $port has an incoming connection"
        @async begin
            try
                args = deserialize(conn)::FileCommand
                @debug args
                try
                    HandleFileCommand(conn, args, base)
                catch err
                    @error "Error While Handling: $err"
                end
            catch err
                @error "Error While Deserializing : $err"
            finally
                close(conn)
            end
        end
    end
end

function main()
    port = parse(Int64, ARGS[1])
    base = args[2]

    start(port, base)
end
end

module Client
using SHA
using SimpleFileServer: Upload, Download, Delete, File
using Logging
@Logging.configure(level=DEBUG)
export make, download, upload, delete
immutable t
    host::AbstractString
    port::Int64
end
@inline function make(host::AbstractString, port::Int64)
    t(host, port)
end
@inline function t_connect(c::t)
    connect(c.host, c.port)
end
function download(client::t, name::ASCIIString, to::IOStream, n_bytes_offset::Nullable{Tuple{Int64, Int64}})
    conn = t_connect(client)
    serialize(conn, Download(name,n_bytes_offset))
    flush(conn)
    s =deserialize(conn)::Int64
    if s == -1
        throw(Base.UVError("Remote Server $(client.host):$(client.port)", Base.UV_ENOENT))
    elseif s == 0
    else
        buffer = Vector{UInt8}( Mmap.PAGESIZE)
        while s != 0
            # TODO is this zero copy
            b = readbytes!(conn, buffer, min(s, Mmap.PAGESIZE) )
            write(to, buffer[1:b])
            s -= b
        end
        flush(to)
    end
    ()
end
function download(client::t, name::ASCIIString, to::ASCIIString, n_bytes_offset::Nullable{Tuple{Int64, Int64}}, append::Bool=false)
    open(to, append ? "a+":"w+") do f
        download(client, name, f, n_bytes_offset)
    end
    Mmap.mmap(to, Vector{UInt8})
end
download(client::t, name::ASCIIString, to::ASCIIString, n_bytes::Integer, offset::Integer, append::Bool=false )= download(client, name, to, Nullable((n_bytes,offset)), append)
download(client::t, name::ASCIIString, to::ASCIIString, append::Bool)= download(client, name, to, Nullable{Tuple{Int64, Int64}}(), append::Bool=false)
function upload(c::t, f::File, m::Vector{UInt8})
    conn = t_connect(c)
    serialize(conn, Upload(f))
    b =write(conn, m)
    flush(conn)
    @info "Uploading wrote $(b) bytes"
    m_err = deserialize(conn)
    if !isnull(m_err)
        throw(m_err.value)
    end
    close(conn)
end

function upload(c::t, name::ASCIIString,m::Vector{UInt8})
    f = File(name, SHA.sha256(m), length(m))
    upload(c, f, m)
    f
end
function upload(c::t, name::ASCIIString, path::ASCIIString)
    open(path, "r") do f
        m = Mmap.mmap(f, Vector{UInt8})
        upload(c, name, m)
    end
end
function delete(c::t, name::ASCIIString)
    conn = t_connect(c)
    serialize(conn, Delete(name))
    flush(conn)
    m_err = deserialize(conn)
    if !isnull(m_err)
        throw(m_err.value)
    end
end
end
end # module


text:
module TestConversions

using Test, DataFrames
using DataStructures: OrderedDict, SortedDict
const ≅ = isequal

@testset "Constructors to Base types" begin
    df = DataFrame()
    df[!, :A] = 1:5
    df[!, :B] = 1.0:5.0
    @test isa(Matrix(df), Matrix{Float64})
    @test isa(Matrix{Any}(df), Matrix{Any})
    @test isa(Matrix{Float64}(df), Matrix{Float64})
    @test isa(Array(df), Matrix{Float64})
    @test isa(Array{Any}(df), Matrix{Any})
    @test isa(Array{Float64}(df), Matrix{Float64})

    df = DataFrame()
    df.A = Vector{Union{Float64, Missing}}(1.0:5.0)
    df.B = Vector{Union{Float64, Missing}}(1.0:5.0)
    @test [df.A df.B] == Matrix(df)
    @test Matrix(df) isa Matrix{Union{Float64, Missing}}
    @test [df.A df.B] == Matrix{Any}(df)
    @test Matrix{Any}(df) isa Matrix{Any}
    @test [df.A df.B] == Matrix(df)
    @test Matrix{Int}(df) isa Matrix{Int}

    df[1, 1] = missing
    @test_throws ArgumentError Matrix{Float64}(df)
    na = Matrix{Union{Float64, Missing}}(df)
    naa = Matrix{Any}(df)
    nai = Matrix{Union{Int, Missing}}(df)
    @test isa(na, Matrix{Union{Float64, Missing}})
    @test isa(naa, Matrix{Union{Any, Missing}})
    @test isa(nai, Matrix{Union{Int, Missing}})
    @test na ≅ [df.A df.B]
    @test naa ≅ [df.A df.B]
    @test nai ≅ [df.A df.B]

    df = DataFrame()
    df.A = Vector{Union{Float64, Missing}}(1.0:5.0)
    df.B = Vector{Union{Float64, Missing}}(1.0:5.0)
    @test [df.A df.B] == Array(df)
    @test Array(df) isa Matrix{Union{Float64, Missing}}
    @test [df.A df.B] == Array{Any}(df)
    @test Array{Any}(df) isa Matrix{Any}
    @test [df.A df.B] == Array(df)
    @test Array{Int}(df) isa Matrix{Int}

    df[1, 1] = missing
    @test_throws ArgumentError Array{Float64}(df)
    na = Array{Union{Float64, Missing}}(df)
    naa = Array{Any}(df)
    nai = Array{Union{Int, Missing}}(df)
    @test isa(na, Matrix{Union{Float64, Missing}})
    @test isa(naa, Matrix{Union{Any, Missing}})
    @test isa(nai, Matrix{Union{Int, Missing}})
    @test na ≅ [df.A df.B]
    @test naa ≅ [df.A df.B]
    @test nai ≅ [df.A df.B]

    a = Union{Float64, Missing}[1.0, 2.0]
    b = Union{Float64, Missing}[-0.1, 3]
    c = Union{Float64, Missing}[-3.1, 7]
    di = Dict("a"=>a, "b"=>b, "c"=>c)

    df = DataFrame(di)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in sort(collect(keys(di)))]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
    df = DataFrame(od)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in keys(od)]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
    df = DataFrame(sd)
    @test isa(df, DataFrame)
    @test names(df) == [x for x in keys(sd)]
    @test df[!, :a] == a
    @test df[!, :b] == b
    @test df[!, :c] == c

    a = [1.0]
    di = Dict("a"=>a, "b"=>b, "c"=>c)
    @test_throws DimensionMismatch DataFrame(di)
end

@testset "conversion tests" begin
    df = DataFrame(a=1:3, b=4:6, c=7:9)
    dfv = @view df[2:3, 2:3]
    @test convert(DataFrame, dfv) == df[2:3, 2:3]
    @test convert(DataFrame, dfv) isa DataFrame
    @test DataFrame(dfv) == df[2:3, 2:3]
    @test DataFrame(dfv) isa DataFrame

    dfr = df[2, 2:3]
    @test convert(NamedTuple, dfr) == (b=5, c=8)
    @test convert(NamedTuple, dfr) isa NamedTuple
    @test NamedTuple(dfr) == (b=5, c=8)
    @test NamedTuple(dfr) isa NamedTuple

    gdf = groupby(df, [:b, :c])
    gk = keys(gdf)[2]
    @test convert(NamedTuple, gk) == (b=5, c=8)
    @test convert(NamedTuple, gk) isa NamedTuple
    @test NamedTuple(gk) == (b=5, c=8)
    @test NamedTuple(gk) isa NamedTuple
end

end # module

ext: jl


text:
# struct to model a cellset

@tm1def mutable struct Cellset
  ID::Union{String,Nothing}
  Edges::Union{Vector,Nothing}
  Elements::Union{Vector,Nothing}
  ElementAttributes::Union{Vector,Nothing}
  Subset::Union{Vector,Nothing}
  DefaultMember::Union{Vector,Nothing}
end

Cellset(id::AbstractString) = Cellset(Dict("ID" => id))

namefield(cellset::Cellset) = cellset.ID

# functions for endpoints

@api_default function get_all_cellsets(api::TM1API; options...)
  #params = Dict("\$expand" => "Cube,Axes,Cells")
  tm1_get_json(api, "Cellsets"; options...)
end

@api_default function get_cellset(api::TM1API, cellset_id::AbstractString; options...)
  params = Dict("\$expand" => "Cube,Axes,Cells")
  tm1_get_json(api, "Cellsets('" * cellset_id * "')"; params = params, options...)
end

@api_default function delete_cellset(api::TM1API, cellset_id::AbstractString; options...)
  tm1_delete(api, "Cellsets('" * cellset_id * "')"; options...)
end


text:
using Multirate
using DSP

function time_firfarrow{Th,Tx}( self::FIRFilter{FIRFarrow{Th}}, x::Vector{Tx} )
    kernel = self.kernel
    xLen   = length( x )
    println( "\nFIRFarrow speed test" )
    @printf( "\tresampling rate  %f\n", resampleRate )
                                                             # @printf( "\tpolynomial order %d\n", polyorder )
    @printf( "\tx type           %s\n", string(Tx) )
    @printf( "\tx length         %d\n", xLen )
    @printf( "\th type           %s\n", string(Th) )
    @printf( "\th length         %d\n", kernel.tapsPer𝜙*kernel.N𝜙 )
    @printf( "\tN𝜙               %d\n", kernel.N𝜙 )
    @printf( "\ttaps per 𝜙       %d\n", kernel.tapsPer𝜙 )
    ( y, elapsed, allocated, z ) = @timed filt( self, x )
    @printf( "\telapsed time (s) %1.3f\n", elapsed )
    @printf( "\tinput samples/s  %1.3e\n", xLen/elapsed )
    @printf( "\toutput samples/s %1.3e\n", length(y)/elapsed )
end

function time_firarbitrary{Th,Tx}( self::FIRFilter{FIRArbitrary{Th}}, x::Vector{Tx} )
    println( "\nFIRArbitrary Speed Test" )
    @printf( "\tresampling rate  %f\n", resampleRate )
    @printf( "\tx type           %s\n", string(Tx) )
    @printf( "\tx length         %d\n", xLen )
    @printf( "\th type           %s\n", string(Th) )
    @printf( "\th length         %d\n", hLen )
    @printf( "\tN𝜙               %d\n", N𝜙 )
    @printf( "\ttaps per 𝜙       %d\n", tapsPer𝜙 )
    ( y, elapsed, allocated, z ) = @timed filt( self, x )
    @printf( "\telapsed time (s) %1.3f\n", elapsed )
    @printf( "\tinput samples/s  %1.3e\n", xLen/elapsed )
    @printf( "\toutput samples/s %1.3e\n", length(y)/elapsed )
end

N𝜙           = 32                                            # Number of polyphase partitions
tapsPer𝜙     = 10                                            # N𝜙 * tapsPer𝜙 will be the length of out protoyTimepe filter taps
resampleRate = 1.0                                           # Can be any arbitrary resampling rate
polyorder    = 4                                             # Our taps will tranformed into
Th           = Float32
cutoffFreq   = min( 0.45/N𝜙, resampleRate/N𝜙 )               # N𝜙 is also the integer interpolation, so set cutoff frequency accordingly
hLen         = tapsPer𝜙*N𝜙                                   # Total number of filter taps
h            = firdes( hLen, cutoffFreq, DSP.kaiser ) .* N𝜙  # Generate filter taps and scale by polyphase interpolation (N𝜙)
xLen         = 10_000_000                                    # Number of signal samples

for resampleRate in ( 1.0, 1/2.123456789 ), Tx in ( Float32, Float64, Complex64, Complex128 )
    farrowfilt = FIRFilter( h, resampleRate, N𝜙, polyorder ) # Construct a FIRFilter{FIRFarrow} object
    arbfilt    = FIRFilter( h, resampleRate, N𝜙 )
    x          = rand( Tx, xLen )
    time_firfarrow( farrowfilt, x )
    time_firarbitrary( arbfilt, x )
end

max_stars_repo_path: examples/Arb-Farrow Speed Comparison.jl

max_issues_repo_path: examples/Arb-Farrow Speed Comparison.jl

max_forks_repo_path: examples/Arb-Farrow Speed Comparison.jl


text:
module WormlikeChain

include("random_utils.jl")

export BeadDefinition
export Chain
export SpecificForce
export ChainSystem
export append
export force_pe
include("interface.jl")

export refcpukernel!
include("cpukernels.jl")

include("bonded_forces_utils.jl")

end


text:
ChainRulesCore.@scalar_rule(airyai(x), airyaiprime(x))
ChainRulesCore.@scalar_rule(airyaiprime(x), x * airyai(x))
ChainRulesCore.@scalar_rule(airybi(x), airybiprime(x))
ChainRulesCore.@scalar_rule(airybiprime(x), x * airybi(x))
ChainRulesCore.@scalar_rule(besselj0(x), -besselj1(x))
ChainRulesCore.@scalar_rule(
    besselj1(x),
    (besselj0(x) - besselj(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(bessely0(x), -bessely1(x))
ChainRulesCore.@scalar_rule(
    bessely1(x),
    (bessely0(x) - bessely(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(dawson(x), 1 - (2 * x * Ω))
ChainRulesCore.@scalar_rule(digamma(x), trigamma(x))
ChainRulesCore.@scalar_rule(erf(x), (2 / sqrt(π)) * exp(-x * x))
ChainRulesCore.@scalar_rule(erfc(x), -(2 / sqrt(π)) * exp(-x * x))
ChainRulesCore.@scalar_rule(erfcinv(x), -(sqrt(π) / 2) * exp(Ω^2))
ChainRulesCore.@scalar_rule(erfcx(x), (2 * x * Ω) - (2 / sqrt(π)))
ChainRulesCore.@scalar_rule(erfi(x), (2 / sqrt(π)) * exp(x * x))
ChainRulesCore.@scalar_rule(erfinv(x), (sqrt(π) / 2) * exp(Ω^2))
ChainRulesCore.@scalar_rule(gamma(x), Ω * digamma(x))
ChainRulesCore.@scalar_rule(
    invdigamma(x),
    inv(trigamma(invdigamma(x))),
)
ChainRulesCore.@scalar_rule(trigamma(x), polygamma(2, x))

# binary
ChainRulesCore.@scalar_rule(
    besselj(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2
    ),
)
ChainRulesCore.@scalar_rule(
    besseli(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    bessely(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        (bessely(ν - 1, x) - bessely(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    besselk(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        -(besselk(ν - 1, x) + besselk(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh1(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        (hankelh1(ν - 1, x) - hankelh1(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh2(ν, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        (hankelh2(ν - 1, x) - hankelh2(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    polygamma(m, x),
    (
        ChainRulesCore.@thunk(error("not implemented")),
        polygamma(m + 1, x),
    ),
)
# todo: setup for common expr
ChainRulesCore.@scalar_rule(
    beta(a, b),
    (Ω*(digamma(a) - digamma(a + b)),
     Ω*(digamma(b) - digamma(a + b)),)
)
ChainRulesCore.@scalar_rule(
    logbeta(a, b),
    (digamma(a) - digamma(a + b),
     digamma(b) - digamma(a + b),)
)

# actually is the absolute value of the logorithm of gamma paired with sign gamma
ChainRulesCore.@scalar_rule(logabsgamma(x), digamma(x), ChainRulesCore.Zero())

ChainRulesCore.@scalar_rule(loggamma(x), digamma(x))


text:
"""
This module implements a memory, according to https://arxiv.org/abs/1703.03129,
that can be used as the last layer in a NN in the Flux framework.
"""

export KNNmemory, query, trainQuery!

"""
    KNNmemory{T}
Structure that contains all the memory data.
"""
mutable struct KNNmemory{T <: Real}
    M::Matrix{T} # keys in the memory
    V::Vector{<:Integer} # values in the memory (labels)
    A::Vector{<:Integer} # age of a given key-value pair
    k::Integer # number of neighbors used in kNN
    α::Real # parameter setting the required distance between the nearest positive and negative sample in the kNN
    encoder # Transformation function that allows to store item in memory in one space and compute the distance in another
			# The final space has to be a hypersphere!
	encodeM # transposed transformation so it can be used on `M`

    """
        KNNmemory{T}(memorySize, keySize, k, labelCount, [α])
    Memory constructor that initializes it with random keys and random labels.
    # Arguments
    - `memorySize::Integer`: number of keys that can be stored in the memoryUpdate!
    - `keySize::Integer`: length of a key
    - `k::Integer`: number of k nearest neighbors to look for
    - `labelCount::Integer`: number of labels that are in the dataset (used for the random initialization)
	- `encoder`: Transformation function that allows to store item in memory in one space and compute the distance in another. Default is identity.
    - `α::Real`: parameter of the memory loss function that determines required distance between clusters
    """
    function KNNmemory{T}(memorySize::Integer, keySize::Integer, k::Integer, labelCount::Integer, encoder = x -> normalizecolumns(x), α::Real = 0.1) where T
        M = rand(T, memorySize, keySize) .* 2 .- 1
        V = rand(0:(labelCount - 1), memorySize)
        A = zeros(Int, memorySize)

        for i = 1:memorySize
            M[i,:] = normalize(M[i,:])
        end

        new(M, V, A, k > memorySize ? memorySize : k, convert(T, α), encoder, m -> collect(encoder(m')'))
    end
end

"""
    findNearestPositiveAndNegative(memory, kLargestIDs, v)
For given set of k nearest neighbours, find the closest two that have the same label `v` and a different label respectively.
`kLargestIDs::Vector{<:Integer}` contains k indices leading to the k most similar keys in the memorySize.
"""
function findNearestPositiveAndNegative(memory::KNNmemory, kLargestIDs::Vector{<:Integer}, v::Integer)
    nearestPositiveID = nothing
    nearestNegativeID = nothing

    # typically this should not result into too many iterations
    for i in 1:memory.k
        if nearestPositiveID == nothing && memory.V[kLargestIDs[i]] == v
            nearestPositiveID = kLargestIDs[i]
        end
        if nearestNegativeID == nothing && memory.V[kLargestIDs[i]] != v
            nearestNegativeID = kLargestIDs[i]
        end
        if nearestPositiveID != nothing && nearestNegativeID != nothing
            break
        end
    end

    #= We assume that there exists such i that memory.V[i] == v
        and also such j that memory.V[j] != v
list
        We also assume that this won't happen very often, otherwise,
        we would need to randomize this selection (possible TODO) =#

    if nearestPositiveID == nothing
        nearestPositiveID = Base.argmax(memory.V .== v)
    end
    if nearestNegativeID == nothing
        nearestNegativeID = Base.argmax(memory.V .!= v)
    end

    return nearestPositiveID, nearestNegativeID
end

"""
    memoryLoss(memory, q, nearestPosAndNegIDs)
Loss generated by the memory based on the lookup of a key-value pair for the key `q` - exactly as in the paper.
`nearestPosAndNegIDs::Tuple` represents ids of the closest key with the same and a different label respectively.
"""
memoryLoss(memory::KNNmemory{T}, q::AbstractArray, nearestPosAndNegIDs::Tuple) where {T} = memoryLoss(memory, q, nearestPosAndNegIDs...)

function memoryLoss(memory::KNNmemory{T}, normalizedQuery::AbstractArray, nearestPositiveID::Integer, nearestNegativeID::Integer) where {T}
    loss = max(dot(normalizedQuery, memory.encoder(memory.M[nearestNegativeID, :])) - dot(normalizedQuery, memory.encoder(memory.M[nearestPositiveID, :])) + memory.α, 0)
end

normalizecolumns(m) = m ./ sqrt.(sum(m .^ 2, dims = 1) .+ eps(eltype(Flux.Tracker.data(m))))

function findInverse(func, invsize, target, precision, learningrate, maxiter)
	loss(x) = sum((func(x) .- target) .^ 2)
	i = 1
	newx = param(rand(Float64, invsize))
	last = newx .+ 1
	opt = ADAM(learningrate)
	while (loss(newx) > precision) & (i < maxiter)
		last = newx
		l = loss(newx)
		Flux.Tracker.back!(l)
		Δ = Flux.Optimise.apply!(opt, newx.data, newx.grad)
        newx.data .-= Δ
		newx.grad .= 0
		i += 1
	end
	return newx
end

"""
    memoryUpdate!(memory, q, v, nearestNeighbourID)
It computes the appropriate update of the memory after a key-value pair was lookedup in it for the key `q` and expected label `v`.
"""
function memoryUpdate!(memory::KNNmemory{T}, q::Vector, v::Integer, nearestNeighbourID::Integer) where {T}
    # If the memory return the correct value for the given key, update the centroid
    if memory.V[nearestNeighbourID] != 1 && memory.V[nearestNeighbourID] == v # TODO: This is a hack to not move anomalies - should be done in a better way!

		target = collect(normalize(Flux.Tracker.data(memory.encoder(q) + memory.encoder(memory.M[nearestNeighbourID, :])))')
		memory.M[nearestNeighbourID, :] = collect(findInverse(memory.encoder, size(memory.M, 2), target', 0.001, 0.01, 1000)')
        memory.A[nearestNeighbourID] = 0

    # If the memory did not return the correct value for the given key, store the key-value pair instead of the oldest element
    else
        oldestElementID = Base.argmax(memory.A .+ rand(1:5))
        memory.M[oldestElementID, :] = Flux.Tracker.data(q)
        memory.V[oldestElementID] = v
        memory.A[oldestElementID] = 0
    end
end

"""
    increaseMemoryAge(memory)
Update the age of all items
"""
function increaseMemoryAge(memory::KNNmemory)
    memory.A .+= 1;
end

"""
    query(memory, q)
Returns the nearest neighbour's value and its confidence level for a given key `q` but does not modify the memory itself.
"""
function query(memory::KNNmemory{T}, q::AbstractArray{T, N} where N) where {T}
    similarity = memory.encodeM(memory.M) * Flux.Tracker.data(memory.encoder(q))
    values = memory.V[Flux.onecold(similarity)]

    function probScorePerQ(index)
        kLargestIDs = collect(partialsortperm(similarity[:, index], 1:memory.k, rev = true))
        probsOfNearestKeys = softmax(similarity[kLargestIDs, index])
        nearestValues = memory.V[kLargestIDs]
        return sum(probsOfNearestKeys[nearestValues .== 1])
    end

    return values, map(probScorePerQ, 1:size(q, 2)) # basicaly returns the nearest value in the memory + sum of probs of anomalies that are in the k-nearest
end

# logc(p, κ) = (p / 2 - 1) * log(κ) - (p / 2) * log(2π) - log(besseli(p / 2 - 1, κ))
logc(p, κ) = (p ./ 2 .- 1) .* log.(κ) .- (p ./ 2) .* log(2π) .- κ .- log.(besselix(p / 2 - 1, κ))

# log likelihood of one sample under the VMF dist with given parameters
# log_vmf(x, μ, κ) = κ * μ' * x .+ log.(c(length(μ), κ))
log_vmf(x, μ, κ) = κ * μ' * x .+ logc(length(μ), κ)

vmf_mix_lkh(x, μs, κ) = vmf_mix_lkh(x, μs, κ, size(μs, 2))
function vmf_mix_lkh(x, μs, κ, μlength::Integer)
    κs = ones(μlength) .* κ # This is quite arbitrary as we don't really know what to use for kappa but it shouldn't matter if it is the same
	l = 0
	for K in 1:μlength
		l += exp.(log_vmf(x, μs[:, K], κs[K]))
	end
	l /= μlength
	return l
end

"""
    prob_query(memory, q)
Returns the nearest neighbour's value and its confidence level for a given key `q` but does not modify the memory itself.
"""
function prob_query(memory::KNNmemory{T}, q::AbstractArray, κ) where {T}
    normq = Flux.Tracker.data(memory.encoder(q))
	encmem = memory.encodeM(memory.M)
    similarity = encmem * normq
    values = memory.V[Flux.onecold(similarity)]

    function probScorePerQ(index)
        kLargestIDs = collect(partialsortperm(similarity[:, index], 1:memory.k, rev = true))
		nearestAnoms = collect(encmem[kLargestIDs, :][memory.V[kLargestIDs] .== 1, :]')
		if length(nearestAnoms) == 0
			return 0
		elseif size(nearestAnoms, 2) == memory.k
			return 1
		end
		nearestNormal = collect(encmem[kLargestIDs, :][memory.V[kLargestIDs] .== 0, :]')
		pxgivena = vmf_mix_lkh(normq[:, index], nearestAnoms, κ)
		pxgivenn = vmf_mix_lkh(normq[:, index], nearestNormal, κ)
		return pxgivena / (pxgivena + pxgivenn)
    end

    return values, map(probScorePerQ, 1:size(q, 2)) # basicaly returns the nearest value in the memory + sum of probs of anomalies that are in the k-nearest
end

"""
    trainQuery!(memory, q, v)
Query 'q' to the memory that does update its content and returns a loss for expected outcome label `v`.
"""
trainQuery!(memory::KNNmemory{T}, q::AbstractArray{T, N} where N, v::Integer) where {T} = trainQuery!(memory, q, [v])

function trainQuery!(memory::KNNmemory{T}, q::AbstractArray{Tt, N} where {Tt, N}, v::Vector{<:Integer}) where {T}
    # Find k nearest neighbours and compute losses
    batchSize = size(q, 2)
    normalizedQuery = Flux.Tracker.data(memory.encoder(q))
    similarity = memory.encodeM(memory.M) * normalizedQuery # computes all similarities of all qs and all keys in the memory at once
    loss::Flux.Tracker.TrackedReal{T} = 0 # loss must be tracked; otherwise flux cannot use it
    nearestNeighbourIDs = zeros(Integer, batchSize)

    for i in 1:batchSize
        kLargestIDs = collect(partialsortperm(similarity[:, i], 1:memory.k, rev = true))
        nearestNeighbourIDs[i] = kLargestIDs[1];
        loss += memoryLoss(memory, normalizedQuery[:, i], findNearestPositiveAndNegative(memory, kLargestIDs, v[i]))
    end

    # Memory update - cannot be done above because we have to compute all losses before changing the memory
    for i in 1:batchSize
		if i % 5 == 1
			println("$i/$batchSize")
		end
        memoryUpdate!(memory, q[:, i], v[i], nearestNeighbourIDs[i])
    end
    increaseMemoryAge(memory)

    return loss / batchSize
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
@testset "transform" begin
    rng = MersenneTwister(123546)
    x = randn(rng, 8)
    XV = [randn(rng, 5) for _ in 1:6]
    XC = ColVecs(randn(rng, 5, 10))
    XR = RowVecs(randn(rng, 11, 3))
    @testset "IdentityTransform($(typeof(x)))" for x in [x, XV, XC, XR]
        @test IdentityTransform()(x) == x
        @test map(IdentityTransform(), x) == x
    end
    test_ADs(() -> transform(SEKernel(), IdentityTransform()))
end


text:
using Random

include("Equation.jl")

macro def(name, definition)
  return quote
      macro $(esc(name))()
          esc($(Expr(:quote, definition)))
      end
  end
end

@def addNodeFields begin
    equation::Union{Nothing,NodeEquation}
end

abstract type Tree end

abstract type Terminal <: Tree end

mutable struct Variable <: Terminal
    @addNodeFields
    variableNumber::UInt
    semantic::Array{<:Real,1}
    Variable(num::UInt; semantic=[]) = new(nothing, num, semantic)
end

mutable struct Constant <: Terminal
    @addNodeFields
    semantic::Real
    Constant(value::Real) = new(nothing, value)
end

mutable struct RandomConstant <: Terminal
    @addNodeFields
    lowerLimit::Real
    upperLimit::Real
    RandomConstant(lowerLimit::Real, upperLimit::Real) = new(nothing, lowerLimit, upperLimit)
end

@def addNonTerminalFields begin
    name::String
    semantic::Union{Nothing,Semantic}
    evalFunction::Function
end

abstract type NonTerminal <: Tree end

mutable struct BinaryNode <: NonTerminal
    @addNodeFields
    @addNonTerminalFields
    child1::Tree
    child2::Tree
    functionEquationChild1::Function
    functionEquationChild2::Function
    function BinaryNode(name::String, evalFunction::Function, child1::Tree, child2::Tree, functionEquationChild1::Function, functionEquationChild2::Function; semantic=nothing)
        return new(nothing, name, semantic, evalFunction, child1, child2, functionEquationChild1, functionEquationChild2);
    end;
end

mutable struct NonBinaryNode <: NonTerminal
    @addNodeFields
    @addNonTerminalFields
    children::Array{Tree,1}
    functionEquationChildren::Array{Function,1}
    function NonBinaryNode(name::String, evalFunction::Function, functionEquationChildren::Array{Function,1}; semantic=nothing, children::Array{Tree,1}=[])
        return new(nothing, name, semantic, evalFunction, children, functionEquationChildren);
    end;
    function NonBinaryNode(type::Int, name::String, evalFunction::Function, functionEquationChildren::Array{Function,1}, children::Array{Tree,1}; semantic=nothing)
        return new(nothing, name, semantic, evalFunction, children, functionEquationChildren);
    end;
end

AddNode(child1::Tree, child2::Tree) = BinaryNode("+", addSemantics, child1, child2, equationAddChild1, equationAddChild2);
SubNode(child1::Tree, child2::Tree) = BinaryNode("-", subSemantics, child1, child2, equationSubChild1, equationSubChild2);
MulNode(child1::Tree, child2::Tree) = BinaryNode("*", mulSemantics, child1, child2, equationMulChild1, equationMulChild2);
DivNode(child1::Tree, child2::Tree) = BinaryNode("/", divSemantics, child1, child2, equationDivChild1, equationDivChild2);

clearEvaluationValues!(node::Constant) = nothing;
clearEvaluationValues!(node::RandomConstant) = error("");
clearEvaluationValues!(node::Variable) = nothing;
function clearEvaluationValues!(node::BinaryNode)
    node.semantic = nothing;
    clearEvaluationValues!(node.child1);
    clearEvaluationValues!(node.child2);
    return nothing;
end
function clearEvaluationValues!(node::NonBinaryNode)
    node.semantic = nothing;
    for childNode in node.children
        clearEvaluationValues!(childNode);
    end;
    return nothing;
end

# These 3 functions clear the equations except in the given path
clearEquations!(node::Terminal) = ( node.equation = nothing; )
function clearEquations!(node::BinaryNode; exceptInPath=nothing)
    if (exceptInPath==nothing)
        node.equation = nothing;
        clearEquations!(node.child1);
        clearEquations!(node.child2);
    else
        if (isempty(exceptInPath))
            clearEquations!(node.child1);
            clearEquations!(node.child2);
        elseif (exceptInPath[1]==1)
            clearEquations!(node.child1, exceptInPath[2:end]);
            clearEquations!(node.child2);
        else
            clearEquations!(node.child1);
            clearEquations!(node.child2, exceptInPath[2:end]);
        end;
    end;
end
function clearEquations!(node::NonBinaryNode; exceptInPath=nothing)
    if (exceptInPath==nothing)
        node.equation = nothing;
        for child in node.children
            clearEquations!(child);
        end;
    else
        if (isempty(exceptInPath))
            for child in node.children
                clearEquations!(child);
            end;
        else
            nextStep = exceptInPath[1];
            for numChild in 1:length(node.children)
                if (numChild==nextStep)
                    clearEvaluationValues!(node.children[numChild], exceptInPath[2:end]);
                else
                    clearEvaluationValues!(node.children[numChild]);
                end;
            end;
        end;
    end;
end

function clearEquations!(node::BinaryNode; exceptInPaths=nothing)
    if (exceptInPaths==nothing)
        node.equation = nothing;
        clearEquations!(node.child1);
        clearEquations!(node.child2);
    else
        if (isempty(exceptInPaths))
            clearEquations!(node.child1);
            clearEquations!(node.child2);
        else
            firstSteps = [path[1] for path in exceptInPaths];
            has1 = anyequal(1,firstSteps);
            has2 = anyequal(2,firstSteps);
            if (has1) && (has2)
                followingSteps = exceptInPaths[firstSteps.==1];
                followingSteps = [path[2:end] for path in followingSteps];
                followingSteps = followingSteps[!isempty.(followingSteps)];
                clearEquations!(node.child1, exceptInPaths=followingSteps);
                followingSteps = exceptInPaths[firstSteps.==2];
                followingSteps = [path[2:end] for path in followingSteps];
                followingSteps = followingSteps[!isempty.(followingSteps)];
                clearEquations!(node.child1, exceptInPaths=followingSteps);
            elseif (has1)
                followingSteps = [path[2:end] for path in exceptInPaths];
                followingSteps = followingSteps[!isempty.(followingSteps)];
                clearEquations!(node.child1, exceptInPaths=followingSteps);
                clearEquations!(node.child2);
            elseif (has2)
                followingSteps = [path[2:end] for path in exceptInPaths];
                followingSteps = followingSteps[!isempty.(followingSteps)];
                clearEquations!(node.child1);
                clearEquations!(node.child2, exceptInPaths=followingSteps);
            else
                error("Don't know which path to go")
            end;
        end;
    end;
end

import Base.string
string(node::Constant) = node.semantic>=0 ? string(node.semantic) : string("(",node.semantic,")");
string(node::RandomConstant) = error("");
string(node::Variable) = string("X", node.variableNumber);
string(node::BinaryNode) = string("(",string(node.child1),node.name,string(node.child2),")")
function string(node::NonBinaryNode)
    if (length(node.children)==2)
        text = string("(",string(node.children[1]),node.name,string(node.children[2]),")")
    else
        text = string(node.name, "(");
        for numChildren = 1:length(node.children)
            text = string(text, string(node.children[numChildren]));
            if (numChildren!=length(node.children))
                text = string(text, ",");
            end;
        end;
        text = string(text, ")");
    end;
    return text;
end;

vectorString(node::Constant; dataInRows=true) = node.semantic>=0 ? string(node.semantic) : string("(",node.semantic,")");
vectorString(node::RandomConstant; dataInRows=true) = error("");
vectorString(node::Variable; dataInRows=true) = dataInRows ? string("X[:,", node.variableNumber,"]") : string("X[", node.variableNumber,",:]");
vectorString(node::BinaryNode; dataInRows=true) = string("(", vectorString(node.child1; dataInRows=dataInRows), " .", node.name, " ", vectorString(node.child2; dataInRows=dataInRows), ")")
function vectorString(node::NonBinaryNode; dataInRows=true)
    if (length(node.children)==2)
        text = string("(", vectorString(node.children[1]; dataInRows=dataInRows), " .", node.name, " ", vectorString(node.children[2]; dataInRows=dataInRows), ")")
    else
        text = string(node.name, ".(");
        for numChildren = 1:length(node.children)
            text = string(text, vectorString(node.children[numChildren]; dataInRows=dataInRows));
            if (numChildren!=length(node.children))
                text = string(text, ",");
            end;
        end;
        text = string(text, ")");
    end;
    return text;
end;

latexString(node::Constant) = string(node);
latexString(node::RandomConstant) = string(node);
latexString(node::Variable) = string("X_{", node.variableNumber,"}");
latexString(node::BinaryNode) = (node.name=="/") ? string("\\frac{",latexString(node.child1),"}{",latexString(node.child2),"}") : string("\\left(",latexString(node.child1),(node.name=="*") ? "\\cdot " : node.name,latexString(node.child2),"\\right)")
latexString(node::NonBinaryNode) = string(node);
function latexString(node::NonBinaryNode)
    if (length(node.children)==2)
        text = string("\\left(",string(node.children[1]),node.name,string(node.children[2]),"\\right)")
    else
        text = string(node.name, "\\left(");
        for numChildren = 1:length(node.children)
            text = string(text, string(node.children[numChildren]));
            if (numChildren!=length(node.children))
                text = string(text, ',');
            end;
        end;
        text = string(text, "\\right)");
    end;
    return text;
end;

evaluateTree(tree::Terminal; checkForErrors=false) = tree.semantic;
evaluateTree(tree::RandomConstant; checkForErrors=false) = error("");
function evaluateTree(tree::BinaryNode; checkForErrors=false)
    if (tree.semantic==nothing)
        tree.semantic = tree.evalFunction( evaluateTree(tree.child1;checkForErrors=checkForErrors), evaluateTree(tree.child2;checkForErrors=checkForErrors), nothing );
    end;
    return tree.semantic;
end;
function evaluateTree(tree::NonBinaryNode; checkForErrors=false)
    if (tree.semantic==nothing)
        evaluationChildren = [evaluateTree(child; checkForErrors=checkForErrors) for child in tree.children];
        tree.semantic = tree.evalFunction( evaluationChildren... );
    end;
    if (checkForErrors)
        @assert(!any(isinf.(tree.semantic)));
        @assert(!any(isnan.(tree.semantic)));
    end;
    return tree.semantic;
end;

reevaluatePath(tree::Terminal, path::Array{Int64,1}, indexPath::Int64; checkForErrors=false) = tree.semantic;
reevaluatePath(tree::RandomConstant, path::Array{Int64,1}, indexPath::Int64; checkForErrors=false) = error("");
function reevaluatePath(tree::BinaryNode, path::Array{Int64,1}, indexPath::Int64=1; checkForErrors=false)
    if (indexPath>length(path))
        semanticChild1 =   evaluateTree(tree.child1;            checkForErrors=checkForErrors);
        semanticChild2 =   evaluateTree(tree.child2;            checkForErrors=checkForErrors);
    elseif (path[indexPath]==1)
        semanticChild1 = reevaluatePath(tree.child1, path, indexPath+1; checkForErrors=checkForErrors);
        semanticChild2 =   evaluateTree(tree.child2;            checkForErrors=checkForErrors);
    elseif (path[indexPath]==2)
        semanticChild1 =   evaluateTree(tree.child1;            checkForErrors=checkForErrors);
        semanticChild2 = reevaluatePath(tree.child2, path, indexPath+1; checkForErrors=checkForErrors);
    else
        error("Don't know which path to go")
    end;
    tree.semantic = tree.evalFunction( semanticChild1, semanticChild2, tree.semantic );
    if (checkForErrors)
        @assert(!any(isinf.(tree.semantic)));
        @assert(!any(isnan.(tree.semantic)));
    end;
    return tree.semantic;
end;

numNodes(node::Terminal) = UInt(1);
numNodes(node::BinaryNode) = UInt(1 + numNodes(node.child1) + numNodes(node.child2));
numNodes(node::NonBinaryNode) = UInt(1 + sum([numNodes(child) for child in node.children]));

height(node::Terminal) = UInt(1);
height(node::BinaryNode) = UInt(1) + max(height(node.child1), height(node.child2));
height(node::NonBinaryNode) = UInt(1) + maximum([height(child) for child in node.children]);

clone(node::Constant) = Constant(node.semantic);
clone(node::RandomConstant) = Constant(rand()*(node.upperLimit-node.lowerLimit) + node.lowerLimit, node.type);
clone(node::Variable) = Variable(node.variableNumber; semantic=((node.semantic==nothing) ? node.semantic : copy(node.semantic)));
clone(node::BinaryNode) = BinaryNode(node.name, node.evalFunction, clone(node.child1), clone(node.child2), node.functionEquationChild1, node.functionEquationChild2; semantic=((node.semantic==nothing) ? node.semantic : copy(node.semantic)));
clone(node::NonBinaryNode) = NonBinaryNode(node.name, node.evalFunction, node.equationsFunction; semantic=((node.semantic==nothing) ? node.semantic : copy(node.semantic)), children=convert(Array{Tree,1},[clone(child) for child in node.children]));

iterateTree(tree::Terminal) = (convert(Array{Tree,1},[tree]), [0x0000000000000001], [0x0000000000000001], [0x0000000000000000], [0x0000000000000001]);
function iterateTree(tree::BinaryNode)
    (nodesChild1, heightsChild1, depthsChild1, numChildrenChild1, numNodesChild1) = iterateTree(tree.child1);
    (nodesChild2, heightsChild2, depthsChild2, numChildrenChild2, numNodesChild2) = iterateTree(tree.child2);
    nodes = [nodesChild1; nodesChild2]; pushfirst!(nodes, tree);
    heights = [heightsChild1; heightsChild2]; pushfirst!(heights, 0x0000000000000000);
    depths = [depthsChild1; depthsChild2]; pushfirst!(depths, 0x0000000000000001);
    numChildren = [numChildrenChild1; numChildrenChild2]; pushfirst!(numChildren, 0x0000000000000001);
    numNodes = [numNodesChild1; numNodesChild2]; pushfirst!(numNodes, 0x0000000000000000);
    heights[1] = maximum(heights[2:end])+1;
    numNodes[1] = length(numNodes);
    depths[2:end] .+= 1;
    return (nodes, heights, depths, numChildren, numNodes);
end
function iterateTree(tree::NonBinaryNode)
    nodes = [tree];
    heights = [0x0000000000000000];
    depths = [0x0000000000000001];
    numChildren = [length(tree.children)];
    numNodes = [0x0000000000000000];
    for child in tree.children
        (nodesThisChild, tiposThisChild, heightsThisChild, depthsThisChild, numChildrenThisChild, numNodesThisChild) = iterateTree(child);
        nodes = [nodes; nodesThisChild];
        # hcat(tipos, tiposEsteHijo);
        heights = [heights; heightsThisChild];
        depths = [depths; depthsThisChild.+1];
        numChildren = [numChildren; numChildrenThisChild];
        numNodes = [numNodes; numNodesThisChild];
    end;
    heights[1] = maximum(heights[2:end])+1;
    numNodes[1] = length(numNodes);
    depths[2:end] .+= 1;
    return (nodes, heights, depths, numChildren, numNodes);
end

calculateEquations!(tree::Terminal, equation::NodeEquation; path=nothing, indexPath=1, checkForErrors=false) = ( tree.equation = equation; );
function calculateEquations!(tree::BinaryNode, equation::NodeEquation; path=nothing, indexPath=1, checkForErrors=false)
    if (checkForErrors)
        @assert (tree.semantic != nothing);
        checkEquation(equation);
    end;
    tree.equation = equation;

    if (path==nothing)
        calculateEquations!(tree.child1, tree.functionEquationChild1( evaluateTree(tree.child2; checkForErrors=checkForErrors), equation); checkForErrors=checkForErrors );
        calculateEquations!(tree.child2, tree.functionEquationChild2( evaluateTree(tree.child1; checkForErrors=checkForErrors), equation); checkForErrors=checkForErrors );
    else
        (indexPath>length(path)) && return
        if (path[indexPath]==1)
            calculateEquations!(tree.child1, tree.functionEquationChild1( evaluateTree(tree.child2; checkForErrors=checkForErrors), equation); path=path, indexPath=indexPath+1, checkForErrors=checkForErrors );
        elseif (path[indexPath]==2)
            calculateEquations!(tree.child2, tree.functionEquationChild2( evaluateTree(tree.child1; checkForErrors=checkForErrors), equation); path=path, indexPath=indexPath+1, checkForErrors=checkForErrors );
        else
            error("Don't know which path to go");
        end;
    end;
end;

setVariableValues!(tree::Terminal, variableValues::Array{<:Array{<:Real,1}}) = ()
setVariableValues!(tree::Variable, variableValues::Array{<:Array{<:Real,1}}) = (tree.semantic = variableValues[tree.variableNumber]; )
setVariableValues!(tree::BinaryNode, variableValues::Array{<:Array{<:Real,1}}) = (setVariableValues!(tree.child1, variableValues); setVariableValues!(tree.child2, variableValues); tree.semantic = nothing; )
setVariableValues!(tree::NonBinaryNode, variableValues::Array{<:Array{<:Real,1}}) = (for child in tree.children setVariableValues!(child, variableValues); end; tree.semantic = nothing; )

findNode(tree::Terminal, node::Tree) = (tree==node) ? Array{Int,1}([]) : nothing
function findNode(tree::BinaryNode, node::Tree)
    (tree==node) && return Array{Int,1}([]);
    path = findNode(tree.child1, node);
    if (path!=nothing)
        pushfirst!(path, 1);
        return path;
    end;
    path = findNode(tree.child2, node);
    if (path!=nothing)
        pushfirst!(path, 2);
        return path;
    end;
    return nothing;
end;
function findNode(tree::NonBinaryNode, node::Tree)
    (tree==node) && return Array{Int,1}([]);
	for numChild in 1:length(tree.children)
        path = findNode(tree.children[numChild], node);
        if (path!=nothing)
            pushfirst!(path, numChild);
            return path;
        end;
    end;
    return nothing;
end;

replaceSubtree!(tree::Terminal, oldSubTree::Tree, newSubTree::Tree) = false;
function replaceSubtree!(tree::BinaryNode, oldSubTree::Tree, newSubTree::Tree)
    replaced = false;
    if (tree.child1==oldSubTree)
        tree.child1 = newSubTree;
        replaced = true;
    elseif (tree.child2==oldSubTree)
        tree.child2 = newSubTree;
        replaced = true;
    else
        replaced = replaceSubtree!(tree.child1, oldSubTree, newSubTree);
        if !replaced
            replaced = replaceSubtree!(tree.child2, oldSubTree, newSubTree);
        end;
    end;
    if (replaced)
        tree.semantic = reevaluatePath(tree, Int64[]);
        return true;
    end;
    return false;
end
function replaceSubtree!(tree::NonBinaryNode, oldSubTree::Tree, newSubTree::Tree)
	for numChild in 1:length(tree.children)
        child = tree.children[numChild]
        if (child==oldSubTree)
            tree.children[numChild] = newSubTree;
            tree.semantic = reevaluatePath(tree, Int64[]);
            return true;
        else
            if (replaceSubtree!(child, oldSubTree, newSubTree))
                tree.semantic = reevaluatePath(tree, Int64[]);
                return true;
            end;
        end;
    end;
    return false;
end

function checkEquations(tree::Tree, mse::Real)
    @assert(!any(isnan.(tree.semantic)));
    @assert(!any(isinf.(tree.semantic)));
    if (tree.equation!=nothing)
        checkEquation(tree.equation);
        @assert(equal(mse, calculateMSEFromEquation(tree.semantic, tree.equation; checkForErrors=true)));
        @assert(equal(mse, calculateMSEFromEquationEfficient(tree.semantic, tree.equation; checkForErrors=true)));
        @assert(equal(mse, calculateMSEFromEquationMoreEfficient(tree.semantic, tree.equation; checkForErrors=true)));
    end;
    if (isa(tree,BinaryNode))
        checkEquations(tree.child1, mse)
        checkEquations(tree.child2, mse)
    elseif (isa(tree,NonBinaryNode))
        for numChild in 1:length(tree.children)
            checkEquations(tree.children[numChild], mse)
        end;
    end;
end;

function buildInitialTree(order::Int64, variables::Array{Variable,1}, constant::Real, maximumNodes::Number)

    function buildTreeOrder(order::Int64, variables::AbstractArray{Variable,1})

        function buildVariableCombinationTreeList(thisOrder::Int64, variablesThisOrder::AbstractArray{Variable,1})
            (thisOrder==0) && return [Constant(eltype(constant)(0.))];
            treeListThisOrder = [];
            for numVariable in 1:length(variablesThisOrder)
                subTrees = buildVariableCombinationTreeList(thisOrder-1, view(variablesThisOrder, numVariable:length(variablesThisOrder)));
                for subTree in subTrees
                    push!(treeListThisOrder, MulNode(subTree, clone(variablesThisOrder[numVariable])) );
                end;
            end;
            return treeListThisOrder;
        end;

        treeList = buildVariableCombinationTreeList(order, variables);
        tree = treeList[1];
        for subTree in view(treeList, 2:length(treeList))
            tree = AddNode( tree, subTree );
        end;
        return tree;
    end;

    outputTree = Constant(constant);

    numNodesCurrent = 1;

    for i in 1:order
        treeOrder = buildTreeOrder(i, variables);
        numNodesTreeOrder = numNodes(treeOrder);
        if (numNodesCurrent+numNodesTreeOrder+1<=maximumNodes)
            outputTree = AddNode(outputTree, treeOrder);
            numNodesCurrent += numNodesTreeOrder + 1;
        else
            # The rest of the trees of higher order are going to have even more nodes,
            #  so exit this function
            return outputTree;
        end;
    end;

    return outputTree;
end;


text:
function C2(t, fft_u, fft_v,
	    A1, A2, matr, conjmatr,
	    sigma, llambda, Ktaubis,
	    epsilon, Ntaumm)

    h1u, h1v = ftau(t, fft_u, fft_v, A1, A2, matr, conjmatr, sigma, llambda)

    h1u_fft = fft(h1u, dims=0)
    h1v_fft = fft(h1v, dims=0)

    dtu0u = h1u_fft[1, :] / Ntaumm
    dtu0v = h1v_fft[1, :] / Ntaumm

    h1u_fft[1, :] .= 0.0
    h1v_fft[1, :] .= 0.0

    h1u = epsilon * ifft(h1u_fft ./ (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t, fft_u, fft_v,
                              dtu0u, dtu0v,
			      A1, A2, matr,
			      conjmatr, sigma,
			      llambda, Ntaumm)

    champu2, champv2 = dtftau(t, fft_u, fft_v,
			      A1, A2, matr,
			      conjmatr, sigma,
			      llambda, Ntaumm)

    champu .= champu1 + champu2
    champv .= champv1 + champv2
    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)
    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0
    dth1u = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth1v = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    h2u, h2v = ftau(t,
		    fft_u .+ h1u,
		    fft_v .+ h1v,
		    A1, A2, matr, conjmatr,
                    sigma, llambda)
    h2u = h2u .- dth1u
    h2v = h2v .- dth1v
    h2u_fft = fft(h2u, dims=1)
    h2v_fft = fft(h2v, dims=1)
    h2u_fft[1, :] .= 0.0
    h2v_fft[1, :] .= 0.0
    h2u = epsilon * ifft(h2u_fft ./ (1im * Ktaubis), dims=1)
    h2v = epsilon * ifft(h2v_fft ./ (1im * Ktaubis), dims=1)

    C2u = fft_u .+ h2u
    C2v = fft_v .+ h2v

    return C2u, C2v

end


text:
using NonlinearSchrodinger
using RecipesBase
using Test

@testset "NonlinearSchrodinger.jl" begin
    # Tests on maximal intensity generating functions
    m = 0.2
    λ = λ_given_m(m, q=6)
    @test λ ≈ 0.9339432380535286im
    λ = λ_maximal(λ, 3, m=m)
    @test λ ≈ [0.9339432380535284im, 0.8929328208629054im, 0.819955288011106im]
    @test_throws ArgumentError λ_maximal(0.6im, 10, m=m)

    @testset "Initial Conditions" begin
        xᵣ = -10=>10
        λ = 0.75im
        λ, T, Ω = params(λ = λ)
        box = Box(xᵣ, T, Nₓ=1, Nₜ = 100, n_periods = 3)
        coeff = [(2.7 + 4.6im)*1e-2]
        ψ₀, A₀ = ψ₀_periodic(coeff, box, Ω)
        @test A₀ ≈ 0.9971509414326399
    end

    @testset "params" begin
        λ, T, Ω = params(a=3/8)
        @test imag(λ) ≈ sqrt(3)/2
        @test T ≈ 2*π
        @test Ω ≈ 1
        λ, T, Ω = params(λ = sqrt(3)/2*im)
        @test T ≈ 2*π
        @test Ω ≈ 1
        λ, T, Ω = params(T = 2*pi)
        @test imag(λ) ≈ sqrt(3)/2
        @test Ω ≈ 1
        λ, T, Ω = params(Ω = 1)
        @test imag(λ) ≈ sqrt(3)/2
        @test T ≈ 2*π
    end

    @testset "Darboux Transformations" begin
        @testset "Soliton 1" begin
            T = 20
            xᵣ = -20=>20
            seed = "0"
            box = Box(xᵣ, T, Nₓ=1000, Nₜ = 1024, n_periods = 1)
            λ = [0.85im]
            xₛ = [0.0]
            tₛ = [0.0]

            calc = Calc(λ, tₛ, xₛ, seed, box)
            solve!(calc)
            λ = 0.85im;
            ψ₀ = 2*imag(λ)./cosh.(2*imag(λ).*box.t)
            @test abs.(calc.ψ[:, end]) ≈ abs.(ψ₀)
        end
        @testset "Soliton 3" begin
            T = 20
            xᵣ = -20=>20
            seed = "0"
            box = Box(xᵣ, T, Nₓ=1000, Nₜ = 1024, n_periods = 1)
            λ = [-0.9 + 0.8im, 0.85im, 0.9 + 0.9im]
            xₛ = [0.0, 0.0, 0.0]
            tₛ = [0.0, 0.0, 0.0]

            calc = Calc(λ, tₛ, xₛ, seed, box)
            solve!(calc)
            peak = PHF(calc)
            @test peak ≈ 5.1

            λ = 0.85im;
            ψ₀ = 2*imag(λ)./cosh.(2*imag(λ).*box.t)
            # Test around edges and exact peak to avoid weird behavior near peak from DT iteration
            @test abs.(calc.ψ[1:250, end]) ≈ abs.(ψ₀[1:250]) atol=1e-3
            @test abs.(calc.ψ[750:end, end]) ≈ abs.(ψ₀[750:end]) atol=1e-3
            @test abs.(calc.ψ[513, end]) ≈ abs.(ψ₀[513]) atol=1e-2
        end
        @testset "AB 1" begin
            T = 20
            xᵣ = 0=>1e-3
            seed = "exp"
            box = Box(xᵣ, T, Nₓ=1, Nₜ = 1024, n_periods = 1)
            λ = [0.9im]
            xₛ = [0.0]
            tₛ = [0.0]

            calc = Calc(λ, tₛ, xₛ, seed, box)
            solve!(calc)
            peak = PHF(calc)
            @test peak ≈ 2.8

            λ = 0.9im;
            a = imag(λ)^2/2
            Ω = 2*sqrt(1 - 2*a)
            ψ₀ = ((1 - 4*a) .+ (sqrt(2*a).*cos.(Ω*calc.box.t)))./
                (sqrt(2*a).*cos.(Ω*calc.box.t) .- 1)

            @test abs.(calc.ψ[:, 1]) ≈ abs.(ψ₀)
        end
        @testset "CN 1" begin
            T = 20
            xᵣ = -10=>10
            box = Box(xᵣ, T, Nₓ=1, Nₜ = 1024, n_periods = 1)
            λ = [0.7im]
            xₛ = [0.0]
            tₛ = [0.0]

            seed = "cn"
            calc = Calc(λ, tₛ, xₛ, seed, box, m = 0.5)
            peak = PHF(calc)
            @test peak ≈ 2.1071067811865474

            solve!(calc)

            # Test the solution is mirrored across the x=0 axis
            @test abs.(calc.ψ[:,2]) ≈ abs.(calc.ψ[:,1])
        end
        @testset "DN 1" begin
            T = 20
            xᵣ = -10=>10
            box = Box(xᵣ, T, Nₓ=1, Nₜ = 1024, n_periods = 1)
            λ = [0.6im]
            xₛ = [0.0]
            tₛ = [0.0]

            seed = "dn"
            calc = Calc(λ, tₛ, xₛ, seed, box, m = 0.5)
            peak = PHF(calc)
            @test peak ≈ 2.2

            solve!(calc)

            # Test the solution is mirrored across the x=0 axis
            @test abs.(calc.ψ[:,2]) ≈ abs.(calc.ψ[:,1])
        end
    end

    # Simulation Tests using Solitons (easiest way to check validity)
    @testset "Soliton Simulations" begin
        λ = 0.75im
        T = 20
        xᵣ = 0=>10
        box = Box(xᵣ, T, dx=1e-3, Nₜ = 256, n_periods = 1)
        ψ₀ = Array{Complex{Float64}}(undef, box.Nₜ)
        ψ₀ .= 2*imag(λ)./cosh.(2*imag(λ).*box.t)

        for algo ∈ (T1A!, T1B!)
            sim = Sim(λ, box, ψ₀, algo)
            @info "Testing Algorithm" algo
            solve!(sim)
            @test abs.(sim.ψ[:, end]) ≈ abs.(ψ₀) atol=1e-2
        end
        for algo ∈ (T2A!, T4A_TJ!, T6A_TJ!, T8A_TJ!, T4A_SF!, T4A_SF!, T6A_SF!, T8A_SF!,
                T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T6A_KLs9!, T8A_Ss15!,
                T8A_KLs17!,
                T2B!, T4B_TJ!, T6B_TJ!, T8B_TJ!, T4B_SF!, T4B_SF!, T6B_SF!, T8B_SF!,
                T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T6B_KLs9!, T8B_Ss15!, T8B_KLs17!)
            sim = Sim(λ, box, ψ₀, algo)
            @info "Testing Algorithm" algo
            solve!(sim)
            @test abs.(sim.ψ[:, end]) ≈ abs.(ψ₀) atol=1e-5
        end
    end
    @testset "Soliton Simulations Hirota" begin
        λ = 0.75im
        T = 20
        xᵣ = 0=>5
        box = Box(xᵣ, T, dx=1e-3, Nₜ = 256, n_periods = 1)
        ψ₀ = Array{Complex{Float64}}(undef, box.Nₜ)
        ψf = Array{Complex{Float64}}(undef, box.Nₜ)
        ψ₀ .= 2*imag(λ)./cosh.(2*imag(λ).*box.t)
        xf = box.x[end]
        α = 0.1
        v = 4*α*imag(λ)^2
        ψf .= 2*imag(λ)./cosh.(2*imag(λ).*(box.t .+ v*xf))

        for algo ∈ (T1A_H!, T2A_H!)
            sim = Sim(λ, box, ψ₀, algo, α = α)
            @info "Testing Algorithm" algo
            solve!(sim)
            # Can't test around the peak due to a minor shift between them
            @test abs.(sim.ψ[1:50, end]) ≈ abs.(ψ₀[1:50]) atol=5e-2
            @test abs.(sim.ψ[175:end, end]) ≈ abs.(ψ₀[175:end]) atol=5e-2
        end
    end
    @testset "Plot Recipes" begin
        λ = 0.75im
        T = 20
        xᵣ = 0=>1
        box = Box(xᵣ, T, dx=1e-3, Nₜ = 256, n_periods = 1)
        ψ₀ = Array{Complex{Float64}}(undef, box.Nₜ)
        ψ₀ .= 2*imag(λ)./cosh.(2*imag(λ).*box.t)
        sim = Sim(λ, box, ψ₀, T1A!)
        solve!(sim)
        compute_IoM!(sim)
        xₛ = Int(ceil(sim.box.Nₓ/500))
        tₛ = Int(ceil(sim.box.Nₜ/500))
        x_ax = sim.box.t[1:tₛ:end]
        y_ax =  sim.box.x[1:xₛ:end]
        ω_ax = sim.box.ω[1:tₛ:end]
        ψ = sim.ψ[1:tₛ:end, 1:xₛ:end]
        ψ̃ = sim.ψ̃[1:tₛ:end, 1:xₛ:end]

        # Test plotting ψ
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), sim, :ψ)
        @test getfield(rec[1], 2)[1] == x_ax
        @test getfield(rec[1], 2)[2] == y_ax
        @test getfield(rec[1], 2)[3] == abs.(ψ')

        # Test plotting ψ̃
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), sim, :ψ̃)
        @test getfield(rec[1], 2)[1] == y_ax
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(:seriestype => :heatmap), sim, :ψ̃)
        @test getfield(rec[1], 2)[1] == ω_ax
        @test getfield(rec[1], 2)[2] == y_ax
        @test getfield(rec[1], 2)[3] == log.(abs.(ψ̃'))

        # Test plotting IoM
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), sim, :IoM)
        @test getfield(rec[1], 2)[1] == y_ax
        @test getfield(rec[2], 2)[1] == y_ax
        @test getfield(rec[3], 2)[1] == y_ax
        @test getfield(rec[4], 2)[1] == y_ax

        # Test Error
        @test_logs (:error, "Unknown Plotting Mode for NonlinearSchrodinger.jl Objects IoMM") RecipesBase.apply_recipe(Dict{Symbol, Any}(), sim, :IoMM)

    end
end


text:
using Documenter, MatrixPencils
DocMeta.setdocmeta!(MatrixPencils, :DocTestSetup, :(using MatrixPencils); recursive=true)

makedocs(
  modules  = [MatrixPencils],
  sitename = "MatrixPencils.jl",
  authors  = "Andreas Varga",
  format   = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"   => "index.md",
     "Library" => [
        "klftools.md",
        "sklftools.md",
        "klfapps.md",
        "sklfapps.md",
        "pregular.md",
        "lstools.md",
        "pmtools.md",
        "pmapps.md",
        "rmtools.md",
        "rmapps.md"
     ],
     "Internal" => [
        "klftools_int.md"
     ],
     "Index" => "makeindex.md"
  ]
)

deploydocs(
  repo = "github.com/andreasvarga/MatrixPencils.jl.git",
  target = "build",
)


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

__precompile__(true)

module IntervalRootFinding

using IntervalArithmetic
using ForwardDiff
using StaticArrays

import Polynomials: roots

## Root finding
export
    derivative, jacobian,  # reexport derivative from ForwardDiff
    Root, is_unique,
    roots, find_roots,
    bisect, newton1d

import Base: ⊆, show

const Interval = IntervalArithmetic.Interval

const derivative = ForwardDiff.derivative
const D = derivative

# Root object:
immutable Root{T}
    interval::T
    status::Symbol
end

show(io::IO, root::Root) = print(io, "Root($(root.interval), :$(root.status))")

is_unique{T}(root::Root{T}) = root.status == :unique

⊆(a::Interval, b::Root) = a ⊆ b.interval   # the Root object has the interval in the first entry
⊆(a::Root, b::Root) = a.interval ⊆ b.interval

# Common functionality:

doc"""Returns the midpoint of the interval x, slightly shifted in case
the midpoint is an exact root"""
function guarded_mid{T}(f::Function, x::Interval{T})
    m = mid(x)

    if f(m) == zero(T)                      # midpoint exactly a root
        α = convert(T, 0.46875)      # close to 0.5, but exactly representable as a floating point
        m = α*x.lo + (one(T)-α)*x.hi   # displace to another point in the interval
    end

    @assert m ∈ x

    return m
end

include("bisect.jl")
include("bisection.jl")
include("newton.jl")
include("krawczyk.jl")
include("complex.jl")
include("branch_and_prune.jl")
include("newton1d.jl")

gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))
const ∇ = gradient

export ∇

function find_roots{T}(f::Function, a::Interval{T}, method::Function = newton;
                    tolerance = eps(T), debug = false, maxlevel = 30)

    method(f, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
end

function find_roots{T}(f::Function, f_prime::Function, a::Interval{T}, method::Function=newton;
                    tolerance=eps(T), debug=false, maxlevel=30)

    method(f, f_prime, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
end

function find_roots(f::Function, a::Real, b::Real, method::Function=newton;
           tolerance=eps(1.0*a), debug=false, maxlevel=30, precision::Int=-1)

    if precision >= 0
        setprecision(Interval, precision) do
            find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
        end

    else  # use current precision

        find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)

    end
end

function find_roots_midpoint(f::Function, a::Real, b::Real, method::Function=newton;
           tolerance=eps(1.0*a), debug=false, maxlevel=30, precision=-1)

    roots = find_roots(f, a, b, method; tolerance=tolerance, debug=debug, maxlevel=maxlevel, precision=precision)

    T = eltype(roots[1].interval)

    midpoints = T[]
    radii = T[]
    root_symbols = Symbol[]  # :unique or :unknown

    if length(roots) == 0
        return (midpoints, radii, root_symbols)  # still empty
    end

    for root in roots
        midpoint, radius = midpoint_radius(root.interval)
        push!(midpoints, midpoint)
        push!(radii, radius)

        push!(root_symbols, root.status)

    end

    (midpoints, radii, root_symbols)

end

function Base.lexcmp{T}(a::Interval{T}, b::Interval{T})
    #@show a, b
    if a.lo < b.lo
        return -1
    elseif a.lo == b.lo
        if a.hi < b.hi
            return -1
        else
            return 0
        end
    end
    return 1

end

Base.lexcmp{T}(a::Root{T}, b::Root{T}) = lexcmp(a.interval, b.interval)

function clean_roots(f, roots)

    # order, remove duplicates, and include intervals X only if f(X) contains 0
    sort!(roots, lt=lexless)
    roots = unique(roots)
    roots = filter(x -> 0 ∈ f(x.interval), roots)

    # merge neighbouring roots if they touch:

    if length(roots) < 2
        return roots
    end

    new_roots = eltype(roots)[]

    base_root = roots[1]

    for i in 2:length(roots)
        current_root = roots[i]

        if isempty(base_root.interval ∩ current_root.interval) ||
                (base_root.status != current_root.status)

            push!(new_roots, base_root)
            base_root = current_root
        else
            base_root = Root(hull(base_root.interval, current_root.interval), base_root.status)
        end
    end

    push!(new_roots, base_root)

    return new_roots

end

end


text:
#helper functions for pre-processing data
using DataFrames

COCO_CLASSES = ["person", "bicycle", "car", "motorcycle",
			"airplane", "bus", "train", "truck", "boat", "traffic light", "fire hydrant",
			"N/A", "stop sign","parking meter", "bench", "bird", "cat", "dog", "horse",
			"sheep", "cow", "elephant", "bear", "zebra", "giraffe", "N/A", "backpack",
			"umbrella", "N/A", "N/A", "handbag", "tie", "suitcase", "frisbee", "skis", "snowboard",
			"sports ball", "kite", "baseball bat", "baseball glove", "skateboard",
			"surfboard", "tennis racket","bottle", "N/A", "wine glass", "cup", "fork", "knife",
			"spoon","bowl", "banana", "apple", "sandwich", "orange", "broccoli", "carrot",
			"hot dog", "pizza","donut", "cake", "chair", "couch", "potted plant", "bed",
			"N/A", "dining table","N/A", "N/A", "toilet", "N/A", "tv", "laptop", "mouse",
			"remote", "keyboard", "cell phone","microwave", "oven", "toaster", "sink", "refrigerator", "N/A",
			"book","clock", "vase", "scissors", "teddy bear", "hair drier",
			"toothbrush"]

#for now, remove books
#office_subset = ["chair", "keyboard", "laptop", "dining table", "potted plant", "cell phone", "bottle"]
office_subset = ["dining table", "microwave", "backpack", "bowl", "couch"]

# dictionary_scenenet_to_office = Dict("chair" => "chair",
# "straightchair" => "chair", "swivelchair" => "chair", "windsorchair" => "chair",
# "cantileverchair" => "chair", "lawnchair" => "chair", "armchair" => "chair",
# "armchair" => "chair",
# "keyboard" => "keyboard", "computerkeyboard" => "keyboard",
# "laptop" => "laptop", "draftingtable" => "dining table",
# "table" => "dining table",
# "plant" => "potted plant",
# "cellphone" => "cell phone", "cellulartelephone" => "cell phone",
# "bottle" => "bottle", "winebottle" => "bottle")

#dictionary_gt_to_office = Dict("chair" => "chair", "sofa" => "couch")
#dictionary_gt_to_office = Dict("chair" => "chair", "microwave" => "microwave", "sofa" => "couch")
dictionary_gt_to_office = Dict("table" => "dining table", "microwave" => "microwave", "backpack" => "backpack", "bowl" => "bowl", "sofa" => "couch")

#helper function for getting ground-truth from dictionary
function get_ground_truth(dict::Array{Any}, num_videos::Int64)
    world_states = Array{Any}(undef, num_videos)
    for v = 1:num_videos
		#println(v)
        world_state = []
        for item = 1:length(dict[v]["labels"])
			scene_net_name = dict[v]["labels"][item]["category_name"]
			#println(scene_net_name)
			office_name = get(dictionary_gt_to_office, scene_net_name, "NA") #office name will be NA if it's not an entry in the dictionary
			#println(office_name)
			label = findfirst(office_subset .== office_name)
            location = dict[v]["labels"][item]["position"]
			#println(location)
            if !isnothing(label)
                push!(world_state, (location[1], location[2], location[3], label))
            end
        end
        world_states[v] = world_state
    end
    return world_states
end

################################################################################
function extract_category(world_states::Any)
    categories_only = deepcopy(world_states)
    for v = 1:num_videos
        for i = 1:length(world_states[v])
            categories_only[v][i] = world_states[v][i][4] #get just the category
        end
    end
    return categories_only
end

################################################################################
function extract_category(data::DataFrame, num_videos::Int64, num_particles::Int64)
    categories_only = deepcopy(data)

	temp = Array{Any}(undef, num_videos)
	for v = 1:num_videos
		temp[v] = last.(data[v, "inferred_best_world_state"])
    end
	categories_only[!,"inferred_best_world_state"] = temp

	for j = 1:num_particles
		temp = Array{Any}(undef, num_videos) #length(online_names) should be the same as num_videos
		str = "world_state_" * string(j)
    	for v = 1:num_videos
        	temp[v] = last.(data[v, str])
		end
		categories_only[!, str] = temp
    end
    return categories_only
end

################################################################################
function threeD2twoD_new(threeDobjects::Vector{Any}, v::Int64, f::Int64,
		params::Video_Params, camera_trajectories::Matrix{Camera_Params})
	paramses = fill(params, length(threeDobjects))
	camera_paramses = fill(camera_trajectories[v, f], length(threeDobjects))
	threeDobjects_2D = map(render, paramses, camera_paramses, threeDobjects)
	threeDobjects_2D = Array{Detection2D}(threeDobjects_2D)
	threeDobjects_2D = filter(p -> within_frame(p), threeDobjects_2D)
	categories = last.(threeDobjects_2D)
end

################################################################################
function threeD2twoD_new(threeDobjects::Vector{Tuple{Float64, Float64, Float64, Int64}}, v::Int64, f::Int64,
		params::Video_Params, camera_trajectories::Matrix{Camera_Params})
	paramses = fill(params, length(threeDobjects))
	camera_paramses = fill(camera_trajectories[v, f], length(threeDobjects))
	threeDobjects_2D = map(render, paramses, camera_paramses, threeDobjects)
	threeDobjects_2D = Array{Detection2D}(threeDobjects_2D)
	threeDobjects_2D = filter(p -> within_frame(p), threeDobjects_2D)
	categories = last.(threeDobjects_2D)
end

################################################################################
#given a string, parse the columns of the dataframe with that string and return
#a vector of world states
function parse_data(data::DataFrame, model_name::String, num_videos::Int64)
    names_list = []
    for i = 1:length(names(data))
        if occursin(model_name, names(data)[i])
            push!(names_list, names(data)[i])
        end
    end

    world_states = Array{Any}(undef, length(names_list)) #length(online_names) should be the same as num_videos
    for i = 1:length(names_list)
        world_states[i] = eval(Meta.parse(data[1, Symbol(names_list[i])])) #1 is for the first row. we only have 1 row
    end
    return world_states
end

################################################################################
#parse the dataframe and return a vector of world states
function new_parse_data(data::DataFrame, num_videos::Int64)
    world_states = Array{Any}(undef, num_videos) #length(online_names) should be the same as num_videos
    for i = 1:num_videos
        world_states[i] = eval(Meta.parse(data[i, "inferred_best_world_state"])) #1 is for the first row. we only have 1 row
    end
    return world_states
end

################################################################################
#parse the dataframe and return a dataframe of world states
function new_parse_data(data::DataFrame, num_videos::Int64, num_particles::Int64)
	world_states = Array{Any}(undef, num_videos) #length(online_names) should be the same as num_videos
	for i = 1:num_videos
        world_states[i] = eval(Meta.parse(data[i, "inferred_best_world_state"])) #1 is for the first row. we only have 1 row
    end
	data[!, "inferred_best_world_state"] = world_states

	for j = 1:num_particles
		temp = Array{Any}(undef, num_videos) #length(online_names) should be the same as num_videos
		str = "world_state_" * string(j)
		for i = 1:num_videos
	        temp[i] = eval(Meta.parse(data[i, str])) #1 is for the first row. we only have 1 row
	    end
		data[!, str] = temp
	end

    return data
end

################################################################################
#Jaccard similarity
function jaccard_similarity(a::Union{Array{Any}, Vector{Int64}}, b::Union{Array{Any}, Vector{Int64}})
    overlap_of_sets = []
    union_of_sets = []
    for i = 1:length(a)
        index = findfirst(b .== a[i])
        if !isnothing(index)
            #remove from b, add to overlap
            entry = splice!(b, index)
            push!(overlap_of_sets, entry)
        end
        #add a value to union
        push!(union_of_sets, a[i])
    end

    #add whatever is left of b
    for j = 1:length(b)
        push!(union_of_sets, b[j])
    end

    #what happens when both a and b were empty?
    if length(union_of_sets) == 0
        return 1
    else
        return length(overlap_of_sets)/length(union_of_sets)
    end
end

################################################################################
#Jaccard similarity for each column of data compared to ground_truth
function jaccard_similarity(ground_truth_categories::Array{Any}, data::DataFrame, num_particles::Int64)
	similarity_data = deepcopy(data)
	similarity_data[!, "inferred_best_world_state"] = map(jaccard_similarity, ground_truth_categories, data[!, "inferred_best_world_state"])
	for j=1:num_particles
		str = "world_state_" * string(j)
		similarity_data[!, str] = map(jaccard_similarity, ground_truth_categories, data[!, str])
	end
	return similarity_data
end

################################################################################
#function for rangling the dataframe.
#returns a dataframe just like t
function add_confidence_interval(df::DataFrame, n_boot::Int64, cil::Float64)
	from = findfirst(names(df).=="world_state_1")
	to = length(names(df))

	data = select(df, from:to)#now we have just the columns we need
	num_particles = convert(Int64, size(data)[2]/2)
	num_videos = size(data)[1]
	lower_ci = Array{Float64}(undef, num_videos)
	upper_ci = Array{Float64}(undef, num_videos)
	for v = 1:num_videos
		weights = Array{Float64}(undef, num_particles)
		world_states = Array{Float64}(undef, num_particles)
		for j = 1:num_particles
			str = "weight_" * string(j)
			weights[j] = data[v, str] #why 1 and not v?
			str = "world_state_" * string(j)
			world_states[j] = data[v, str] #why 1 and not v
		end
		new_df = DataFrame(weights = weights, values = world_states)
		bs = bootstrap(highest_weighted_world_state, new_df, BasicSampling(n_boot))
		bci = confint(bs, PercentileConfInt(cil))
		lower_ci[v] = bci[1][2]
		upper_ci[v] = bci[1][3]
	end

	#make the dataframe I want to return
	return df[!, "inferred_best_world_state"], lower_ci, upper_ci
end

################################################################################
#this one is for getting sim for the NN
function jacccard_sim_2D(num_videos::Int64, num_frames::Int64, params::Video_Params, dict::Any,
    ground_truth_world_states::Vector{Any}, threshold::Float64, top_n = Inf)

    receptive_fields = make_receptive_fields(params)

    objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold, top_n)

    sim_NN = zeros(num_videos)
    for v = 1:num_videos
        for f = 1:num_frames
            gt_categories = threeD2twoD_new(ground_truth_world_states[v], v, f, params, camera_trajectories)
            NN_categories = last.(objects_observed[v,f])
            sim_NN[v] = sim_NN[v] + jaccard_similarity(gt_categories, NN_categories)
        end
        sim_NN[v] = sim_NN[v]/num_frames
    end
    return sim_NN
end

################################################################################
#this one is for getting sim and confidence interval for a DataFrame
function jacccard_sim_2D(num_videos::Int64, num_frames::Int64, params::Video_Params, dict::Any,
    ground_truth_world_states::Vector{Any}, df::DataFrame, num_bootstraps::Int64, cil::Float64)

    receptive_fields = make_receptive_fields(params)

    objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold)

    sim = zeros(num_videos)
    for v = 1:num_videos
        for f = 1:num_frames
			gt_categories = threeD2twoD_new(ground_truth_world_states[v], v, f, params, camera_trajectories)
			model_categories = threeD2twoD_new(df[v, "inferred_best_world_state"], v, f, params, camera_trajectories)
            sim[v] = sim[v] + jaccard_similarity(gt_categories, model_categories)
        end
        sim[v] = sim[v]/num_frames
    end

	#let's make confidence intervals
	from = findfirst(names(df).=="world_state_1")
	to = length(names(df))
	data = select(df, from:to)#now we have just the columns we need
	num_particles = convert(Int64, size(data)[2]/2)
	lower_ci = Array{Float64}(undef, num_videos)
	upper_ci = Array{Float64}(undef, num_videos)
	sim_values = zeros(num_videos, num_particles)

	#bootstrapping happens for each video
	for v = 1:num_videos
		weights = Array{Float64}(undef, num_particles)
		for j = 1:num_particles
			str = "weight_" * string(j)
			weights[j] = df[v, str]
			for f = 1:num_frames
				gt_categories = threeD2twoD_new(ground_truth_world_states[v], v, f, params, camera_trajectories)
				str = "world_state_" * string(j)
				model_categories = threeD2twoD_new(df[v, str], v, f, params, camera_trajectories)
				sim_values[v, j] = sim_values[v, j] + jaccard_similarity(gt_categories, model_categories)
			end
			sim_values[v,j] = sim_values[v,j]/num_frames
		end
		new_df = DataFrame(weights = weights, values = sim_values[v,:])
		#println("new_df ", new_df)
		bs = bootstrap(highest_weighted_world_state, new_df, BasicSampling(num_bootstraps))
		#println("bs ", bs)
		bci = confint(bs, PercentileConfInt(cil))
		#println("bci ", bci)
		lower_ci[v] = bci[1][2]
		upper_ci[v] = bci[1][3]
	end
    return sim, lower_ci, upper_ci
end

################################################################################
# #some data is a dataframe with just one row. it has num_particles values and num_particles weights
# #returns the value of the highest weighted particle
# function highest_weighted_world_state(some_data::DataFrame)
# 	num_particles = convert(Int64, size(some_data)[2]/2)
# 	world_states = Array{Float64}(undef, num_particles)
# 	weights = Array{Float64}(undef, num_particles)
# 	for j = 1:num_particles
# 		str = "world_state_" * string(j)
# 		world_states[j] = some_data[1, str]
# 		str = "weight_" * string(j)
# 		weights[j] = some_data[1, str]
# 	end
# 	best_world_state = world_states[findmax(weights)[2]]
# 	return best_world_state
# end

#some data is a dataframe with two columns (sim and weights) and num_particles rows.
#returns the value of the highest weighted particle
function highest_weighted_world_state(some_data::DataFrame)
	row = findmax(some_data[!, "weights"])[2]
	return some_data[row, "values"]
end

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

function within_frame(x::Float64, y::Float64)
    x >= 0 && x <= 256 && y >= 0 && y <= 256 #hard-codded frame size
end

function within_frame(p::Detection2D)
    p[1] >= 0 && p[1] <= 256 && p[2] >= 0 && p[2] <= 256 #hard-codded frame size
end

################################################################################
#for making ground-truth demo
function write_gt_to_dict(dict::Array{Any,1}, camera_trajectories::Matrix{Camera_Params}, gt_objects::Vector{Any})
    params = Video_Params()
	num_videos, num_frames = size(camera_trajectories)

    #add inferences about the objects in the scenes
    for v = 1:num_videos
		#println("v ", v)
        #for each object
        gt = gt_objects[v]

        #add per frame
        for f = 1:num_frames
			#println("f ", f)
            labels = []
            centers = []
            for i in 1:length(gt)
                r = gt[i]
                x, y = get_image_xy(camera_trajectories[v,f], params, Coordinate(r[1],r[2],r[3]))
                if within_frame(x, y)
					#println(r)
                    push!(labels, r[4])
                    push!(centers, (x, y))
                end
            end
            d_f = dict[v]["views"][f]
            a = Dict("labels" => labels, "centers" => centers)
            merge!(d_f, Dict("ground_truth" => a))
        end
    end

    return dict
end

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
for cid in ls("^c")
 !isa(eval(cid), CORE) && continue
 println("$cid = $(length(contradictionsof(eval(cid))))")
end
println()

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
"""
    struct BlockEnsemble{T}
        recordings::Vector{T}
    end

Gathers all blocks relevant for an entire single recording.

    BlockEnsemble(
        se::SamplingEnsemble,
        ranges,
        ρρ=0.0,
        ll_hist_len=0
    )

Base constructor.
"""
struct BlockEnsemble{T}
    recordings::Vector{T}

    function BlockEnsemble(
            se::SamplingEnsemble,
            ranges,
            ρρ=0.0,
            ll_hist_len=0
        )
        num_rec = num_recordings(se)
        ρρ = _vec_me(ρρ, num_rec)
        ll_hist_len = _vec_me(ll_hist_len, num_rec)
        rec = map(1:num_rec) do i
            BlockCollection(se.recordings[i], ranges[i], ρρ[i], ll_hist_len[i])
        end
        new{eltype(rec)}(rec)
    end
end

OBS.num_recordings(aob::BlockEnsemble) = length(aob.recordings)

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

        IMPUTATION

===============================================================================#
"""
    draw_proposal_path!(be::BlockEnsemble)

Call `draw_proposal_path!` separately for each recording
"""
draw_proposal_path!(be::BlockEnsemble) = draw_proposal_path!.(be.recordings)

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

        ACCEPT/REJECT DECISION

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

"""
    accept_reject_proposal_path!(be::BlockEnsemble, mcmciter)

Call `accept_reject_proposal_path!` separately for each recording.
"""
function accept_reject_proposal_path!(be::BlockEnsemble, mcmciter)
    for rec in be.recordings
        accept_reject_proposal_path!(rec, mcmciter)
    end
end

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

        SWAPS

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

"""
    swap_paths!(be::BlockEnsemble)

Call `swap_paths!` separately for each recording.
"""
swap_paths!(be::BlockEnsemble) = swap_paths!.(be.recordings)

"""
    swap_XX!(be::BlockEnsemble)

Call `swap_XX!` separately for each recording.
"""
swap_XX!(be::BlockEnsemble) = swap_XX!.(be.recordings)

"""
    swap_WW!(be::BlockEnsemble)

Call `swap_WW!` separately for each recording.
"""
swap_WW!(be::BlockEnsemble) = swap_WW!.(be.recordings)

"""
    swap_PP!(be::BlockEnsemble)

Call `swap_PP!` separately for each recording.
"""
swap_PP!(be::BlockEnsemble) = swap_PP!.(be.recordings)

"""
    swap_ll!(be::BlockEnsemble)

Call `swap_ll!` separately for each recording.
"""
swap_ll!(be::BlockEnsemble) = swap_ll!.(be.recordings)

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

        UTILITY

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

"""
    loglikhd!(be::BlockEnsemble)

Call `loglikhd!` separately for each recording.
"""
loglikhd!(be::BlockEnsemble) = loglikhd!.(be.recordings)

"""
    loglikhd°!(be::BlockEnsemble)

Call `loglikhd°!` separately for each recording.
"""
loglikhd°!(be::BlockEnsemble) = loglikhd°!.(be.recordings)

@doc raw"""
    fetch_ll(be::BlockEnsemble)

Retreive the log-likelihood for all accepted paths
!!! warning
    The function uses only internal fields `ll` for this computation, which
    means that for the call to this function to make sense the log-likelihood
    must have been previously computed and stored in the field `ll`. If it
    hasn't been done, then you must first call `loglikhd!(be)`.
"""
fetch_ll(be::BlockEnsemble) = mapreduce(rec->fetch_ll(rec), +, be.recordings)

@doc raw"""
    fetch_ll°(be::BlockEnsemble)

Retreive the log-likelihood for all proposed paths
!!! warning
    The function uses only internal fields `ll` for this computation, which
    means that for the call to this function to make sense the log-likelihood
    must have been previously computed and stored in the field `ll`. If it
    hasn't been done, then you must first call `loglikhd°!(be)`.
"""
fetch_ll°(be::BlockEnsemble) = mapreduce(rec->fetch_ll°(rec), +, be.recordings)

"""
    save_ll!(be::BlockEnsemble, i::Int)

Call `save_ll!` separately for each recording.
"""
save_ll!(be::BlockEnsemble, i::Int) = save_ll!.(be.recordings, i)

"""
    ll_of_accepted(be::BlockEnsemble, i)

Return an array of log-likelihoods (one for each block in each recording) of the
paths that were accepted at the `i`th iteration.
"""
ll_of_accepted(be::BlockEnsemble, i) = ll_of_accepted.(be.recordings, i)

"""
    accpt_rate(be::BlockEnsemble, range)

Compute the acceptance rate over the `range` of MCMC accept/reject history for
each block for each recording.
"""
function accpt_rate(be::BlockEnsemble, range)
    map(be.recordings) do rec
        accpt_rate(rec, range)
    end
end

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

        SETTING ARTIFICIAL OBSERVATIONS

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

"""
    GP.set_obs!(be::BlockEnsemble)

Call `save_ll!` separately for each recording.
"""
GP.set_obs!(be::BlockEnsemble) = set_obs!.(be.recordings)

"""
    GP.recompute_guiding_term!(be::BlockEnsemble, [::Val{:_only}])

Call `recompute_guiding_term!` separately for each recording.
"""
function GP.recompute_guiding_term!(be::BlockEnsemble)
    recompute_guiding_term!.(be.recordings)
end

function GP.recompute_guiding_term!(be::BlockEnsemble, v::Val)
    recompute_guiding_term!.(be.recordings, v)
end

"""
    find_W_for_X!(be::BlockEnsemble)

Call `find_W_for_X!` separately for each recording.
"""
find_W_for_X!(be::BlockEnsemble) = find_W_for_X!.(be.recordings)

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

        SETTING PARAMETERS

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

"""
    is_critical_update(be::BlockEnsemble, pnames)

Call `is_critical_update` separately for each recording.
"""
GP.is_critical_update(be::BlockEnsemble, pnames) = map(
    rec->GP.is_critical_update(rec, pnames),
    be.recordings
)

"""
    set_proposal_law!(
        be::BlockEnsemble,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bc, pnames);
        skip=0
    )

Call `set_proposal_law!` separately for each recording.
"""
function set_proposal_law!(
        be::BlockEnsemble,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bc, pnames);
        skip=0
    )
    for i in eachindex(be.recordings)
        set_proposal_law!(
            be.recordings[i], θ°, pnames.recordings[i], critical_change[i];
            skip=skip
        )
    end
end


text:
#= src/psvd.jl
=#

mutable struct PartialSVD{T,Tr<:Real} <: Factorization{T}
  U::Matrix{T}
  S::Vector{Tr}
  Vt::Matrix{T}

  function PartialSVD{T,Tr}(U::Matrix{T}, S::Vector{Tr}, Vt::Matrix{T}) where {T,Tr<:Real}
    size(U,2) == size(Vt,1) == length(S) || throw(DimensionMismatch("$(size(U)), $(length(S)), $(size(Vt)) not compatible"))
    new{T,Tr}(U, S, Vt)
  end
end

PartialSVD(U::AbstractMatrix{T}, S::AbstractVector{Tr}, Vt::AbstractMatrix{T}) where {T,Tr<:Real} =
  PartialSVD{T,Tr}(Matrix(U), Vector(S), Matrix(Vt))

conj!(A::PartialSVD) = PartialSVD(conj!(A.U), A.S, conj!(A.Vt))
conj(A::PartialSVD) = PartialSVD(conj(A.U), A.S, conj(A.Vt))

function convert(::Type{PartialSVD{T}}, A::PartialSVD) where T
  Tr = real(T)
  PartialSVD(
    convert(Array{T}, A.U), convert(Array{Tr}, A.S), convert(Array{T}, A.Vt))
end
convert(::Type{Factorization{T}}, A::PartialSVD) where {T} = convert(PartialSVD{T}, A)
convert(::Type{Factorization{T}}, A::PartialSVD{T,Tr}) where {T,Tr<:Real} = convert(PartialSVD{T}, A)
convert(::Type{Array}, A::PartialSVD) = Matrix(A)
convert(::Type{Array{T}}, A::PartialSVD) where {T} = convert(Array{T}, Matrix(A))

copy(A::PartialSVD) = PartialSVD(copy(A.U), copy(A.S), copy(A.Vt))

adjoint!(A::PartialSVD) = PartialSVD(A.Vt', A.S, A.U')
adjoint(A::PartialSVD) = PartialSVD(A.Vt', copy(A.S), A.U')
transpose!(A::PartialSVD) = PartialSVD(transpose(A.Vt), A.S, transpose(A.U))
transpose(A::PartialSVD) = PartialSVD(transpose(A.Vt), copy(A.S), transpose(A.U))

Matrix(A::PartialSVD) = (A[:U]*Diagonal(A[:S]))*A[:Vt]

function getindex(A::PartialSVD, d::Symbol)
  if     d == :S   return A.S
  elseif d == :U   return A.U
  elseif d == :V   return A.Vt'
  elseif d == :Vt  return A.Vt
  elseif d == :k   return length(A.S)
  else             throw(KeyError(d))
  end
end

ishermitian(::PartialSVD) = false
issym(::PartialSVD) = false

isreal(::PartialSVD{T}) where {T} = T <: Real

ndims(::PartialSVD) = 2

size(A::PartialSVD) = (size(A.U,1), size(A.Vt,2))
size(A::PartialSVD, dim::Integer) =
  dim == 1 ? size(A.U,1) : (dim == 2 ? size(A.Vt,2) : 1)

# BLAS/LAPACK multiplication/division routines

  ## left-multiplication

  mul!(y::AbstractVector{T}, A::PartialSVD{T}, x::AbstractVector{T}) where {T} =
    mul!(y, A[:U], scalevec!(A[:S], A[:Vt]*x))
  mul!(C::AbstractMatrix{T}, A::PartialSVD{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:U], lmul!(Diagonal(A[:S]), A[:Vt]*B))

  ## right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::PartialSVD{T}) where {T} =
    mul!(C, rmul!(A*B[:U], Diagonal(B[:S])), B[:Vt])

  for Adj in (:Adjoint, :Transpose)
    @eval begin
      function mul!(C::AbstractMatrix{T}, A::PartialSVD{T}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        tmp = A[:Vt] * Bc
        lmul!(Diagonal(A[:S]), tmp)
        mul!(C, A[:U], tmp)
      end
      function mul!(y::AbstractVector{T}, Ac::$Adj{T,<:PartialSVD{T}}, x::AbstractVector{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:U]) * x
        scalevec!(A[:S], tmp)
        mul!(y, $Adj(A[:Vt]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:PartialSVD{T}}, B::AbstractMatrix{T}) where T
        A = parent(Ac)
        tmp = $Adj(A[:U]) * B
        lmul!(Diagonal(A[:S]), tmp)
        mul!(C, $Adj(A[:Vt]), tmp)
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:PartialSVD{T}}, Bc::$Adj{T,<:AbstractMatrix{T}}) where T
        A = parent(Ac)
        tmp = $Adj(A[:U]) * Bc
        lmul!(Diagonal(A[:S]), tmp)
        mul!(C, $Adj(A[:Vt]), tmp)
      end
      ## right-multiplication
      function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::$Adj{T,<:PartialSVD{T}}) where T
        B = parent(Bc)
        tmp = A * $Adj(B[:Vt])
        rmul!(tmp, Diagonal(B[:S]))
        mul!(C, tmp, $Adj(B[:U]))
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::PartialSVD{T}) where T
        tmp = Ac * B[:U]
        rmul!(tmp, Diagonal(B[:S]))
        mul!(C, tmp, B[:Vt])
      end
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, Bc::$Adj{T,<:PartialSVD{T}}) where T
        tmp = Ac * $Adj(B[:Vt])
        rmul!(tmp, Diagonal(B[:S]))
        mul!(C, tmp, $Adj(B[:U]))
      end
    end
  end

  ## left-division (pseudoinverse left-multiplication)
  ldiv!(y::AbstractVector{T}, A::PartialSVD{T}, x::AbstractVector{T}) where {T} =
    mul!(y, A[:Vt]', iscalevec!(A[:S], A[:U]'*x))
  ldiv!(C::AbstractMatrix{T}, A::PartialSVD{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:Vt]', iscale!(A[:S], A[:U]'*B))

  # standard operations

  ## left-multiplication
  function *(A::PartialSVD{TA}, B::AbstractVector{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialSVD{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1))
    mul!(CT, AT, BT)
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function *(Ac::$Adj{TA,<:PartialSVD{TA}}, B::AbstractVector{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialSVD{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1))
        mul!(CT, $Adj(AT), BT)
      end
    end
  end

  function *(A::PartialSVD{TA}, B::AbstractMatrix{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialSVD{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end
  function *(A::AbstractMatrix{TA}, B::PartialSVD{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = (T == TA ? A : convert(Array{T}, A))
    BT = convert(PartialSVD{T}, B)
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end
  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function *(A::PartialSVD{TA}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialSVD{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:PartialSVD{TA}}, B::AbstractMatrix{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialSVD{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:PartialSVD{TA}}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialSVD{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
      function *(A::AbstractMatrix{TA}, Bc::$Adj{TB,<:PartialSVD{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialSVD{T}, B)
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, B::PartialSVD{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialSVD{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, Bc::$Adj{TB,<:PartialSVD{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialSVD{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
    end
  end

## left-division
function \(A::PartialSVD{TA}, B::AbstractVector{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialSVD{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2))
  ldiv!(CT, AT, BT)
end
function \(A::PartialSVD{TA}, B::AbstractMatrix{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialSVD{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2), size(B,2))
  ldiv!(CT, AT, BT)
end

# factorization routines

function psvdfact(
    A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
  opts = isempty(args) ? opts : copy(opts; args...)
  m, n = size(A)
  if m >= n
    V = idfact(:n, A, opts)
    Q,R = qr!(getcols(:n, A, V[:sk]))
    Ũ, σ, Ṽ = svd!(R*V)
    k = psvdrank(σ, opts)
    if k < V[:k]
      U  = Q*view(Ũ,:,1:k)
      S  = σ[1:k]
      Vt = Ṽ'[1:k,:]
    else
      U  = Q*Ũ
      S  = σ
      Vt = Ṽ'
    end
  else
    V = idfact(:c, A, opts)
    Q,R = qr!(getcols(:c, A, V[:sk]))
    Ũ, σ, Ṽ = svd!(V'*R')
    k = psvdrank(σ, opts)
    if k < V[:k]
      U  = Ũ[:,1:k]
      S  = σ[1:k]
      Vt = view(Ṽ',1:k,:)*Q'
    else
      U  = Ũ
      S  = σ
      Vt = Ṽ'*Q'
    end
  end
  PartialSVD(U, S, Vt)
end

function psvdvals(
    A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
  opts = isempty(args) ? opts : copy(opts; args...)
  m, n = size(A)
  if m >= n
    V = idfact(:n, A, opts)
    Q,R = qr!(getcols(:n, A, V[:sk]))
    s = svdvals!(R*V)
  else
    V = idfact(:c, A, opts)
    Q,R = qr!(getcols(:c, A, V[:sk]))
    s = svdvals!(V'*R')
  end
  k = psvdrank(s, opts)
  k < V[:k] && return s[1:k]
  s
end

for f in (:psvdfact, :psvdvals)
  @eval $f(A, args...; kwargs...) = $f(LinOp(A), args...; kwargs...)
end

function psvd(A, args...; kwargs...)
  F = psvdfact(A, args...; kwargs...)
  F.U, F.S, F.Vt'
end

function psvdrank(s::Vector{T}, opts::LRAOptions) where T<:Real
  k = length(s)
  ptol = max(opts.atol, opts.rtol*s[1])
  @inbounds for i = 2:k
    s[i] <= ptol && return i - 1
  end
  k
end