text:
module Test_fditools

using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Testing examples for EFDSYN
println("Test_fdiutils")
@testset "fdiutils" begin

# setup the system with additive faults
# [Gu(s) Gd(s) Gf(s)], where Gf(s) = [ Gu(s) I]
s = rtf('s'); # define the Laplace variable s
Gu = [(s+1)/(s+2); (s+2)/(s+3)] # enter Gu(s)
Gd = [(s-1)/(s+2); 0]; # enter Gd(s)
# build state space model of [Gu(s) Gd(s) Gf(s)] and set input groups
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 1:2)
display(sysf)

# tests FDFilter
filter = rss(3,1,3)
Q = FDFilter(filter,1,2)
Q2 = FDFilter(filter; outputs = [1],controls = [2,3])
@test iszero(Q.sys-Q2.sys,atol=1.e-7) &&
            Q.controls == Q2.controls && Q.outputs == Q2.outputs
Q3 = FDFilter(filter,outputs = 1,controls = 2:3)
@test iszero(Q.sys-Q3.sys,atol=1.e-7) &&
            Q.controls == Q3.controls && Q.outputs == Q3.outputs
display(Q)

# tests FDIFilter
filter = [rss(3,1,3), rss(2,2,3)]
Q = FDIFilter(filter,1,2)
Q2 = FDIFilter(filter; outputs = [1],controls = [2,3])
@test iszero(Q.sys[1]-Q2.sys[1],atol=1.e-7) &&  iszero(Q.sys[2]-Q2.sys[2],atol=1.e-7) &&
            Q.controls == Q2.controls && Q.outputs == Q2.outputs
Q3 = FDIFilter(filter,outputs = 1,controls = 2:3)
@test iszero(Q.sys[1]-Q3.sys[1],atol=1.e-7) &&  iszero(Q.sys[2]-Q3.sys[2],atol=1.e-7) &&
            Q.controls == Q3.controls && Q.outputs == Q3.outputs
display(Q)

# tests FDFilterIF
filteri = rss(3,1,6)
R = FDFilterIF(filteri,controls=1:2,disturbances=[3],aux=6,noise=5,faults=[4])
R1 = FDFilterIF(filteri,2,1,1,1,1)
@test iszero(R.sys-R1.sys,atol=1.e-7) &&
       R.controls == R1.controls && R.disturbances == R1.disturbances &&
       R.faults == R1.faults && R.noise == R1.noise &&
       R.aux == R1.aux

R = FDFilterIF(filteri,aux=[6],noise=[5],faults=[4])
R1 = FDFilterIF(filteri,0,0,1,1,1;moff = 3)
@test iszero(R.sys-R1.sys,atol=1.e-7) &&
       R.controls == R1.controls && R.disturbances == R1.disturbances &&
       R.faults == R1.faults && R.noise == R1.noise &&
       R.aux == R1.aux

R2 = FDFilterIF(filteri[:,4:end],0,0,1,1,1)
@test iszero(R.sys-R2.sys,atol=1.e-7) &&
       R.controls == R2.controls && R.disturbances == R2.disturbances &&
       R.faults == R2.faults && R.noise == R2.noise &&
       R.aux == R2.aux

# tests FDIFilterIF
filteri = [rss(3,1,6), rss(2,2,6)]
R = FDIFilterIF(filteri,controls=1:2,disturbances=[3],aux=6,noise=5,faults=[4])
R1 = FDIFilterIF(filteri,2,1,1,1,1)
@test iszero(R.sys[1]-R1.sys[1],atol=1.e-7) && iszero(R.sys[2]-R1.sys[2],atol=1.e-7) &&
       R.controls == R1.controls && R.disturbances == R1.disturbances &&
       R.faults == R1.faults && R.noise == R1.noise &&
       R.aux == R1.aux
display(R)

R = FDIFilterIF(filteri,aux=[6],noise=[5],faults=[4])
R1 = FDIFilterIF(filteri,0,0,1,1,1;moff = 3)
@test iszero(R.sys[1]-R1.sys[1],atol=1.e-7) && iszero(R.sys[2]-R1.sys[2],atol=1.e-7) &&
       R.controls == R1.controls && R.disturbances == R1.disturbances &&
       R.faults == R1.faults && R.noise == R1.noise &&
       R.aux == R1.aux

sysf1 = fdimodset(rss(3,1,6,stable=true),c = 1:2, d = 3, f = 4:5, n = 6, aux = 2:6)
Q = FDFilter(rss(3,1,3,stable=true),1,2)
R = fdIFeval(Q, sysf1)
syse = [sysf1.sys;eye(2,11)]
@test iszero(Q.sys*syse[:,1:2]-R.sys[:,R.controls],atol=1.e-7) &&
      iszero(Q.sys*syse[:,3]-R.sys[:,R.disturbances],atol=1.e-7) &&
      iszero(Q.sys*syse[:,4:5]-R.sys[:,R.faults],atol=1.e-7) &&
      iszero(Q.sys*syse[:,6]-R.sys[:,R.noise],atol=1.e-7) &&
      iszero(Q.sys*syse[:,7:end]-R.sys[:,R.aux],atol=1.e-7)
@test fditspec(R) == fditspec_(R.sys[:,R.faults])
@time scond, β, γ = fdiscond(R)
@test β == fdhinfminus(R.sys[:,R.faults])[1] && γ == fdhinfmax(R.sys[:,R.faults])[1]
freq = 0
@time scond, β, γ = fdiscond(R,freq)
@test β == fdhinfminus(R.sys[:,R.faults],freq)[1] && γ == fdhinfmax(R.sys[:,R.faults],freq)[1]
freq = [0, 3, 5];
@time scond, β, γ = fdiscond(R,freq)
@test β == fdhinfminus(R.sys[:,R.faults],freq)[1] && γ == fdhinfmax(R.sys[:,R.faults],freq)[1]

sysf1 = fdimodset(rss(3,1,6,stable=true),c = 1:2, d = 3, f = 4:5, n = 6, aux = 2:6)
Q = FDIFilter([rss(3,1,3,stable=true),rss(2,2,3,stable=true)],1,2)
R = fdIFeval(Q, sysf1)
syse = [sysf1.sys;eye(2,11)]
@test all(iszero.(Q.sys .* [syse]-R.sys,atol=1.e-7))

@time S = fdisspec(R,[0,1]; FDGainTol = 1.e-5)
@test S ==  ([1  1 ; 1  1] .> 0)
@time S1 = fditspec(R; FDfreq = [0,1], FDtol = 1.e-5)
@test S == S1
@time S2 = fditspec(R)
@test S2 == ([1  1 ; 1  1] .> 0)

@test fditspec_(sysf.sys[:,sysf.faults],block=true,FDtol = 1.e-3) == ([1  1  1] .> 0)
@test fditspec_(sysf.sys[:,sysf.faults],block=true,FDfreq=0)[:,:,1] == ([1  1  1] .> 0)
@test fditspec_(sysf.sys[:,sysf.faults],block=true,FDfreq=0)[:,:,1] == ([1  1  1] .> 0)
@test fditspec_(sysf.sys[:,sysf.faults],block=false) == ([1  1  0; 1  0  1] .> 0)
@test fditspec_(sysf.sys[:,sysf.faults],block=false,FDfreq=0)[:,:,1] == ([1  1  0; 1  0  1] .> 0)
@test fditspec_(sysf.sys[:,sysf.faults],block=false,FDfreq=[0,1,Inf])[:,:,1] == ([1  1  0; 1  0  1] .> 0)

sysf1 = fdimodset(rss(3,1,6,stable=true),c = 1:2, d = 3, f = 4:5, n = 6, aux = 2:6)
sysf2 = fdimodset(rss(2,2,6,stable=true),c = 1:2, d = 3, f = 4:5, n = 6, aux = 2:6)

@time S, gains = fdisspec_(sysf.sys[:,sysf.faults],[0,1,Inf]; block=false)
@test S[:,:,1] == ([1  1  0; 1  0  1] .> 0) &&
     gains[:,:,1] ≈  [0.5  1.0  0.0; 2/3  0.0  1.0] && gains[:,:,3] ≈  [1.0 1.0 0.0; 1.0 0.0 1.0] && gains[:,:,2] ≈  [0.6324555320336759 1.0 0.0; 0.7071067811865475 0.0 1.0]

@time S1, gains1 = fdisspec_(sysf.sys[:,sysf.faults],[0,1,Inf]; block=true)
@test S1[:,:,1] == Bool[1 1 1] && gains1[:,:,1] ≈  [0.8333333333333334 1.0 1.0]

# some non-trivial examples for frequency values equal to poles
# S1 = [1 0]
s = rtf('s')
g = dss([1/(s^2+1) 1/s]);
S1, g1 = fdisspec_(g, [0,1]; block=true,stabilize = true)
S2 = fditspec_(g,FDfreq = [0,1],block=true,poleshift = true)
@test S1 == S2

S1, g1 = fdisspec_(g, [0,1]; block=false,stabilize = true)
S2 = fditspec_(g,FDfreq =[0,1],block=false,poleshift = true)
@test S1 == S2

# S1 = [1 0]
s = rtf('s')
g = dss([1/s 1/(s+2)]);
S1, g1 = fdisspec_(g,block=true,stabilize = true)
S2 = fditspec_(g,FDfreq = 0,block=true,poleshift = true)
@test S1 == S2

S1, g1 = fdisspec_(g,block=false,stabilize = true)
S2 = fditspec_(g,FDfreq = 0,block=false,poleshift = true)
@test S1 == S2

# S1 = [1 1]
s = rtf('s')
g = dss([1/s 1/s/(s+2)]);
S1, g1 = fdisspec_(g,block=true,stabilize = true)
S2 = fditspec_(g,FDfreq = 0,block=true,poleshift = true)
@test S1 == S2

S1, g1 = fdisspec_(g,block=false,stabilize = true)
S2 = fditspec_(g,FDfreq = 0,block=false,poleshift = true)
@test S1 == S2

s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)
# new information in smat and zero local H-minus index values

smat, gains = fdisspec_(dss(Gf), [0, 2]; FDGainTol = 1.e-3)
@test smat[:,:,1] == Bool[0 0; 1 1] && smat[:,:,2] == Bool[1 0; 0 1] &&
      gains[:,:,1] ≈ [0.0 0.0; 0.4444444444444444 1.0] && gains[:,:,2] ≈ [0.894427190999916 0.0; 0.0 1.0]
ssmat = fditspec_(dss(Gf), atol = 1.e-7,FDtol = 1.e-3, FDfreq = [0, 2])
@test smat == ssmat
swmat = fditspec_(dss(Gf), atol=1.e-7,FDtol = 1.e-3)
@test swmat == Bool[1 0; 1 1]

s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [1/s; 1/(s^2+4)]; # enter Gf(s)

smat, gains = fdisspec_(dss(Gf), [0, 2]; FDGainTol = 1.e-3, stabilize = true)
@test smat[:,:,1] == trues(2,1) && smat[:,:,2] == trues(2,1) && all(maximum(gains,dims=1) .> 0.001)
ssmat = fditspec_(dss(Gf), atol = 1.e-7,FDtol = 1.e-3, FDfreq = [0, 2],poleshift = true)
@test smat == ssmat
swmat = fditspec_(dss(Gf), atol=1.e-7,FDtol = 1.e-3)
@test swmat == trues(2,1)

s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s; s^2+4]; # enter Gf(s)

smat, gains = fdisspec_(dss(Gf), [0, 2]; FDGainTol = 1.e-3,stabilize = true)
@test smat[:,:,1] == eye(2)[:,2:2] && smat[:,:,2] == eye(2)[:,1:1] &&
      gains[1,1,1] < 1.e-7 && gains[2,1,2] < 1.e-7
ssmat = fditspec_(dss(Gf), atol = 1.e-7,FDtol = 1.e-3, FDfreq = [0, 2],poleshift = true)
@test smat == ssmat
swmat = fditspec_(dss(Gf), atol=1.e-7,FDtol = 1.e-3)
@test swmat ==trues(2,1)

end # fditools

end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 3, Zadanie 6. metoda Newtona, funkcja `f_1` — edge case dla `x_0 ∈ (1, ∞)`
#

include("iterative-methods.jl")
using .IterativeMethods

fl = IterativeMethods.FL

a = fl(10)

r, fr, k, err =  IterativeMethods.mstycznych(
    x -> exp(1-x) - 1,
    x -> -1.0 * exp(1 - x),
    a,
    fl(10^(-5)),
    fl(10^(-5)),
    100
)

println(
    "r = ", r, "\n",
    "f(r) = ", fr, "\n",
    "k = ", k, "\n",
    err == 0 ? "" : string("error code: ", err)
)

ext: jl

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
export BestResponsePolicy

using Flux:onehot

struct BestResponsePolicy{E, S, A, X, P<:AbstractPolicy} <: AbstractCFRPolicy
    cfr_reach_prob::Dict{S, Vector{Pair{E, Float64}}}
    best_response_action_cache::Dict{S,A}
    best_response_value_cache::Dict{E,Float64}
    best_responder::X
    policy::P
end

"""
    BestResponsePolicy(policy, env, best_responder)

- `policy`, the original policy to be wrapped in the best response policy.
- `env`, the environment to handle.
- `best_responder`, the player to choose best response action.
"""
function BestResponsePolicy(policy, env, best_responder; state_type=String, action_type=Int)
    # S = typeof(get_state(env))  # TODO: currently it will break the OpenSpielEnv. Can not get information set for chance player
    # A = eltype(get_actions(env))  # TODO: for chance players it will return ActionProbPair
    S = state_type
    A = action_type
    E = typeof(env)

    p = BestResponsePolicy(
        Dict{S, Vector{Pair{E, Float64}}}(),
        Dict{S, A}(),
        Dict{E, Float64}(),
        best_responder,
        policy
    )

    e = copy(env)
    @assert e == env  "The copy method doesn't seem to be implemented for environment: $env"
    @assert hash(e) == hash(env) "The hash method doesn't seem to be implemented for environment: $env"
    RLBase.reset!(e)  # start from the root!
    init_cfr_reach_prob!(p, e)
    p
end

function (p::BestResponsePolicy)(env::AbstractEnv)
    if get_current_player(env) == p.best_responder
        best_response_action(p, env)
    else
        p.policy(env)
    end
end

function init_cfr_reach_prob!(p, env, reach_prob=1.0)
    if !get_terminal(env)
        if get_current_player(env) == p.best_responder
            push!(get!(p.cfr_reach_prob, get_state(env), []), env => reach_prob)

            for a in get_legal_actions(env)
                init_cfr_reach_prob!(p, child(env, a), reach_prob)
            end
        elseif get_current_player(env) == get_chance_player(env)
            for a::ActionProbPair in get_actions(env)
                init_cfr_reach_prob!(p, child(env, a), reach_prob * a.prob)
            end
        else  # opponents
            for a in get_legal_actions(env)
                init_cfr_reach_prob!(p, child(env, a), reach_prob * get_prob(p.policy, env, a))
            end
        end
    end
end

function best_response_value(p, env)
    get!(p.best_response_value_cache, env) do
        if get_terminal(env)
            get_reward(env, p.best_responder)
        elseif get_current_player(env) == p.best_responder
                a = best_response_action(p, env)
                best_response_value(p, child(env, a))
        elseif get_current_player(env) == get_chance_player(env)
            v = 0.
            for a::ActionProbPair in get_actions(env)
                v += a.prob * best_response_value(p, child(env, a))
            end
            v
        else
            v = 0.
            for a in get_legal_actions(env)
                v += get_prob(p.policy, env, a) * best_response_value(p, child(env, a))
            end
            v
        end
    end
end

function best_response_action(p, env)
    get!(p.best_response_action_cache, get_state(env)) do
        best_action, best_action_value = nothing, typemin(Float64)
        for a in get_legal_actions(env)
            # for each information set (`get_state(env)` here), we may have several paths to reach it
            # here we sum the cfr reach prob weighted value to find out the best action
            v = sum(p.cfr_reach_prob[get_state(env)]) do (e, reach_prob)
                reach_prob * best_response_value(p, child(e, a))
            end
            if v > best_action_value
                best_action, best_action_value = a, v
            end
        end
        best_action
    end
end

RLBase.update!(p::BestResponsePolicy, args...) = nothing

function RLBase.get_prob(p::BestResponsePolicy, env::AbstractEnv)
    if get_current_player(env) == p.best_responder
        onehot(p(env), get_actions(env))
    else
        get_prob(p.policy, env)
    end
end


text:
using Printf

export run_sim

# ============================================================================ #
"""
run_sim(d::LiveDemo, amp::Real=1.0, dur::Real=4.0, nrep::Integer=10)
"""
function run_sim(d::LiveDemo, amp::Real=1.0, dur::Real=4.0, nrep::Integer=10)
    stim = SineWave(d, amp, 4.0)
    return run_sim(d.net, stim, dur, nrep)
end
"""
run_sim(d::LiveDemo, stim::Stimulus, dur::Real=4.0, nrep::Integer=10)
"""
function run_sim(d::LiveDemo, stim::Stimulus, dur::Real=4.0, nrep::Integer=10)
    return run_sim(d.net, stim, dur, nrep)
end
# ---------------------------------------------------------------------------- #
"""
run_sim(net::LIFNetwork, stim::Stimulus, dur::Real, nrep::Integer,
    tstart::Time=0.0, record_vm::Bool=false)
"""
function run_sim(net::LIFNetwork, stim::Stimulus, dur::Real,
    nrep::Integer, tstart::Time=0.0, record_vm::Bool=false)

    #dur should be in seconds
    if dur > 100.0
        error("Simulation duration is really long... check your units!")
    end

    ncell = length(net)

    evt_ts = tstart:dur:tstart + dur*(nrep-1)

    fstim(id::Integer, t::Float64) = getstim(stim, id, t)
    # begin
        # if id == 1
        #     return amp*((0.5*sin(2.0*pi*(t*1e-3)*4.0))+0.5)
        # else
        #     return 0.0
        # end
    # end

    ts = [Vector{Float64}() for x in 1:ncell]

    vm = [zeros(1) for x in 1:ncell]
    if record_vm
        rec = collect(1:ncell)
    else
        rec = Int[]
    end

    @inbounds for kr = 1:nrep

        reset!(stim)
        Networks.reset!(net)
        d = simulate(net, fstim, dur*1e3, rec)

        @inbounds for k = 1:ncell
            append!(ts[k], d[:ts][k] .+ evt_ts[kr])
            if record_vm
                append!(vm[k], d[:vm][:,k])
            end
        end
    end

    return ts, vm, evt_ts
end
# =========================================================================== #


text:
@testset "Utilities" begin
  a, b, c = P2(0,0), P2(1,0), P2(0,1)
  @test signarea(a, b, c) == T(0.5)
  a, b, c = P2(0,0), P2(0,1), P2(1,0)
  @test signarea(a, b, c) == T(-0.5)

  p1, p2, p3 = P2(0,0), P2(1,1), P2(0.25,0.5)
  s = Segment(P2(0.5,0.0), P2(0.0,1.0))
  @test sideof(p1, s) == :LEFT
  @test sideof(p2, s) == :RIGHT
  @test sideof(p3, s) == :ON

  p1, p2, p3 = P2(0.5,0.5), P2(1.5,0.5), P2(1,1)
  c = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
  @test sideof(p1, c) == :INSIDE
  @test sideof(p2, c) == :OUTSIDE
  @test sideof(p3, c) == :INSIDE

  # drop units from unitful value and type
  @test Meshes.dropunits(1.0u"mm") == Float64
  @test Meshes.dropunits(typeof(1.0u"mm")) == Float64

  # return the same type in case of no units
  @test Meshes.dropunits(1.0) == Float64
  @test Meshes.dropunits(Float64) == Float64
end


text:
# MPO for the triangular AF Ising
T = mpo_triangular_AF_ising()

T_arr = reshape(T.data, (2, 2, 2, 2))

@timedtestset "test act ($ix)" for ix in 1:10

    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    psi_arr = reshape(psi.data, (5, 2, 5))
    Tpsi = act(T, psi)
    Tpsi_arr = reshape(ein"acmd,bme->abcde"(T_arr, psi_arr), (20, 10))
    @test Tpsi_arr ≈ Tpsi.data

end

# test rrule of act
@timedtestset "test rrule for act ($ix)" for ix in 1:10
    function f_with_OMEinsum(psi_arr::Array{ComplexF64, 3})
        Tpsi_arr = reshape(ein"acmd,bme->abcde"(T_arr, psi_arr), (10, 2, 10))
        result = ein"abc,abc->"(conj(Tpsi_arr), Tpsi_arr) |> real
        return result[1]
    end
    function f(psi::TensorMap{ComplexSpace, 2, 1})
        Tpsi = act(T, psi)
        return dot(Tpsi, Tpsi) |> real
    end

    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    psi_arr = reshape(psi.data, (5, 2, 5))
    @test f(psi) ≈ f_with_OMEinsum(psi_arr)

    grad_von_OMEinsum = gradient(f_with_OMEinsum, psi_arr)[1]
    grad = reshape(gradient(f, psi)[1].data, (5, 2, 5))
    @test grad_von_OMEinsum ≈ grad
end

@timedtestset "test transf_mat constructions ($ix)" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    phi = TensorMap(rand, ComplexF64, ℂ^6*ℂ^2, ℂ^6)
    v = TensorMap(rand, ComplexF64, ℂ^6, ℂ^5)
    lop = transf_mat(psi, phi)
    lop_T = transf_mat_T(psi, phi)

    psi_arr = toarray(psi)
    phi_arr = toarray(phi)
    v_arr = toarray(v)
    lop_arr = ein"aec,bed->abcd"(conj.(psi_arr), phi_arr)

    lop_v = ein"abcd,dc->ba"(lop_arr, v_arr)
    lopT_v = ein"abcd,ba->dc"(conj.(lop_arr), v_arr)

    @test lop_v ≈ lop(v) |> toarray
    @test lopT_v ≈ lop_T(v) |> toarray
end

@timedtestset "dominant vr and vl should be positive definite" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    v = TensorMap(rand, ComplexF64, ℂ^5, ℂ^5)
    lop = transf_mat(psi, psi)
    lop_T = transf_mat_T(psi, psi)

    wr, vr = eigsolve(lop, v, 1)
    wl, vl = eigsolve(lop_T, v, 1)

    @test wr[1] ≈ wl[1]

    vr, vl = vr[1], vl[1]

    # vr, vl is only positive definite up to a phase
    u1, _, v1 = tsvd(vr)
    u2, _, v2 = tsvd(vl)
    @test u1.data ≈ (u1[1,1]/v1[1,1]') * v1.data'
    @test u2.data ≈ (u2[1,1]/v2[1,1]') * v2.data'

end

@timedtestset "test ovlp and its rrule" for ix in 1:10
    function ovlp_von_arr(psi_arr::Array{ComplexF64, 3}, phi_arr::Array{ComplexF64, 3})
        lop_arr = ein"aec,bed->abcd"(conj.(psi_arr), phi_arr)
        lop_arr = reshape(lop_arr, (30, 30))
        w_von_arr = eigvals(lop_arr) |> last
        return w_von_arr
    end

    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    phi = TensorMap(rand, ComplexF64, ℂ^6*ℂ^2, ℂ^6)
    psi_arr = toarray(psi)
    phi_arr = toarray(phi)

    @test ovlp(psi, phi) ≈ ovlp_von_arr(psi_arr, phi_arr)

    psi_pushback1 = gradient(arr->real(ovlp_von_arr(arr, phi_arr)), psi_arr)[1]
    psi_pushback2 = gradient(arr->real(ovlp(arr_to_TensorMap(arr), phi)), psi_arr)[1]
    @test psi_pushback1 ≈ psi_pushback2

    phi_pushback1 = gradient(arr->real(ovlp_von_arr(psi_arr, arr)), phi_arr)[1]
    phi_pushback2 = gradient(arr->real(ovlp(psi, arr_to_TensorMap(arr))), phi_arr)[1]
    @test phi_pushback1 ≈ phi_pushback2

end

@timedtestset "test Gamma-Lambda conversion" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    norm_psi, Γ, Λ = lambda_gamma(psi)

    @test ovlp(psi, psi) ≈ norm_psi^2

    @tensor L[-1; -2] := Γ'[-1, 3, 4] * Λ'[3, 2] * Λ[2, 1] * Γ[1, 4, -2]
    @tensor R[-1; -2] := Γ'[3, -2, 4] * Λ'[2, 3] * Λ[1, 2] * Γ[-1, 4, 1]
    Id = id(ℂ^5)
    @test L ≈ Id
    @test R ≈ Id
end

@timedtestset "test mps add" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    phi = TensorMap(rand, ComplexF64, ℂ^6*ℂ^2, ℂ^6)

    psi_plus_phi = mps_add(psi, phi)
    @test toarray(psi_plus_phi)[1:5, :, 1:5] ≈ toarray(psi)
    @test toarray(psi_plus_phi)[6:11, :, 6:11] ≈ toarray(phi)
    @test toarray(psi_plus_phi)[1:5, :, 6:11] ≈ zeros(5, 2, 6)
    @test toarray(psi_plus_phi)[6:11, :, 1:5] ≈ zeros(6, 2, 5)
end

@timedtestset "test right_canonical_QR" for ix in 1:10
    chi = 10
    psi = TensorMap(rand, ComplexF64, ℂ^chi*ℂ^2, ℂ^chi)
    Y, psi_R = right_canonical_QR(psi)

    # test the iterated equation
    @tensor rhs[-1, -2; -3] := Y'[-1, 1] * psi_R[1, -2, -3]
    lhs = psi * Y'
    @test lhs ≈ (lhs[1] / rhs[1]) * rhs

    # test the right-canonical relation
    psi_Rtmp = permute(psi_R, (1,), (2, 3))
    @test psi_Rtmp * psi_Rtmp' ≈ id(ℂ^chi)
    # test fidelity
    @test ln_fidelity(psi, psi_R) > -1e-14
end

@timedtestset "test left_canonical_QR" for ix in 1:10
    chi = 10
    psi = TensorMap(rand, ComplexF64, ℂ^chi*ℂ^2, ℂ^chi)
    X, psi_L = left_canonical_QR(psi)

    # test the iterated equation
    @tensor lhs[-1, -2; -3] := X[-1, 1] * psi[1, -2, -3]
    rhs = psi_L * X
    @test lhs ≈ (lhs[1] / rhs[1]) * rhs

    # test the left canonical relation
    @test psi_L' * psi_L ≈ id(ℂ^chi)

    # test the fidelity
    @test ln_fidelity(psi, psi_L) > -1e-14
end


text:
# Return a random DNA/RNA sequence of the given length.
function random_seq(n::Integer, nts, probs, outtype = String)
    cumprobs = cumsum(probs)
    x = Vector{Char}(undef, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return outtype(x)
end

function random_seq(::Type{A}, n::Integer) where {A<:Alphabet}
    # TODO: Resolve the use of symbols(A()).
    nts = symbols(A())
    probs = Vector{Float64}(undef, length(nts))
    fill!(probs, 1 / length(nts))
    return LongSequence{A}(random_seq(n, nts, probs))
end

function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

function random_dna_symbols(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs, Vector{DNA})
end

function random_rna_symbols(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs, Vector{RNA})
end

function random_rna_symbols(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs, Vector{RNA})
end

function random_aa_symbols(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
     'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ], probs, Vector{AminoAcid})
end

function random_dna_kmer(len)
    return random_dna(len, [0.25, 0.25, 0.25, 0.25])
end

function random_rna_kmer(len)
    return random_rna(len, [0.25, 0.25, 0.25, 0.25])
end

function dna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    complementer = Dict(zip("-ACGTSWYRKMDVHBN", "-TGCASWRYMKHBDVN"))
    for (i, c) in enumerate(seq)
        seqc[i] = complementer[c]
    end
    return String(seqc)
end

function rna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    complementer = Dict(zip("-ACGUSWYRKMDVHBN", "-UGCASWRYMKHBDVN"))
    for (i, c) in enumerate(seq)
        seqc[i] = complementer[c]
    end
    return String(seqc)
end


text:
"""
    Hierachical drawing of directed graphs inspired by the Sugiyama framework.
    In particular see Chapter 13 of 'Hierachical Drawing Algorithms' from
    the 'Handbook of Graph Drawing and Visualization' and the article
      K. Sugiyama, S. Tagawa, and M. Toda. Methods for visual understanding
      of hierarchical system structures. IEEE Transaction on Systems, Man,
      and Cybernetics, 11(2):109–125, 1981.

    The method as implemented here has 4 steps:
    1. Cycle removal [if needed]
    2. Layer assignment + break up long edges
    3. Vertex ordering [to reduce crossings]
    4. Vertex coordinates [to straighten edges]
    5. Draw the tree

    Will return the Compose.jl visualization of the tree, and optionally
    save to file (SVG) if a filename is provided.

    Arguments:
    adj_list        Directed graph in adjacency list format
    labels          Label for each vertex

    Optional arguments:
    filename        Output filename
    cycles          If false, assume no cycles. Default true.
    ordering        Vertex ordering method to use. Options are:
                        :optimal        Uses JuMP (integer program)
                        :barycentric    Sugiyama heuristic
    coord           Vertex coordinate method to use. Options are:
                        :optimal        Uses JuMP (linear program)

    xsep            Controls the minimum vertex horizontal spacing
    ysep            Controls the minimum vertex vertical spacing
    scale           Scales the output figure size
    labelpad        Padding around text in vertices
    background      Defaults to nothing. If not-nothing, that value
                    will be used for the background color
"""
function layout_tree{T}(adj_list::AdjList{T},
                        labels::Vector;
                        filename    = "",

                        cycles      = true,
                        ordering    = :optimal,
                        coord       = :optimal,
                        xsep        = 3,
                        ysep        = 20,
                        scale       = 0.05,
                        labelpad    = 1.2,
                        background  = nothing)
    # Calculate the original number of vertices
    n = length(adj_list)

    # 1     Cycle removal
    if cycles
        # Need to remove cycles first
        error("Cycle removal not implemented!")
    end

    # 2     Layering
    # 2.1   Assign a layer to each vertex
    layers = _layer_assmt_longestpath(adj_list)
    num_layers = maximum(layers)
    # 2.2  Create dummy vertices for long edges
    adj_list, layers = _layer_assmt_dummy(adj_list, layers)
    orig_n, n = n, length(adj_list)

    # 3     Vertex ordering [to reduce crossings]
    # 3.1   Build initial permutation vectors
    layer_verts = [L => Int[] for L in 1:num_layers]
    for i in 1:n
        push!(layer_verts[layers[i]], i)
    end
    # 3.2  Reorder permutations to reduce crossings
    if ordering == :barycentric
        layer_verts = _ordering_barycentric(adj_list, layers, layer_verts)
    elseif ordering == :optimal
        layer_verts = _ordering_ip(adj_list, layers, layer_verts)
    end

    # 4     Vertex coordinates [to straighten edges]
    # 4.1   Place y coordinates in layers
    locs_y = zeros(n)
    for L in 1:num_layers
        for (x,v) in enumerate(layer_verts[L])
            locs_y[v] = (L-1)*ysep
        end
    end
    # 4.2   Get widths of each label, if there are any
    widths  = ones(n); widths[orig_n+1:n]  = 0
    heights = ones(n); heights[orig_n+1:n] = 0
    # Note that we will convert these sizes into "absolute" units
    # and then work in these same units throughout. The font size used
    # here is just arbitrary, and unchanging. This hack arises because it
    # is meaningless to ask for the size of the font in "relative" units
    # but we don't want to collapse to absolute units until the end.
    if length(labels) == orig_n
        extents = text_extents("sans",10pt,labels...)
        for (i,(width,height)) in enumerate(extents)
            widths[i]  = width.value
            heights[i] = height.value
        end
    end
    locs_x = _coord_ip(adj_list, layers, layer_verts, orig_n, widths, xsep)
    # 4.3   Summarize vertex info
    max_x, max_y = maximum(locs_x), maximum(locs_y)
    max_w, max_h = maximum(widths), maximum(heights)

    # 5     Draw the tree
    # 5.1   Create the vertices
    verts = [_tree_textrect(locs_x[i], locs_y[i], labels[i], widths[i], heights[i]) for i in 1:orig_n]
    # 5.2   Create the arrows
    arrows = Any[]
    for L in 1:num_layers, i in layer_verts[L], j in adj_list[i]
        push!(arrows, _arrow_tree(
                locs_x[i], locs_y[i], i<=orig_n ? max_h : 0,
                locs_x[j], locs_y[j], j<=orig_n ? max_h : 0))
    end
    # 5.3   Assemble composition
    # We need to decide the true font size now. The logic is as follows:
    # - All the text has the same height (max_h)
    # - The image is max_h+max_y units tall
    # - The real image will be scale*max_y inches tall
    # - Therefore one unit = ...
    ratio = (scale*max_y)/(max_y+max_h)
    #   inchs.
    # - Now we want to scale the text size, in inches, so that
    #   the height of the text (in inches) is approximately the height
    #   we already have in 'units'. We'll fudge font size = font height,
    #   then use the padding factor to scale the text. Bigger padding,
    #   smaller font relative to its box.
    fsize = (1.0/labelpad) * max_h * (scale*max_y)/(max_y+max_h) * inch
    # Determine the background, if we want one
    bg = background == nothing ? [] : [rectangle(), fill(background)]
    c = compose(
        context(units=UnitBox(-max_w/2,-max_h/2,max_x+max_w,max_y+max_h)),
        font("sans"), fontsize(fsize),
        bg..., verts..., arrows...
    )
    # 5.4   Draw it
    if filename != ""
        draw(SVG(filename, scale*max_x*inch, scale*max_y*inch), c)
    end
    return c
end

@doc """
    Creates a rectangle with text inside it for use in tree drawing.

    Arguments:
    x,y             Center of rectangle/text (in context units)
    label           The text inside the rectangle
    width, height   The size of the rectangle
""" ->
_tree_textrect(x, y, label, width, height) = compose(
        context(x - width/2, y - height/2, width, height),
        (context(), text(0.5w, 0.5h, label, hcenter, vcenter), fill("black")),
        (context(), rectangle(), fill("white"), stroke("black"))
    )

@doc """
    Creates an arrow between two rectangles in the tree

    Arguments:
    o_x, o_y, o_h   Origin x, y, and height
    d_x, d_y, d_h   Destination x, y, and height
""" ->
function _arrow_tree(o_x, o_y, o_h, d_x, d_y, d_h)
    x1, y1 = o_x, o_y + o_h/2
    x2, y2 = d_x, d_y - d_h/2
    Δx, Δy = x2 - x1, y2 - y1
    θ = atan2(Δy, Δx)
    # Put an arrow head only if destination isn't dummy
    head = d_h != 0 ? _arrow_heads(θ, x2, y2, 2) : []
    compose(context(), stroke("black"),
        line([(x1,y1),(x2,y2)]), head...)
end

@doc """
    Creates an arrow head given the angle of the arrow and its destination.

    Arguments:
    θ               Angle of arrow (radians)
    dest_x, dest_y  End of arrow
    λ               Length of arrow head tips
    ϕ               Angle of arrow head tips relative to angle of arrow
""" ->
_arrow_heads(θ, dest_x, dest_y, λ, ϕ=0.125π) = [ line([
    (dest_x - λ*cos(θ+ϕ), dest_y - λ*sin(θ+ϕ)),
    (dest_x, dest_y),
    (dest_x - λ*cos(θ-ϕ), dest_y - λ*sin(θ-ϕ))
]) ]


text:
import AoC2019
using LightGraphs, MetaGraphs

test_input = """
9 ORE => 2 A
8 ORE => 3 B
7 ORE => 5 C
3 A, 4 B => 1 AB
5 B, 7 C => 1 BC
4 C, 1 A => 1 CA
2 AB, 3 BC, 4 CA => 1 FUEL
"""

example_13312 = """
157 ORE => 5 NZVS
165 ORE => 6 DCFZ
44 XJWVT, 5 KHKGT, 1 QDVJ, 29 NZVS, 9 GPVTF, 48 HKGWZ => 1 FUEL
12 HKGWZ, 1 GPVTF, 8 PSHF => 9 QDVJ
179 ORE => 7 PSHF
177 ORE => 5 HKGWZ
7 DCFZ, 7 PSHF => 2 XJWVT
165 ORE => 2 GPVTF
3 DCFZ, 7 NZVS, 5 HKGWZ, 10 PSHF => 8 KHKGT
"""

function create_mapping(lines)
    mapping = Dict()
    names_set = Set{String}()
    min_produce = Dict("ORE" => 1)
    for line in lines
      if line == ""
        continue
      end
      entries = split(line, " => ")
      @assert length(entries) == 2
      inputs = split(entries[1], ", ")
      outputs = entries[2]
      out_count = split(outputs, " ")[1]
      out_name = split(outputs, " ")[2]
      push!(names_set, out_name)
      min_produce[out_name] = parse(UInt8, out_count)
      for inp in inputs
        inp_count = split(inp, " ")[1]
        inp_name = split(inp, " ")[2]
        mapping[(inp_name, out_name)] = parse(UInt8, inp_count)
        push!(names_set, inp_name)
      end
    end
    names=collect(names_set)
    min_produce_list = map(x -> min_produce[x], names)

    return (mapping=mapping, names=names, min_produce = min_produce_list)
end

lines = AoC2019.read_lines(".julia/dev/AoC2019/src/input14.txt")

function generate_fuel(lines, fuel_count=1)
  (mapping, names, min_produce) = create_mapping(lines)
  #println(stderr, names)
  g = MetaDiGraph(length(names))

  for m in keys(mapping)
    #add_edge!(g, findfirst(isequal())
    inp_ind = findfirst(isequal(m[1]), names)
    out_ind = findfirst(isequal(m[2]), names)

    add_edge!(g, inp_ind, out_ind)
    set_props!(g, Edge(inp_ind, out_ind), Dict(:conv =>mapping[m]))

  end

  #start out with Fuel
  fuel_ind = findfirst(isequal("FUEL"), names)
  ore_ind = findfirst(isequal("ORE"), names)
  queue = Array{Tuple{Int64, Int64}, 1}()
  push!(queue, (fuel_count, fuel_ind))

  produced = zeros(Int64, length(names))
  remainder = zeros(Int64, length(names))

  while length(queue) > 0
    current_node = popfirst!(queue)
    id = current_node[2]
    current_required = current_node[1]
    extra_required = current_required - remainder[id]
    neighbors = inneighbors(g, current_node[2])

    producing = convert(Int64, ceil(extra_required / min_produce[id])) * min_produce[id]
    remainder[id] = (producing - extra_required)
    #println("Current: $(names[id]), current_required: $current_required, extra_required: $extra_required, producing: $producing")

    produced[id] += producing

    for n in neighbors
      local conversion_rate = get_prop(g, Edge(n, current_node[2]), :conv)
      #remainder = current_required - producing
      local new_val = producing / min_produce[id] * conversion_rate
      push!(queue, (new_val, n))

    end
  end
  queue_converted = map(x -> (x[1], names[x[2]]), queue)
  #for (i, name) in enumerate(names)
  #  println("Resource $name generated $(produced[i]) times")
  #end
  #for (i, name) in enumerate(names)
  #  println("Resource $name wasted $(remainder[i]) times")
  #end
  return produced[ore_ind]
end

#ore_required = generate_fuel(split(test_input, "\n"))
ore_required = generate_fuel(lines)

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
using Base.Test
using Church
using Distributions
const iter = 10^6
#Geometric distribution

geom() = @If(bernoulli(), 1+geom(), 1)
test_geom() = begin
    a = geom()
    normal(a, 3; condition=15)
    as = zeros(iter)
    for i = 1:iter
        resample()
        as[i] = value(a)
    end
    as
end
as = test_geom()

ms = mean(as)
vs = sqrt(var(as))

#Analytic
n = 30
xs = 1:n
prior = 0.5.^xs
likelihood = pdf(Normal(15, 3), xs)
posterior = prior .* likelihood
norm_posterior = posterior/sum(posterior)
ma = sum(norm_posterior.*xs)
va = sqrt(sum(norm_posterior.*xs.^2) - ma^2)

println((ms, ma))
println((vs, va))

@test ma*0.95 < ms < ma*1.05
@test va*0.95 < vs < va*1.05

gc_church()

#Normal distribution
test_norm() = begin
    a = normal(0, 2)
    b = normal(a, 1; condition=3)
    as = zeros(iter)
    for i = 1:iter
        resample()
        as[i] = value(a)
    end
    as
end
as = test_norm()
ms = mean(as)
vs = var(as)

#Analytic Moments
va = 1/(1/2^2 + 1/1)
ma = (3/1)/(1/2^2+1/1)

@test 0.95*ma < ms < 1.05*ma
@test 0.95*va < vs < 1.05*va

test_mixture() = begin
    mss = zeros(iter, 2)
    n = 10
    data = [randn(n) + 10, randn(n) - 10]
    ms = [normal(0, 10) for i = 1:2]
    ks = [bernoulli() for i = 1:length(data)]
    for i = 1:length(data)
        normal(ms[ks[i]+1]; condition=data[i])
    end
    for i = 1:iter
        mss[i, 1] = value(ms[1])
        mss[i, 2] = value(ms[2])
        resample()
    end
    mss
end
(a, b) = test_mixture()[end, :]

@test sign(a) != sign(b)
@test abs(a) > 8
@test abs(b) > 8
gc_church()


text:
struct FVLAE{E,D,C,F,G,XD,ZD,V,XDIST<:Val} <: AbstractVLAE
    e::E
    d::D
    c::C # critic
    f::F # extracts latent variables
    g::G # concatenates latent vars with the rest
    xdim::XD # (h,w,c)
    zdim::ZD # scalar
    var::V # dense or conv last layer (ignore for dense fvlae)
    xdist::XDIST # gaussian or bernoulli
end
FVLAE(e,d,c,f,g,xdim,zdim,var,xdist::Symbol=:gaussian) =
	FVLAE(e,d,c,f,g,xdim,zdim,var,Val(xdist))
Flux.@functor FVLAE
(m::FVLAE)(x) = reconstruct(m, x)

function Base.show(io::IO, m::FVLAE)
    msg = """$(nameof(typeof(m))):
        encoder =
            e - $(repr(m.e))
            f - $(repr(m.f))
        decoder =
            d - $(repr(m.d))
            g - $(repr(m.g))
        critic =
            $(repr(m.c))
    """
    print(io, msg)
end

"""
	FVLAE(zdim::Int, hdim::Int, ks, ncs, stride, datasize; discriminator_nlayers::Int=3,
		layer_depth::Int=1, var=:dense, activation="relu", xdist=:gaussian)

    FVLAE(zdim::Int, hdims, discriminator_hdim::Int, datasize; discriminator_nlayers::Int=2,
    	activation=relu, layer_depth=1)

Factor VLAE constructor.
"""
function FVLAE(zdim::Int, hdim::Int, ks, ncs, strd, datasize; discriminator_nlayers::Int=3,
		var=:dense, xdist=:gaussian, kwargs...)
	(xdist in [:gaussian, :bernoulli]) ? nothing :
		error("xdist must be either :gaussian or :bernoulli")

	# get encoder, decoder and latent extractors and reshapers
	e,d,f,g = basic_model_constructor(zdim, ks, ncs, strd, datasize; var=var, xdist=xdist,
		kwargs...)

	# now construct the critic (discriminator)
	c = discriminator_constructor(zdim*length(ks), hdim, discriminator_nlayers; kwargs...)

    return FVLAE(e,d,c,f,g,datasize[1:3],zdim,var,xdist)
end
function FVLAE(zdim::Int, hdims, d_hdim::Int, datasize; discriminator_nlayers::Int=2, kwargs...)
    # get encoder, decoder and latent extractors and reshapers
    e,d,f,g = basic_model_constructor(zdim, hdims, datasize[1]; kwargs...)

	# now construct the critic (discriminator)
	c = HierarchicalAD.discriminator_constructor(zdim*length(hdims), d_hdim,
		discriminator_nlayers; kwargs...)

    return HierarchicalAD.FVLAE(e,d,c,f,g,(datasize[1],),zdim,nothing,:gaussian)
end

function logpdf(m::FVLAE{E,D,C,F,G,X,Z,V,XD}, x::AbstractArray{T,2}, zs...) where
		{E,D,C,F,G,X,Z,V,XD<:Val{:gaussian}} where T
    μx, σx = _decoded_mu_var(m, zs...)
    _x = (m.var == :dense) ? vectorize(x) : x
    return  gauss_logpdf(_x, μx, σx)
end
function logpdf(m::FVLAE{E,D,C,F,G,X,Z,V,XD}, x::AbstractArray{T,4}, zs...) where
	{E,D,C,F,G,X,Z,V,XD<:Val{:gaussian}} where T
    μx, σx = _decoded_mu_var1(m, zs...)
    _x = (m.var == :dense) ? vectorize(x) : x
    return  gauss_logpdf(_x, μx, σx)
end
function logpdf(m::FVLAE{E,D,C,F,G,X,Z,V,XD}, x::AbstractArray{T,4}, zs...) where
	{E,D,C,F,G,X,Z,V,XD<:Val{:bernoulli}} where T
    _x = _decoder_out(m, zs...)
    return  bernoulli_logpdf(_x, x)
end

dloss(d,z::AbstractArray{T,2},zp::AbstractArray{T,2}) where T  =
	-0.5f0*Flux.mean(log.(d(z) .+ eps(Float32)) .+ log.(1 .- d(zp) .+ eps(Float32)))
function total_correlation(d,z::AbstractArray{T,2}) where T
	dz = d(z)
	Flux.mean(log.(dz .+ eps(Float32)) .- log.(1 .- dz .+ eps(Float32)))
end

function factor_aeloss(m, x, γ::Float32)
	# encoder pass - KL divergence
	μzs_σzs = _encoded_mu_vars(m, x)
	zs = map(y->rptrick(y...), μzs_σzs)
	kldl = sum(map(y->Flux.mean(kld(y...)), μzs_σzs))

	# decoder pass - logpdf
	lpdf = Flux.mean(logpdf(m, x, zs...))
	elbo = -kldl + lpdf

	# now the discriminator loss
	_zs = cat(zs...,dims=1)
	tc = total_correlation(m.c, _zs)

	-elbo + γ*tc
end

reconstruction_loss(m,x)  =
	-Flux.mean(logpdf(m, x, encode_all(m, x)...))
kld_loss(m, x) =
    kldl = sum(map(y->Flux.mean(kld(y...)), _encoded_mu_vars(m, x)))
tc_loss(m, x) =
	total_correlation(m.c, cat(encode_all(m, x)...,dims=1))

function init_perm_mat(x::AbstractArray{T,2}) where T
    n = size(x,2)
    fill!(similar(x, n, n), 0), n
end
function init_perm_mat(x::SubArray)
    n = size(x,1)
    fill!(similar(x.parent, n, n), 0), n
end
function permutation_mat(x)
	perm_mat, n = init_perm_mat(x)
	for (j,i) in zip(1:n, sample(1:n, n, replace=false))
		perm_mat[j,i] = 1
	end
	perm_mat
end
# permutation for individual latents
permute_cols(x) where T = x*permutation_mat(x)
# permutation for inidivdual dims
permute_rows(x::AbstractArray{T,2}) where T =
    vcat(map(x->reshape(permutation_mat(x)*x,1,size(x,1)), eachrow(x))...)
Flux.Zygote.@nograd permutation_mat # or less?

function factor_closs_per_dim(m, x)
	zs = _encode_all(m, x; mean=false)
	zps = map(zs) do z
		zp = permute_rows(z)
		zp
	end
	_zs = cat(zs..., dims=1)
	_zps = cat(zps..., dims=1)
	dloss(m.c, _zs, _zps)
end
function factor_closs_per_latent(m, x)
	zs = _encode_all(m, x; mean=false)
	zps = map(zs) do z
		zp = permute_cols(z)
		zp
	end
	_zs = cat(zs..., dims=1)
	_zps = cat(zps..., dims=1)
	dloss(m.c, _zs, _zps)
end

#### THIS IS FOR IMAGE DATA ####
"""
	train_fvlae(zdim, hdim, batchsize, ks, ncs, stride, nepochs, tr_x, val_x;
	λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,4), layer_depth=1, lr=0.001f0, var=:dense,
	activation="relu", discriminator_nlayers=3, xdist=:gaussian, pad=0,
	initial_convergence_threshold=0)

Train a factored VLAE.
"""
function train_fvlae(zdim, hdim, batchsize, ks, ncs, strd, nepochs, tr_x::AbstractArray{T,4},
	val_x::AbstractArray{T,4}; tr_y=nothing, val_y=nothing, factors=nothing,
	disentangle_per_latent=true,
	λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,4),
	layer_depth=1, lr=0.001f0, var=:dense,
	activation="relu", discriminator_nlayers=3, xdist=:gaussian, pad=0,
	convergence_threshold=0f0, initial_convergence_epochs=10,
    max_retrain_tries=10, early_stopping=true, kwargs...) where T
    # this is to ensure that the model converges to something meaningful
	hist, rdata, model, zs, aeopt, copt = nothing, nothing, nothing, nothing, nothing, nothing

    ntries = 0

    while isnothing(hist) && ntries<=max_retrain_tries
        if ntries > 0
            @info "Restarting training automatically..."
        end

        # construct the model
        model = gpu(FVLAE(zdim, hdim, ks, ncs, strd, size(tr_x), layer_depth=layer_depth, var=var,
            activation=activation, nlayers=discriminator_nlayers, xdist=xdist, pad=pad))
        aeopt = ADAM(lr)
        copt = ADAM(lr)

        # train it
        hist, rdata, zs = train!(model, nepochs, batchsize, tr_x, val_x, aeopt, copt;
        	tr_y=tr_y, val_y=val_y, factors=factors,
        	disentangle_per_latent=disentangle_per_latent,
            λ=λ, γ=γ, epochsize = epochsize,
            convergence_threshold=convergence_threshold,
            initial_convergence_epochs=initial_convergence_epochs,
            early_stopping=early_stopping)
        ntries += 1
    end

    # if the training was still unsuccesfull
    if isnothing(hist)
        @info "Gave up the training after $ntries tries."
        return model, hist, rdata, zs, aeopt, copt
    end

    return model, hist, rdata, zs, aeopt, copt
end

"""
	train!(model::FVLAE, nepochs, batchsize, tr_x::AbstractArray{T,4},
		val_x::AbstractArray{T,4}, aeopt, copt;
		tr_y=nothing, val_y=nothing, factors=nothing, disentangle_per_latent=true,
		λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,4), convergence_threshold=0f0,
    	initial_convergence_epochs=2)
"""
function train!(model::FVLAE, nepochs, batchsize, tr_x::AbstractArray{T,4},
	val_x::AbstractArray{T,4}, aeopt, copt;
	tr_y=nothing, val_y=nothing, factors=nothing, disentangle_per_latent=true,
	disentanglement_repeats::Int=1,
	λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,4), convergence_threshold=0f0,
    initial_convergence_epochs=10, early_stopping=true,
    kwargs...) where T
	if !isnothing(tr_y) && !(isnothing(val_y)) && isnothing(factors)
		error("specify factors you want to disentangle")
	end

	# data an initial reconstruction loss
	gval_x = gpu(val_x[:,:,:,1:min(1000, size(val_x,4))]);
	local data_itr
	@suppress begin # suppress the batchsize warning
		data_itr = Flux.Data.DataLoader(sample_tensor(tr_x, epochsize), batchsize=batchsize)
	end
    init_rloss = batched_loss(x->reconstruction_loss(model,x), gval_x, batchsize)
    control_rloss = init_rloss

    # params
    # this is according to their code right
    aeps = Flux.params(model.e, model.d, model.f, model.g)
    cps = Flux.params(model.c)

    # losses
    aeloss(x) = factor_aeloss(model, gpu(x), γ) + λ*sum(l2, aeps)
    closs(x) = disentangle_per_latent ?
    	factor_closs_per_latent(model, gpu(x)) : factor_closs_per_dim(model, gpu(x))

    # rest of the setup
    rdata = []
    hist = MVHistory()
    nl = length(model.e)
    zs = [[] for _ in 1:nl]

	println("Training in progress...")
	for epoch in 1:nepochs
		for x in data_itr
			# first update the autoencoder
			param_update!(aeloss, aeps, x, aeopt)

			# then update the critic/discriminator
			param_update!(closs, cps, x, copt)
		end

		 # sometimes the model is stuck in some local minima and cant get out - better restart the training
        rloss = batched_loss(x->reconstruction_loss(model, x), gval_x, batchsize)
        if epoch == initial_convergence_epochs
            if abs((init_rloss - rloss)/init_rloss) < convergence_threshold
                @info "\nInitial improvement after $(initial_convergence_epochs) epochs is $(init_rloss) => $rloss, model probably stuck in local minima, terminating training. Try restarting the model."
                return nothing, nothing, nothing
            end
        end

		# logging
		Flux.Zygote.ignore() do
			# compute loss values
			ael=round(batched_loss(aeloss, gval_x, batchsize), digits=2)
			cl=round(batched_loss(closs, gval_x, batchsize), digits=4)
			rl=round(batched_loss(x->reconstruction_loss(model, x), gval_x, batchsize), digits=2)
			kldl=round(batched_loss(x->kld_loss(model, x), gval_x, batchsize), digits=2)
			tcl=round(batched_loss(x->tc_loss(model, x), gval_x, batchsize), digits=4)
			println("Epoch $(epoch)/$(nepochs), validation loss: AE=$ael | C=$cl | R=$rl | KLD=$kldl | TC=$tcl")

			push!(hist, :autoencoder_loss, epoch, ael)
			push!(hist, :critic_loss, epoch, cl)
			push!(hist, :reconstruction_loss, epoch, rl)
			push!(hist, :kld, epoch, kldl)
			push!(hist, :total_correlation, epoch, tcl)
			push!(rdata, cpu(reconstruct(model, gval_x)))

			# also, compute the disentanglement metric
			dmls = []
			dmds = []
			if !isnothing(factors)
				zs_all = cpu(encode_all(model, cat(tr_x, val_x, dims=4)))
				y = (isnothing(tr_y) || isnothing(val_y)) ? nothing : vcat(tr_y, val_y)
				# in case factors is just one vector
				if typeof(factors[1]) == Symbol || length(factors[1]) == 1
					factors = [factors]
				end
				for factor_vec in factors
					dml = disentanglement_per_latent(zs_all, y, factor_vec, disentanglement_repeats)
					dmd = disentanglement_per_dim(vcat(zs_all...), y, factor_vec, disentanglement_repeats)
					println("Disentanglement of $(factor_vec): $dml (per latent) | $dmd (per dim)")
					push!(dmls, dml)
					push!(dmds, dmd)
				end
				println("")
			end
			push!(hist, :disentanglement_per_latent, epoch, dmls)
			push!(hist, :disentanglement_per_dim, epoch, dmds)

			#
			for i in 1:nl
				z = encode(model, gval_x, i)
				push!(zs[i], cpu(z))
			end
		end

		# early stopping
		if early_stopping && control_rloss < rloss
			@info "Stoppping early after $epoch epochs due to no improvement."
			return hist, rdata, zs
		else
			control_rloss = rloss
		end
	end

	return hist, rdata, zs
end

#### THIS IS FOR TAB DATA ####
"""
	train_fvlae(zdim, hdim, batchsize, discriminator_hdim, nepochs, tr_x, val_x;
	λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,2), layer_depth=1, lr=0.001f0,
	activation="relu", discriminator_nlayers=3,
	initial_convergence_threshold=0, initial_convergence_epochs=10,
    max_retrain_tries=10)

Train a factored dense VLAE.
"""
function train_fvlae(zdim, hdims, batchsize, d_hdim, nepochs, tr_x::AbstractArray{T,2},
	val_x::AbstractArray{T,2}; λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,2),
	layer_depth=1, lr=0.001f0, activation="relu", discriminator_nlayers=3,
	initial_convergence_threshold=0f0, initial_convergence_epochs=nepochs+1,
    max_retrain_tries=10, verb=true, kwargs...) where T

    # this is to ensure that the model converges to something meaningful
	hist, rdata, model, zs, aeopt, copt = nothing, nothing, nothing, nothing, nothing, nothing

    ntries = 0

    while isnothing(hist) && ntries<=max_retrain_tries
        if ntries > 0
            @info "Restarting training automatically..."
        end

        # construct the model
        model = FVLAE(zdim, hdims, d_hdim, size(tr_x); activation=activation,
        	discriminator_nlayers=discriminator_nlayers,
        	layer_depth=layer_depth)
        aeopt = ADAM(lr)
        copt = ADAM(lr)

        # train it
        hist, rdata, zs = train!(model, nepochs, batchsize, tr_x, val_x, aeopt, copt;
            λ=λ, γ=γ, epochsize = epochsize,
            initial_convergence_threshold=initial_convergence_threshold,
            initial_convergence_epochs=initial_convergence_epochs,
            verb=verb)
        ntries += 1
    end

    # if the training was still unsuccesfull
    if isnothing(hist)
        @info "Gave up the training after $ntries tries."
        return model, hist, rdata, zs, aeopt, copt
    end

    return model, hist, rdata, zs, aeopt, copt
end

"""
	train!(model::FVLAE, nepochs, batchsize, tr_x::AbstractArray{T,2},
		val_x::AbstractArray{T,2}, aeopt, copt;
		λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,2), initial_convergence_threshold=0f0,
    	initial_convergence_epochs=50)
"""
function train!(model::FVLAE, nepochs, batchsize, tr_x::AbstractArray{T,2},
	val_x::AbstractArray{T,2}, aeopt, copt;
	λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,2), initial_convergence_threshold=0f0,
    initial_convergence_epochs=50, verb=true) where T

	# data an initial reconstruction loss
	subtr_x = tr_x[:,1:min(1000, size(tr_x,2))]
	subval_x = val_x[:,1:min(1000, size(val_x,2))]
	local data_itr
	@suppress begin # suppress the batchsize warning
		data_itr = Flux.Data.DataLoader(sample_tensor(tr_x, epochsize), batchsize=batchsize)
    end
    init_rloss = batched_loss(x->reconstruction_loss(model,x), subval_x, batchsize)

    # params
    aeps = Flux.params(model.e, model.d, model.f, model.g)
    cps = Flux.params(model.c)

    # losses
    aeloss(x) = factor_aeloss(model, x, γ) + λ*sum(l2, aeps)
    closs(x) = factor_closs(model, x)

    # rest of the setup
    rdata = []
    hist = MVHistory()
    nl = length(model.e)
    zs = [[] for _ in 1:nl]

    if verb
		println("Training in progress...")
	end
	for epoch in 1:nepochs
		for x in data_itr
			# first update the autoencoder
			param_update!(aeloss, aeps, x, aeopt)

			# then update the critic/discriminator
			param_update!(closs, cps, x, copt)
		end

		 # sometimes the model is stuck in some local minima and cant get out - better restart the training
        if epoch == initial_convergence_epochs
            rloss = batched_loss(x->reconstruction_loss(model, x), subval_x, batchsize)
            if abs((init_rloss - rloss)/init_rloss) < initial_convergence_threshold
                @info "Initial improvement after $(initial_convergence_epochs) epochs is $(init_rloss) => $rloss, model probably stuck in local minima, terminating training. Try restarting the model."
                return nothing, nothing, nothing
            end
        end

		# logging
		Flux.Zygote.ignore() do
			ael=round(batched_loss(aeloss, subval_x, batchsize), digits=4)
			cl=round(batched_loss(closs, subval_x, batchsize), digits=4)
			rl=round(batched_loss(x->reconstruction_loss(model, x), subval_x, batchsize), digits=4)
			kldl=round(batched_loss(x->kld_loss(model, x), subval_x, batchsize), digits=2)
			tcl=round(batched_loss(x->tc_loss(model, x), subval_x, batchsize), digits=4)

			if verb
				println("Epoch $(epoch)/$(nepochs), validation loss: AE=$ael | C=$cl | R=$rl | KLD=$kldl | TC=$tcl")
			end
			for i in 1:nl
				z = encode(model, subval_x, i)
				push!(zs[i], z)
			end
			push!(hist, :autoencoder_loss, epoch, ael)
			push!(hist, :critic_loss, epoch, cl)
			push!(hist, :reconstruction_loss, epoch, rl)
			push!(hist, :kld, epoch, kldl)
			push!(hist, :total_correlation, epoch, tcl)
			push!(rdata, reconstruct(model, subval_x))
		end
	end

	return hist, rdata, zs
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
#!/usr/bin/julia

#cont-subtraction.jl
#continuous testing of julia subtraction directives
include("../../unum.jl")
using Unums

error = 0
total = 0

while true
  x1 = exp(100 * randn());  x2 = x1 + exp(100 * randn());

  u1 = convert(Unum{4,6}, x1); u2 = convert(Unum{4,6}, x2)
  u3 = Unums.__diff_exact(u2, -u1, Unums.decode_exp(u2), Unums.decode_exp(u1))
  x3 = convert(Float64, u3)

  if (abs(x3 - (x2 - x1)) > (abs(x3) * 0.000000001))
    println("$x3 != $(x2 - x1) = $x2 - $x1, diff: $(x3 - (x2 - x1))")
    println("u1:", u1)
    println("u2:", u2)
    println("u3:", u3)
    println("u1:  $(bits(u1, " "))")
    println("u2:  $(bits(u2, " "))")
    println("u3:  $(bits(u3, " "))")
    error += 1
  end
  total += 1
  #println("error rate: $(error/total)")
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
const edge_table = (0,265,515,778,1030,1295,1541,1804,2060,2309,2575,2822,3082,
                    3331,3593,3840,400,153,915,666,1430,1183,1941,1692,2460,2197,
                    2975,2710,3482,3219,3993,3728,560,825,51,314,1590,1855,1077,
                    1340,2620,2869,2111,2358,3642,3891,3129,3376,928,681,419,170,
                    1958,1711,1445,1196,2988,2725,2479,2214,4010,3747,3497,3232,
                    1120,1385,1635,1898,102,367,613,876,3180,3429,3695,3942,2154,
                    2403,2665,2912,1520,1273,2035,1786,502,255,1013,764,3580,3317,
                    4095,3830,2554,2291,3065,2800,1616,1881,1107,1370,598,863,85,
                    348,3676,3925,3167,3414,2650,2899,2137,2384,1984,1737,1475,
                    1226,966,719,453,204,4044,3781,3535,3270,3018,2755,2505,2240,
                    2240,2505,2755,3018,3270,3535,3781,4044,204,453,719,966,1226,
                    1475,1737,1984,2384,2137,2899,2650,3414,3167,3925,3676,348,
                    85,863,598,1370,1107,1881,1616,2800,3065,2291,2554,3830,
                    4095,3317,3580,764,1013,255,502,1786,2035,1273,1520,2912,
                    2665,2403,2154,3942,3695,3429,3180,876,613,367,102,1898,
                    1635,1385,1120,3232,3497,3747,4010,2214,2479,2725,2988,1196,
                    1445,1711,1958,170,419,681,928,3376,3129,3891,3642,2358,2111,
                    2869,2620,1340,1077,1855,1590,314,51,825,560,3728,3993,3219,
                    3482,2710,2975,2197,2460,1692,1941,1183,1430,666,915,153,400,
                    3840,3593,3331,3082,2822,2575,2309,2060,1804,1541,1295,1030,
                    778,515,265,0)

const tri_table = ((-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,9,4,10,9,2,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,3,11,1,3,10,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,9,4,3,11,9,11,10,9,-1,-1,-1,-1,-1,-1),
                    (4,12,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,12,3,9,12,1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,10,1,3,4,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,12,3,2,10,12,10,9,12,-1,-1,-1,-1,-1,-1),
                    (4,11,2,12,11,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,11,2,1,9,11,9,12,11,-1,-1,-1,-1,-1,-1),
                    (4,10,1,4,12,10,12,11,10,-1,-1,-1,-1,-1,-1),
                    (10,9,11,11,9,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,8,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,4,1,8,4,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,9,5,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,2,10,5,8,2,8,4,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,9,5,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,5,8,4,1,5,2,3,11,-1,-1,-1,-1,-1,-1),
                    (10,3,11,10,1,3,9,5,8,-1,-1,-1,-1,-1,-1),
                    (3,11,10,3,10,8,3,8,4,8,10,5,-1,-1,-1),
                    (9,5,8,4,12,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,5,8,12,3,5,3,1,5,-1,-1,-1,-1,-1,-1),
                    (10,1,2,9,5,8,3,4,12,-1,-1,-1,-1,-1,-1),
                    (5,8,12,10,5,12,10,12,3,10,3,2,-1,-1,-1),
                    (4,11,2,4,12,11,8,9,5,-1,-1,-1,-1,-1,-1),
                    (2,12,11,2,5,12,2,1,5,8,12,5,-1,-1,-1),
                    (5,8,9,10,1,12,10,12,11,12,1,4,-1,-1,-1),
                    (5,8,12,5,12,10,10,12,11,-1,-1,-1,-1,-1,-1),
                    (10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,6,5,1,9,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,6,5,2,6,1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,6,5,9,4,6,4,2,6,-1,-1,-1,-1,-1,-1),
                    (2,3,11,10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,9,2,3,11,5,10,6,-1,-1,-1,-1,-1,-1),
                    (6,3,11,6,5,3,5,1,3,-1,-1,-1,-1,-1,-1),
                    (3,11,6,4,3,6,4,6,5,4,5,9,-1,-1,-1),
                    (10,6,5,3,4,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,12,3,1,9,12,5,10,6,-1,-1,-1,-1,-1,-1),
                    (1,6,5,1,2,6,3,4,12,-1,-1,-1,-1,-1,-1),
                    (3,2,6,3,6,9,3,9,12,5,9,6,-1,-1,-1),
                    (11,4,12,11,2,4,10,6,5,-1,-1,-1,-1,-1,-1),
                    (5,10,6,1,9,2,9,11,2,9,12,11,-1,-1,-1),
                    (6,5,1,6,1,12,6,12,11,12,1,4,-1,-1,-1),
                    (6,5,9,6,9,11,11,9,12,-1,-1,-1,-1,-1,-1),
                    (10,8,9,6,8,10,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,4,1,10,6,4,6,8,4,-1,-1,-1,-1,-1,-1),
                    (1,8,9,1,2,8,2,6,8,-1,-1,-1,-1,-1,-1),
                    (2,6,4,4,6,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,8,9,10,6,8,11,2,3,-1,-1,-1,-1,-1,-1),
                    (11,2,3,10,6,1,6,4,1,6,8,4,-1,-1,-1),
                    (9,1,3,9,3,6,9,6,8,11,6,3,-1,-1,-1),
                    (3,11,6,3,6,4,4,6,8,-1,-1,-1,-1,-1,-1),
                    (8,10,6,8,9,10,4,12,3,-1,-1,-1,-1,-1,-1),
                    (10,6,8,10,8,3,10,3,1,3,8,12,-1,-1,-1),
                    (3,4,12,1,2,9,2,8,9,2,6,8,-1,-1,-1),
                    (12,3,2,12,2,8,8,2,6,-1,-1,-1,-1,-1,-1),
                    (10,6,9,9,6,8,11,2,4,11,4,12,-1,-1,-1),
                    (6,8,1,6,1,10,8,12,1,2,1,11,12,11,1),
                    (12,11,1,12,1,4,11,6,1,9,1,8,6,8,1),
                    (12,11,6,8,12,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,6,11,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,1,2,6,11,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,9,4,2,10,9,6,11,7,-1,-1,-1,-1,-1,-1),
                    (2,7,6,3,7,2,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,7,6,2,3,7,4,1,9,-1,-1,-1,-1,-1,-1),
                    (10,7,6,10,1,7,1,3,7,-1,-1,-1,-1,-1,-1),
                    (6,10,9,6,9,3,6,3,7,4,3,9,-1,-1,-1),
                    (3,4,12,11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,1,9,12,3,1,11,7,6,-1,-1,-1,-1,-1,-1),
                    (1,2,10,3,4,12,6,11,7,-1,-1,-1,-1,-1,-1),
                    (6,11,7,2,10,3,10,12,3,10,9,12,-1,-1,-1),
                    (7,4,12,7,6,4,6,2,4,-1,-1,-1,-1,-1,-1),
                    (1,9,12,1,12,6,1,6,2,6,12,7,-1,-1,-1),
                    (4,12,7,1,4,7,1,7,6,1,6,10,-1,-1,-1),
                    (7,6,10,7,10,12,12,10,9,-1,-1,-1,-1,-1,-1),
                    (6,11,7,5,8,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,4,1,5,8,4,7,6,11,-1,-1,-1,-1,-1,-1),
                    (2,10,1,6,11,7,9,5,8,-1,-1,-1,-1,-1,-1),
                    (11,7,6,2,10,8,2,8,4,8,10,5,-1,-1,-1),
                    (7,2,3,7,6,2,5,8,9,-1,-1,-1,-1,-1,-1),
                    (2,3,6,6,3,7,4,1,5,4,5,8,-1,-1,-1),
                    (9,5,8,10,1,6,1,7,6,1,3,7,-1,-1,-1),
                    (8,4,10,8,10,5,4,3,10,6,10,7,3,7,10),
                    (4,12,3,8,9,5,11,7,6,-1,-1,-1,-1,-1,-1),
                    (6,11,7,5,8,3,5,3,1,3,8,12,-1,-1,-1),
                    (1,2,10,5,8,9,3,4,12,6,11,7,-1,-1,-1),
                    (10,3,2,10,12,3,10,5,12,8,12,5,6,11,7),
                    (9,5,8,4,12,6,4,6,2,6,12,7,-1,-1,-1),
                    (6,2,12,6,12,7,2,1,12,8,12,5,1,5,12),
                    (1,6,10,1,7,6,1,4,7,12,7,4,9,5,8),
                    (7,6,10,7,10,12,5,8,10,8,12,10,-1,-1,-1),
                    (11,5,10,7,5,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,11,7,5,10,11,1,9,4,-1,-1,-1,-1,-1,-1),
                    (11,1,2,11,7,1,7,5,1,-1,-1,-1,-1,-1,-1),
                    (9,4,2,9,2,7,9,7,5,7,2,11,-1,-1,-1),
                    (2,5,10,2,3,5,3,7,5,-1,-1,-1,-1,-1,-1),
                    (4,1,9,2,3,10,3,5,10,3,7,5,-1,-1,-1),
                    (1,3,5,5,3,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,4,3,9,3,5,5,3,7,-1,-1,-1,-1,-1,-1),
                    (11,5,10,11,7,5,12,3,4,-1,-1,-1,-1,-1,-1),
                    (1,9,3,3,9,12,5,10,11,5,11,7,-1,-1,-1),
                    (4,12,3,1,2,7,1,7,5,7,2,11,-1,-1,-1),
                    (7,5,2,7,2,11,5,9,2,3,2,12,9,12,2),
                    (10,7,5,10,4,7,10,2,4,12,7,4,-1,-1,-1),
                    (9,12,2,9,2,1,12,7,2,10,2,5,7,5,2),
                    (4,12,7,4,7,1,1,7,5,-1,-1,-1,-1,-1,-1),
                    (7,5,9,12,7,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,11,7,8,9,11,9,10,11,-1,-1,-1,-1,-1,-1),
                    (1,8,4,1,11,8,1,10,11,7,8,11,-1,-1,-1),
                    (11,7,8,2,11,8,2,8,9,2,9,1,-1,-1,-1),
                    (11,7,8,11,8,2,2,8,4,-1,-1,-1,-1,-1,-1),
                    (2,3,7,2,7,9,2,9,10,9,7,8,-1,-1,-1),
                    (3,7,10,3,10,2,7,8,10,1,10,4,8,4,10),
                    (8,9,1,8,1,7,7,1,3,-1,-1,-1,-1,-1,-1),
                    (8,4,3,7,8,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,12,11,7,9,11,9,10,9,7,8,-1,-1,-1),
                    (3,1,8,3,8,12,1,10,8,7,8,11,10,11,8),
                    (2,9,1,2,8,9,2,11,8,7,8,11,3,4,12),
                    (12,3,2,12,2,8,11,7,2,7,8,2,-1,-1,-1),
                    (9,10,7,9,7,8,10,2,7,12,7,4,2,4,7),
                    (1,10,2,12,7,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,9,1,8,1,7,4,12,1,12,7,1,-1,-1,-1),
                    (8,12,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,7,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,9,12,8,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,2,10,12,8,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,2,10,9,4,2,12,8,7,-1,-1,-1,-1,-1,-1),
                    (11,2,3,7,12,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,11,4,1,9,7,12,8,-1,-1,-1,-1,-1,-1),
                    (3,10,1,3,11,10,7,12,8,-1,-1,-1,-1,-1,-1),
                    (7,12,8,3,11,4,11,9,4,11,10,9,-1,-1,-1),
                    (8,3,4,7,3,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (8,1,9,8,7,1,7,3,1,-1,-1,-1,-1,-1,-1),
                    (3,8,7,3,4,8,1,2,10,-1,-1,-1,-1,-1,-1),
                    (2,7,3,2,9,7,2,10,9,9,8,7,-1,-1,-1),
                    (11,8,7,11,2,8,2,4,8,-1,-1,-1,-1,-1,-1),
                    (11,8,7,2,8,11,2,9,8,2,1,9,-1,-1,-1),
                    (1,4,8,1,8,11,1,11,10,7,11,8,-1,-1,-1),
                    (8,7,11,8,11,9,9,11,10,-1,-1,-1,-1,-1,-1),
                    (7,9,5,12,9,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,7,12,4,1,7,1,5,7,-1,-1,-1,-1,-1,-1),
                    (9,7,12,9,5,7,10,1,2,-1,-1,-1,-1,-1,-1),
                    (10,5,7,10,7,4,10,4,2,12,4,7,-1,-1,-1),
                    (7,9,5,7,12,9,3,11,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,4,1,12,1,7,12,1,5,7,-1,-1,-1),
                    (5,12,9,5,7,12,1,3,10,3,11,10,-1,-1,-1),
                    (11,10,4,11,4,3,10,5,4,12,4,7,5,7,4),
                    (9,3,4,9,5,3,5,7,3,-1,-1,-1,-1,-1,-1),
                    (1,5,3,5,7,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,10,1,3,4,5,3,5,7,5,4,9,-1,-1,-1),
                    (2,10,5,2,5,3,3,5,7,-1,-1,-1,-1,-1,-1),
                    (9,2,4,9,7,2,9,5,7,7,11,2,-1,-1,-1),
                    (11,2,1,11,1,7,7,1,5,-1,-1,-1,-1,-1,-1),
                    (5,7,4,5,4,9,7,11,4,1,4,10,11,10,4),
                    (11,10,5,7,11,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,10,6,8,7,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,4,5,10,6,12,8,7,-1,-1,-1,-1,-1,-1),
                    (6,1,2,6,5,1,8,7,12,-1,-1,-1,-1,-1,-1),
                    (12,8,7,9,4,5,4,6,5,4,2,6,-1,-1,-1),
                    (10,6,5,11,2,3,8,7,12,-1,-1,-1,-1,-1,-1),
                    (7,12,8,2,3,11,1,9,4,5,10,6,-1,-1,-1),
                    (8,7,12,6,5,11,5,3,11,5,1,3,-1,-1,-1),
                    (4,5,9,4,6,5,4,3,6,11,6,3,12,8,7),
                    (8,3,4,8,7,3,6,5,10,-1,-1,-1,-1,-1,-1),
                    (10,6,5,1,9,7,1,7,3,7,9,8,-1,-1,-1),
                    (4,7,3,4,8,7,2,6,1,6,5,1,-1,-1,-1),
                    (7,3,9,7,9,8,3,2,9,5,9,6,2,6,9),
                    (10,6,5,11,2,7,2,8,7,2,4,8,-1,-1,-1),
                    (2,7,11,2,8,7,2,1,8,9,8,1,10,6,5),
                    (5,1,11,5,11,6,1,4,11,7,11,8,4,8,11),
                    (8,7,11,8,11,9,6,5,11,5,9,11,-1,-1,-1),
                    (7,10,6,7,12,10,12,9,10,-1,-1,-1,-1,-1,-1),
                    (4,7,12,1,7,4,1,6,7,1,10,6,-1,-1,-1),
                    (1,12,9,1,6,12,1,2,6,6,7,12,-1,-1,-1),
                    (7,12,4,7,4,6,6,4,2,-1,-1,-1,-1,-1,-1),
                    (2,3,11,10,6,12,10,12,9,12,6,7,-1,-1,-1),
                    (1,12,4,1,7,12,1,10,7,6,7,10,2,3,11),
                    (12,9,6,12,6,7,9,1,6,11,6,3,1,3,6),
                    (7,12,4,7,4,6,3,11,4,11,6,4,-1,-1,-1),
                    (6,9,10,6,3,9,6,7,3,4,9,3,-1,-1,-1),
                    (10,6,7,10,7,1,1,7,3,-1,-1,-1,-1,-1,-1),
                    (2,6,9,2,9,1,6,7,9,4,9,3,7,3,9),
                    (2,6,7,3,2,7,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,4,7,2,7,11,4,9,7,6,7,10,9,10,7),
                    (11,2,1,11,1,7,10,6,1,6,7,1,-1,-1,-1),
                    (1,4,9,6,7,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (11,6,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,6,11,8,6,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (12,6,11,12,8,6,9,4,1,-1,-1,-1,-1,-1,-1),
                    (6,12,8,6,11,12,2,10,1,-1,-1,-1,-1,-1,-1),
                    (11,8,6,11,12,8,10,9,2,9,4,2,-1,-1,-1),
                    (12,2,3,12,8,2,8,6,2,-1,-1,-1,-1,-1,-1),
                    (1,9,4,2,3,8,2,8,6,8,3,12,-1,-1,-1),
                    (10,8,6,10,3,8,10,1,3,3,12,8,-1,-1,-1),
                    (8,6,3,8,3,12,6,10,3,4,3,9,10,9,3),
                    (3,6,11,3,4,6,4,8,6,-1,-1,-1,-1,-1,-1),
                    (9,3,1,9,6,3,9,8,6,11,3,6,-1,-1,-1),
                    (10,1,2,6,11,4,6,4,8,4,11,3,-1,-1,-1),
                    (10,9,3,10,3,2,9,8,3,11,3,6,8,6,3),
                    (2,4,6,4,8,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,9,8,1,8,2,2,8,6,-1,-1,-1,-1,-1,-1),
                    (10,1,4,10,4,6,6,4,8,-1,-1,-1,-1,-1,-1),
                    (10,9,8,6,10,8,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (6,9,5,6,11,9,11,12,9,-1,-1,-1,-1,-1,-1),
                    (6,1,5,6,12,1,6,11,12,12,4,1,-1,-1,-1),
                    (1,2,10,9,5,11,9,11,12,11,5,6,-1,-1,-1),
                    (11,12,5,11,5,6,12,4,5,10,5,2,4,2,5),
                    (3,6,2,3,9,6,3,12,9,5,6,9,-1,-1,-1),
                    (1,5,12,1,12,4,5,6,12,3,12,2,6,2,12),
                    (1,3,6,1,6,10,3,12,6,5,6,9,12,9,6),
                    (10,5,6,3,12,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,6,11,4,6,3,4,5,6,4,9,5,-1,-1,-1),
                    (6,11,3,6,3,5,5,3,1,-1,-1,-1,-1,-1,-1),
                    (4,11,3,4,6,11,4,9,6,5,6,9,1,2,10),
                    (6,11,3,6,3,5,2,10,3,10,5,3,-1,-1,-1),
                    (9,5,6,9,6,4,4,6,2,-1,-1,-1,-1,-1,-1),
                    (1,5,6,2,1,6,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (9,5,6,9,6,4,10,1,6,1,4,6,-1,-1,-1),
                    (10,5,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,12,8,5,10,12,10,11,12,-1,-1,-1,-1,-1,-1),
                    (1,9,4,5,10,8,10,12,8,10,11,12,-1,-1,-1),
                    (2,11,12,2,12,5,2,5,1,8,5,12,-1,-1,-1),
                    (4,2,5,4,5,9,2,11,5,8,5,12,11,12,5),
                    (5,12,8,10,12,5,10,3,12,10,2,3,-1,-1,-1),
                    (10,8,5,10,12,8,10,2,12,3,12,2,1,9,4),
                    (12,8,5,12,5,3,3,5,1,-1,-1,-1,-1,-1,-1),
                    (12,8,5,12,5,3,9,4,5,4,3,5,-1,-1,-1),
                    (3,10,11,3,8,10,3,4,8,8,5,10,-1,-1,-1),
                    (10,11,8,10,8,5,11,3,8,9,8,1,3,1,8),
                    (4,8,11,4,11,3,8,5,11,2,11,1,5,1,11),
                    (2,11,3,9,8,5,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,10,2,5,2,8,8,2,4,-1,-1,-1,-1,-1,-1),
                    (5,10,2,5,2,8,1,9,2,9,8,2,-1,-1,-1),
                    (5,1,4,8,5,4,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (5,9,8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (10,11,9,11,12,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,1,10,4,10,12,12,10,11,-1,-1,-1,-1,-1,-1),
                    (1,2,11,1,11,9,9,11,12,-1,-1,-1,-1,-1,-1),
                    (4,2,11,12,4,11,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,3,12,2,12,10,10,12,9,-1,-1,-1,-1,-1,-1),
                    (4,1,10,4,10,12,2,3,10,3,12,10,-1,-1,-1),
                    (1,3,12,9,1,12,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (4,3,12,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,9,3,9,11,11,9,10,-1,-1,-1,-1,-1,-1),
                    (10,11,3,1,10,3,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (3,4,9,3,9,11,1,2,9,2,11,9,-1,-1,-1),
                    (2,11,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (2,4,9,10,2,9,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,10,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (1,4,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1),
                    (-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1))

"""
`marching_cubes(sdf::SignedDistanceField, [iso = 0.0,] [MT = HomogenousMesh{Point{3,Float64},Face{3,Int,0}}])`

Construct a `HomogenousMesh` from a `SignedDistanceField` using the
marching cubes algorithm. This method is faster than Marching Tetrahedra
and generates fewer vertices and faces (about 1/4 as many).
However it may generate non-manifold meshes, while Marching
Tetrahedra guarentees a manifold mesh.
"""
function marching_cubes{ST,FT,M<:AbstractMesh}(sdf::SignedDistanceField{3,ST,FT},
                               iso=0.0,
                               MT::Type{M}=HomogenousMesh)
    nx, ny, nz = size(sdf)
    h = HyperRectangle(sdf)
    w = widths(h)
    s = Point{3,Float64}(w[1]/nx, w[2]/ny, w[3]/nz)

    # arrays for vertices and faces
    vts = Point{3,Float64}[]
    fcs = Face{3,Int,0}[]
    vertlist = Vector{Point{3,Float64}}(12)
    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = sdf[xi,yi,zi] < iso ? 1 : 0
        sdf[xi+1,yi,zi] < iso && (cubeindex |= 2)
        sdf[xi+1,yi+1,zi] < iso && (cubeindex |= 4)
        sdf[xi,yi+1,zi] < iso && (cubeindex |= 8)
        sdf[xi,yi,zi+1] < iso && (cubeindex |= 16)
        sdf[xi+1,yi,zi+1] < iso && (cubeindex |= 32)
        sdf[xi+1,yi+1,zi+1] < iso && (cubeindex |= 64)
        sdf[xi,yi+1,zi+1] < iso && (cubeindex |= 128)

        # Cube is entirely in/out of the surface
        edge_table[cubeindex+1] == 0 && continue

        points = (Point{3,Float64}(xi-1,yi-1,zi-1).*s,
                  Point{3,Float64}(xi,yi-1,zi-1).*s,
                  Point{3,Float64}(xi,yi,zi-1).*s,
                  Point{3,Float64}(xi-1,yi,zi-1).*s,
                  Point{3,Float64}(xi-1,yi-1,zi).*s,
                  Point{3,Float64}(xi,yi-1,zi).*s,
                  Point{3,Float64}(xi,yi,zi).*s,
                  Point{3,Float64}(xi-1,yi,zi).*s)

        # Find the vertices where the surface intersects the cube
        if (edge_table[cubeindex+1] & 1 != 0)
          vertlist[1] =
             vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi])
        end
        if (edge_table[cubeindex+1] & 2 != 0)
          vertlist[2] =
             vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi])
        end
        if (edge_table[cubeindex+1] & 4 != 0)
          vertlist[3] =
             vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi])
        end
        if (edge_table[cubeindex+1] & 8 != 0)
          vertlist[4] =
             vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi])
        end
        if (edge_table[cubeindex+1] & 16 != 0)
          vertlist[5] =
             vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1])
        end
        if (edge_table[cubeindex+1] & 32 != 0)
          vertlist[6] =
             vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1])
        end
        if (edge_table[cubeindex+1] & 64 != 0)
          vertlist[7] =
             vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1])
        end
        if (edge_table[cubeindex+1] & 128 != 0)
          vertlist[8] =
             vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1])
        end
        if (edge_table[cubeindex+1] & 256 != 0)
          vertlist[9] =
             vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1])
        end
        if (edge_table[cubeindex+1] & 512 != 0)
          vertlist[10] =
             vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1])
        end
        if (edge_table[cubeindex+1] & 1024 != 0)
          vertlist[11] =
             vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1])
        end
        if (edge_table[cubeindex+1] & 2048 != 0)
          vertlist[12] =
             vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1])
        end

        # Create the triangle
        for i = 1:3:13
            tri_table[cubeindex+1][i] == -1 && break
            push!(vts, vertlist[tri_table[cubeindex+1][i  ]])
            push!(vts, vertlist[tri_table[cubeindex+1][i+1]])
            push!(vts, vertlist[tri_table[cubeindex+1][i+2]])
            fct = length(vts)
            push!(fcs, Face{3,Int,0}(fct, fct-1, fct-2))
        end
    end
    MT(vts,fcs)
end

# Linearly interpolate the position where an isosurface cuts
# an edge between two vertices, each with their own scalar value
function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)

    abs(iso - valp1) < eps && return p1
    abs(iso - valp2) < eps && return p2
    abs(valp1-valp2) < eps && return p1
    mu = (iso - valp1) / (valp2 - valp1)
    p = p1 + mu * (p2 - p1)

    return p
end


text:
using Pkg
Pkg.activate(".")

filenames = (train_image="train-images-idx3-ubyte",
             train_label="train-labels-idx1-ubyte",
             test_image="t10k-images-idx3-ubyte",
             test_label="t10k-labels-idx1-ubyte")

url = "http://yann.lecun.com/exdb/mnist/"

using CodecZlib

for fname in filenames
    fnamegz = fname * ".gz"
    isfile(fname) && continue
    @info "Processing $fname"
    if !isfile("$fname.gz")
        @info "Downloading $fnamegz"
        download(url*fnamegz, fnamegz)
    end
    open(fname, "w") do io_out
        open(fnamegz) do io_in
            data = read(GzipDecompressorStream(io_in))
            write(io_out, data)
        end
    end
end
@info "Done!"


text:
using SpectralDistances, Distances, DSP

@testset "Sinkhorn convolutional" begin
    @info "Testing Sinkhorn convolutional"
    m,n,T = 8, 10, Float64
    for m = 8:2:12, n = 8:2:12, T = (Float64, )
        # @show m,n,T

        a = zeros(T, m, n)
        a[2,2] = 1
        b = zeros(T, m, n)
        b[3,3] = 1
        d1 = sinkhorn_convolutional(a,b,τ=1e2)[1]
        b = zeros(T, m, n)
        b[4,4] = 1
        d2 = sinkhorn_convolutional(a,b,τ=1e2)[1]
        @test 2sqrt(d1) ≈ sqrt(d2) rtol = 0.1
        b = zeros(T, m, n)
        b[6,6] = 1
        d3 = sinkhorn_convolutional(a,b,τ=1e2)[1]
        @test 2sqrt(d2) ≈ sqrt(d3) rtol = 0.1

        # Do the same as above but with different log stabilization
        a = zeros(T, m, n)
        a[2,2] = 1
        b = zeros(T, m, n)
        b[3,3] = 1
        d4 = sinkhorn_convolutional(a,b,τ=1e5)[1]
        b = zeros(T, m, n)
        b[4,4] = 1
        d5 = sinkhorn_convolutional(a,b,τ=1e5)[1]
        @test 2sqrt(d4) ≈ sqrt(d5) rtol = 0.1
        b = zeros(T, m, n)
        b[6,6] = 1
        d6 = sinkhorn_convolutional(a,b,τ=1e5)[1]
        @test 2sqrt(d5) ≈ sqrt(d6) rtol = 0.1

        @test d1 ≈ d4 # Test that the result is the same with and without log stabilization.
        @test d2 ≈ d5
        @test d3 ≈ d6
    end
end

@testset "Invariant axis" begin
    @info "Testing Invariant axis"
    w,h = 5,5
    for w = [5, 6, 20], h = [5,6,20]
        @show w,h
        β = 0.05
        D = [abs2((i-j)/(h-1)) for i = 1:h, j = 1:h]
        B = [zeros(h,w) for _ in 1:2]

        B[1][1,1] = 1
        B[1][1,5] = 1
        B[2][2,1] = 1
        B[2][2,5] = 1

        B = s1.(B)
        c,V,U = sinkhorn_convolutional(B[1], B[2], β=β, τ=1e20)
        c1,V1,U1 = sinkhorn_convolutional(B[1], B[1], β=β, τ=1e20)
        c2,V2,U2 = sinkhorn_convolutional(B[2], B[2], β=β, τ=1e20)
        c3 = sqrt((c - 0.5(c1+c2))) # The size of the summed over dimension needs to multiply here
        v,u = sum(B[1], dims=2)[:], sum(B[2], dims=2)[:] # Integrate over the dimension to be sensitive to (because we need to be invariant to this dimension in this step)
        @assert sum(u) ≈ 1
        @assert sum(v) ≈ 1
        Γ = discrete_grid_transportplan(v, u)    # Solve 1D OT between integrated measures
        invariant_cost = sqrt(dot(Γ, D))

        @test c3 ≈ invariant_cost rtol=0.15

        if isinteractive()
            Plots.heatmap(Γ, layout=5, title="Transport plan")
            Plots.heatmap!(V, sp=2, title="V")
            Plots.heatmap!(U, sp=3, title="U")
            Plots.plot!(v, sp=4, lab=invariant_cost)
            Plots.plot!(u, sp=5, lab=c3)
        end

        invariant_cost = dot(Γ, D)
        di = ConvOptimalTransportDistance(β=β, invariant_axis=1)
        d = ConvOptimalTransportDistance(β=β, invariant_axis=0)
        cdist = d(B[1], B[2]) - 0.5*(d(B[1], B[1]) + d(B[2], B[2]))
        @test cdist ≈ c3^2
        ci = (di(B[1], B[2]) - 0.5*(di(B[1], B[1]) + di(B[2], B[2])))

        @test ci ≈ cdist-invariant_cost rtol=0.15 atol=sqrt(eps())

        @test di(B[1], B[2]) ≈ d(B[1], B[2])-invariant_cost rtol=0.15

    end
end

@testset "Misc convolutional" begin
    @info "Testing Misc convolutional"

    @test -0.72 < SpectralDistances.default_dynamic_floor(rand(100,100)) < -0.6
    a1 = zeros(Float64, 10,10)
    a1[2,2] = 1
    a2 = zeros(Float64, 10,10)
    a2[6,6] = 1
    A = [a1,a2]
    β = 0.001
    λ = [0.5, 0.5]
    b = barycenter_convolutional(A,λ,β=β)
    @test maximum(b) == b[4,4]
    @test sum(b) ≈ 1

    d = ConvOptimalTransportDistance(β=0.01)
    l = barycentric_coordinates(d, A, a1)[1]
    @test l ≈ [1, 0]
    @test sum(l) ≈ 1
    l = barycentric_coordinates(d, A, a2)[1]
    @test l ≈ [0, 1]

    l = barycentric_coordinates(d, A, b)[1]
    @test l ≈ [0.5, 0.5]

    A1 = spectrogram(sin.(1:10000) .+ 0.1randn(10000), 256, window=hanning)
    A2 = spectrogram(sin.(1:10000) .+ 0.1randn(10000), 256, window=hanning)
    A = [A1, A2]

    @test -5 < SpectralDistances.default_dynamic_floor(A1) < -4
    @test -5 < SpectralDistances.default_dynamic_floor(A) < -4

    o = similar(A1.power)
    @test normalize_spectrogram!(o,A1) == normalize_spectrogram(A1)
    A1n = normalize_spectrogram(A1)
    A1f = mask_filter(A1n)
    @test sum(A1f) < sum(A1n)

    B = barycenter_convolutional(A,λ,β=β, dynamic_floor=-3)
    @test B isa typeof(A1)

    d = ConvOptimalTransportDistance(β=β, dynamic_floor=-3.0)
    B2 = barycenter(d, A)
    @test power(B) == power(B2)
    @test time(B) == time(B2)

    @test d(A1,A2) > d(A1,A1)
    @test d(A1,A2) > d(A2,A2)
    @test d(A1,A2) - 0.5*(d(A1,A1) + d(A2,A2)) > 0

    A3 = spectrogram(sin.(LinRange(0.8,1.2,300_000) .* (1:300_000)) .+ 0.01randn(300_000), 256, window=hanning)

    d = ConvOptimalTransportDistance(β=0.05, dynamic_floor=-3.0)
    @time D = distance_profile(d, A1, A3, tol=1e-6, stride=15)
    isinteractive() && plot(D)
    @test 40 <= argmin(D) <= 70

    dslide = SlidingConvOptimalTransportDistance(d)
    @test dslide(A1, A3; tol=1e-6, stride=15) ≈ minimum(D)
    @test SlidingConvOptimalTransportDistance(β=0.05, dynamic_floor=-3.0).d.β == SlidingConvOptimalTransportDistance(ConvOptimalTransportDistance(β=0.05, dynamic_floor=-3.0)).d.β

end

# @test B.power == barycenter_convolutional([A1.power, A2.power],λ,β=β)

# w = SpectralDistances.BCWorkspace(A, β)
# @btime SpectralDistances.barycenter_convolutional($w,$A,$λ)
# 1.198 ms (46 allocations: 1.26 MiB)
# 585.176 μs (46 allocations: 1.26 MiB) # add @avx
# 469.132 μs (16 allocations: 313.13 KiB) # operate inplace
# 387.881 μs (0 allocations: 0 bytes) # views and generator in sum

##
# m,n = 1280,1280
# x = zeros(m,n)
# x[10,10] = 1
# β = 0.0003 # The threshold for imfilter to be faster is around 0.0003 on size 128
#
# t = LinRange(0, 1, m)
# Y, X = meshgrid(t, t)
# xi1 = @. exp(-(X - Y)^2 / β)
#
# t = LinRange(0, 1, n)
# Y, X = meshgrid(t, t)
# xi2 = @. exp(-(X - Y)^2 / β)
#
# K1(x) = xi1 * x * xi2
#
#
# using ImageFiltering
#
# t = LinRange(-0.5, 0.5, m)
# k1 = @. exp(-t^2 / β)
# t = LinRange(-0.5, 0.5, n)
# k2 = @. exp(-t^2 / β)
#
# function kfg(x)
#     l = round(Int,5x)
#     l = isodd(l) ? l : l+1
#     KernelFactors.gaussian(x,l)
# end
#
# kk1 = kfg(sqrt(0.5β)*m)
# kk2 = kfg(sqrt(0.5β)*n)
# kk1,kk2 = kk1 ./ maximum(kk1), kk2 ./ maximum(kk2)
# # kernf = kernelfactors((centered(k1), centered(k2)))
#
# # kernf = KernelFactors.gaussian((sqrt(β)*m, sqrt(β)*m))
# kernf = kernelfactors((kk1, kk2))
# kernp = broadcast(*, kernf...)
#
# CPU1(Algorithm.FIR())
# K2(x) = imfilter(x, kernf)
#
# plot(heatmap(K1(x)), heatmap(K2(x)), heatmap(K2(x)-K1(x)))
#
# using FFTW
# FFTW.set_num_threads(2)
#
#
# using ImageFiltering.ComputationalResources
# @btime $xi1 * $x * $xi2
# @btime imfilter!($(similar(x)), $x, $kernf)
#
# using LoopVectorization
# function filter2davx!(out::AbstractMatrix, A::AbstractMatrix, kerns)
#     for kern in kerns
#         @avx for J in CartesianIndices(out)
#             tmp = zero(eltype(out))
#             for I ∈ CartesianIndices(kern)
#                 tmp += A[I + J] * kern[I]
#             end
#             out[J] = tmp
#         end
#     end
#     out
# end
#
# out = similar(x)
#
# out = ImageFiltering.OffsetArray(similar(x, size(x).-2), 1, 1)
#
# ka = convert(AbstractArray, kernf[1])
# kb = convert(AbstractArray, kernf[2])
#
# @btime filter2davx!($out, $x, $((ka, kb)));
#
# plot(heatmap(K1(x)), heatmap(K2(x)), heatmap(K2(x)-K1(x)), heatmap(out-K1(x)))


text:
i = 0
while true
    global i += 1
    x, y = parse.(Int, split(readline()))
    (x == y == 0) && break
    if x < y
        println(x, " ", y)
    else
        println(y, " ", x)
    end
end


text:
# "Hall of Shame" time matrices that broke previous tests
function breaking_tstruct(ts::Int64, nx::Int64, fs::Float64)
  Δ = round(Int64, sμ/fs)
  t = rand([
    [1 ts+Δ; nx-1 7000; nx 57425593],
    [1 ts+Δ; nx-2 7000; nx-1 12345; nx 57425593],
    [1 ts+Δ; nx-1 5001; nx 57425593],
    [1 ts+Δ; 10 13412300; 11 123123123; nx-3 3030303; nx-2 -30000000; nx-1 12300045; nx 5700],
    [1 ts+Δ; 10 13412300; 11 123123123; nx-3 303030300000; nx-1 12300045; nx 57425593],
    [1 ts+Δ; 10 13412300; 11 123123123; nx-3 303030300000; nx-1 12300045; nx 57425593],
    [1 ts+Δ; 10 13412300; 11 123123123; nx-3 3030303; nx-2 -30000000; nx-1 12300045; nx 57425593]])
  δ = div(Δ, 2) + 1
  i = findall(t[:,2] .< δ)
  t[i, 2] .= δ

  return t
end
#=
Time matrices beyond program limitations:

(1) Scenario  gap of length |δ| ≤ Δ/2 for any Δ
    Example   δ = div(Δ,2); t = [1 ts+Δ; nx-1 δ; nx 57425593]
    Test      write and reread
    Outcome   last two windows combine
    Cause     gap length in penultimate window too short
=#

function breaking_seis()
  S = SeisData(randSeisData(), randSeisEvent(), randSeisData(2, c=1.0, s=0.0)[2])

  # Test a channel with every possible dict type
  S.misc[1] = breaking_dict

  # Test a channel with no notes
  S.notes[1] = []

  # Need a channel with a very long name to test in show.jl
  S.name[1] = "The quick brown fox jumped over the lazy dog"

  # Need a channel with a non-ASCII filename
  S.name[2] = "Moominpaskanäköinen"
  S.misc[2]["whoo"] = String[]        # ...and an empty String array in :misc
  S.misc[2]["♃♄♅♆♇"] = rand(3,4,5,6)  # ...and a 4d array in :misc

  #= Here we test true, full Unicode support;
    only 0xff can be a separator in S.notes[2] =#
  S.notes[2] = Array{String,1}(undef,6)
  S.notes[2][1] = String(Char.(0x00:0xfe))
  for i = 2:1:6
    uj = randperm(rand(1:n_unicode))
    S.notes[2][i] = join(unicode_chars[uj])
  end

  # Test short data, loc arrays
  S.loc[1] = GenLoc()
  S.loc[2] = GeoLoc()
  S.loc[3] = UTMLoc()
  S.loc[4] = XYLoc()

  # Responses
  S.resp[1] = GenResp()
  S.resp[2] = PZResp()
  S.resp[3] = MultiStageResp(6)
  S.resp[3].stage[1] = CoeffResp()
  S.resp[3].stage[2] = PZResp()
  S.resp[3].gain[1] = 3.5e15
  S.resp[3].fs[1] = 15.0
  S.resp[3].stage[1].a = randn(Float64, 120)
  S.resp[3].stage[1].b = randn(Float64, 120)
  S.resp[3].i[1] = "{counts}"
  S.resp[3].o[1] = "m/s"

  S.x[4] = rand(Float64,4)
  S.t[4] = vcat(S.t[4][1:1,:], [4 0])

  # Some IDs that I can search for
  S.id[1] = "UW.VLL..EHZ"
  S.id[2] = "UW.VLM..EHZ"
  S.id[3] = "UW.TDH..EHZ"

  # Breaking time structures in previous tests
  for i in 1:3
    if S.fs[i] > 0
      S.t[i] = breaking_tstruct(S.t[i][1,2], length(S.x[i]), S.fs[i])
    end
  end
  return S
end