text:
#
#   UNFINISHED !!!
#
#

#
# Intersection of elements from left and right stream
#

#
# this could be done with an inner join,
# computation of intersection, and resorting
#
# Here we are trying to solve the problem using heaps for left and
# right stream. Each stream has a (bpEnd-sorted) heap that contains
# all currently open segments of that stream
#
# Whenever we encounter a new segment (on left or right stream),
# we take the (bpEnd-sorted) heap from the other stream,
# pop off any segments that have already ended,
# and then report the intersection of the new segment with all
# remaining elements on the heap.
# i.e. lelt.Chrom, max(lelt.bpStart, relt.bpStart),
#                  min(lelt.bpEnd, relt.bpEnd)
#
#
# Need to keep track of which stream produced new segment,
# (:left or :right) and which idx of other streams heap is
# to be reported next
#

export gorintersect

using Tables
using DataStructures: heappush!, heappop!

struct Intersection{L,R}
    left::L
    right::R
end

Base.IteratorEltype(::Type{Intersection{L,R}}) where {L,R} = Base.HasEltype()

#Base.eltype(::Type{Intersection{L,R}}) where {L,R} =
#    Base.return_types(mergex, (eltype(L), eltype(R)))[1]

Base.IteratorSize(::Type{Intersection{L,R}}) where {L,R} = Base.SizeUnknown()

Tables.istable(i::Intersection{L,R}) where {L,R} = true
Tables.rowaccess(i::Intersection{L,R}) where {L,R} = true
Tables.schema(i::Intersection{L,R}) where {L,R} = Tables.Schema(eltype(i))
Tables.rows(i::Intersection{L,R}) where {L,R} = i

gorintersect(left, right) = Intersection{typeof(left), typeof(right)}(left, right)

gorintersect(right) = left -> gorintersect(left, right)

function Base.iterate(i::Intersection{L,R}) where {L,R}

    lelt_s = iterate(i.left)
    lelt_s == nothing && return nothing

    relt_s = iterate(i.right)

    lh = Vector{Tuple{String, Int64}, eltype(i.left)}()
    rh = Vector{Tuple{String, Int64}, eltype(i.right)}()

    state = (lelt_s, lh, relt_s, rh, :left, 1)
    iterate(i, state)

end


text:
using Test
using PlayingCards
using TexasHoldem
import TexasHoldem
TH = TexasHoldem

include("tester_bots.jl")

@testset "Game: Play (invalid number of players)" begin
    @test_throws AssertionError("Invalid number of players") Game(ntuple(i->Player(Bot5050(), i), 1))
    @test_throws AssertionError("Invalid number of players") Game(ntuple(i->Player(Bot5050(), i), 11))
end

@testset "Game: Play (BotCheckCall)" begin
    play!(Game(ntuple(i->Player(BotCheckCall(), i), 3)))
end

@testset "Game: Play (BotCheckFold)" begin
    play!(Game(ntuple(i->Player(BotCheckFold(), i), 3)))
end

@testset "Game: Play (BotBetSB) - breaks lowest allowable bet" begin
    game = Game((Player(BotBetSB(),1), Player(BotCheckCall(),2)))
    # Cannot raise small blind on pre-flop! The big blind is the raise min raise, so we must raise at least big-blind
    @test_throws AssertionError("Raise must be between [\$4.0, \$200.0]") play!(game)
end

@testset "Game: Play (BotBetBB)" begin
    play!(Game((Player(BotBetBB(), 1), Player(BotCheckCall(), 2),)))
end

struct NoActionBot <: AbstractAI end

TH.player_option!(game::Game, player::Player{NoActionBot}, ::AbstractGameState, ::PlayerOptions) = nothing

@testset "Game: Play (NoActionBot)" begin
    game = Game((Player(BotCheckCall(), 1), Player(NoActionBot(), 2),))
    @test_throws AssertionError("Must take exactly 1 action.") play!(game)
end

@testset "Non-valid option using BotCheckOnCallRaiseFold" begin
    game = Game((Player(BotCheckCall(), 1), Player(BotCheckOnCallRaiseFold(), 2),))
    @test_throws ErrorException("Cannot check. Available options: CallRaiseFold") play!(game)
end
@testset "Non-valid option using BotCheckOnCallAllInFold" begin
    game = Game((Player(BotCheckOnCallAllInFold(), 1), Player(BotRaiseAlmostAllIn(), 2)))
    @test_throws ErrorException("Cannot check. Available options: CallAllInFold") play!(game)
end
@testset "Non-valid option using BotCheckOnCallFold" begin
    game = Game((Player(BotCheckOnCallFold(), 1), Player(BotRaiseAllIn(), 2)))
    @test_throws ErrorException("Cannot check. Available options: CallFold") play!(game)
end
@testset "Non-valid option using BotCallOnCheckRaiseFold" begin
    game = Game((Player(BotCheckCall(), 1), Player(BotCallOnCheckRaiseFold(), 2),))
    # We catch this incorrect option error before it's completed,
    # so we can't error in `validate_action`
    # @test_throws ErrorException("Cannot call. Available options: CheckRaiseFold") play!(game)
    @test_throws AssertionError("Cannot contribute \$0.0 to the pot!") play!(game)
end
@testset "Non-valid option using BotRaiseOnCallFold" begin
    game = Game((Player(BotRaiseOnCallFold(), 1), Player(BotRaiseAllIn(), 2)))
    @test_throws AssertionError("Cannot raise 0.") play!(game)
end

@testset "Player limps all in" begin
    # Note that the order of the players seems to be important here.
    game = Game((Player(BotCheckCall(), 2; bank_roll = 11), Player(BotLimpAllIn(), 1; bank_roll = 1)))
    play!(game)
end

n_check_actions = [0]
n_call_actions = [0]
@testset "N-actions (3 players)" begin
    play!(Game(ntuple(i->Player(BotNActions(), i), 3)))
    @test n_call_actions[1] == 2 # dealer + small blind
    # 1 check pre-flop (big blind)
    # 3 checks (flop)
    # 3 checks (turn)
    # 3 checks (river)
    @test n_check_actions[1] == 10
end

n_check_actions = [0]
n_call_actions = [0]
@testset "N-actions (4 players)" begin
    play!(Game(ntuple(i->Player(BotNActions(), i), 4)))
    @test n_call_actions[1] == 3 # dealer + small blind + 1 other
    # 1 check pre-flop (big blind)
    # 4 checks (flop)
    # 4 checks (turn)
    # 4 checks (river)
    @test n_check_actions[1] == 13
end

n_check_actions = [0]
n_call_actions = [0]
@testset "N-actions (custom 1)" begin
    # Inspired by:
    # [ Info: Initial bank roll summary: (260.0, 0.0, 37.0, 0.0, 203.0)
    # [ Info: Buttons (dealer, small, big, 1ˢᵗToAct): (1, 3, 5, 1)
    play!(Game((
        Player(BotNActions(), 1),
        Player(BotNActions(), 2; bank_roll = 0),
        Player(BotNActions(), 3),
        Player(BotNActions(), 4; bank_roll = 0),
        Player(BotNActions(), 5),
    ); dealer_id = 1))

    @test n_call_actions[1] == 2 # dealer + small blind
    # 1 check pre-flop (big blind)
    # 3 checks (flop)
    # 3 checks (turn)
    # 3 checks (river)
    @test n_check_actions[1] == 10
end

n_check_actions = [0]
n_call_actions = [0]
@testset "N-actions (custom 2)" begin
    # Inspired by:
    # [ Info: Initial bank roll summary: (195.0, 0.0, 256.0, 4.0, 45.0)
    # [ Info: Buttons (dealer, small, big, 1ˢᵗToAct): (5, 1, 3, 4)
    play!(Game((
        Player(BotNActions(), 1; bank_roll = 195.0),
        Player(BotCheckCall(), 2; bank_roll = 0),
        Player(BotPreFlopRaise(7.0), 3; bank_roll = 256.0),
        Player(BotCheckCall(), 4; bank_roll = 4.0),
        Player(BotCheckFold(), 5; bank_roll = 45.0),
    ); dealer_id = 1))

    @test n_call_actions[1] == 2 # dealer + small blind
    # 1 check pre-flop (big blind)
    # 3 checks (flop)
    # 3 checks (turn)
    # 3 checks (river)
    @test n_check_actions[1] == 3
end

ext: jl


text:
##### OVERRIDES FOR EFFICIENCY / CORRECTNESS

function add_vertices!(g::AbstractSimpleWeightedGraph, n::Integer)
    T = eltype(g)
    U = weighttype(g)
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    emptycols = spzeros(U, nv(g) + n, n)
    g.weights = hcat(g.weights, emptycols[1:nv(g), :])
    g.weights = vcat(g.weights, emptycols')
    return true
end

function degree_matrix(g::AbstractSimpleWeightedGraph, T::DataType=weighttype(g); dir::Symbol=:out)
    if is_directed(g)
        if dir == :out
            d = vec(sum(g.weights, dims=1))
        elseif dir == :in
            d = vec(sum(g.weights, dims=2))
        elseif dir == :both
            d = vec(sum(g.weights, dims=1)) + vec(sum(g.weights, dims=2))
        else
            throw(DomainError(dir, "invalid argument, only accept :in, :out and :both"))
        end
    else
        d = vec(sum(g.weights, dims=1))
    end
    return SparseMatrixCSC(T.(diagm(0=>d)))
end

function adjacency_matrix(g::AbstractSimpleWeightedGraph, T::DataType=weighttype(g); dir::Symbol=:out)
    if dir == :out
        return SparseMatrixCSC(T.(copy(g.weights))')
    else
        return T.(copy(g.weights))
    end
end

function laplacian_matrix(g::AbstractSimpleWeightedGraph, T::DataType=weighttype(g); dir::Symbol=:out)
    degree_matrix(g, T; dir=dir) - adjacency_matrix(g, T; dir=dir)
end

function pagerank(g::SimpleWeightedDiGraph, α=0.85, n::Integer=100, ϵ=1.0e-6)
    A = weights(g)
    S = vec(sum(A, dims=1))
    S = 1 ./ S
    S[findall(S .== Inf)] .= 0.0
    M = A'  # need a separate line due to bug #17456 in julia
    # scaling the adjmat to stochastic adjacency matrix
    M = (Diagonal(S) * M)'
    N = Int(nv(g))
    # solution vector
    x = fill(1.0 / N, N)
    # personalization vector
    p = fill(1.0 / N, N)
    # temporary to hold the results of SpMV
    y = zeros(Float64, N)
    # adjustment for leaf nodes in digraph
    dangling_weights = p
    is_dangling = findall(S .== 0)
    # save some flops by precomputing this
    pscaled = (1 .- α) .* p
    for _ in 1:n
        xlast = x
        # in place SpMV to conserve memory
        mul!(y, M, x)
        # using broadcast to avoid temporaries
        x = α .* (y .+ sum(x[is_dangling]) .* dangling_weights) .+ pscaled
        # l1 change in solution convergence criterion
        err = sum(abs, (x .- xlast))
        if (err < N * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.")
end

savegraph(fn::AbstractString, g::AbstractSimpleWeightedGraph, gname::AbstractString="graph"; compress=true) =
    savegraph(fn, g, gname, SWGFormat(), compress=compress)

savegraph(fn::AbstractString, d::Dict{T, U}; compress=true) where T <: AbstractString where U <: AbstractSimpleWeightedGraph =
    savegraph(fn, d, SWGFormat(), compress=compress)

# It is possible that this is suboptimal, but it is the most trivial extension of the implementation used in Graphs.jl
function cartesian_product(g::G, h::G) where G <: AbstractSimpleWeightedGraph
    z = G(nv(g) * nv(h))
    id(i, j) = (i - 1) * nv(h) + j
    for e in edges(g)
        i1, i2 = Tuple(e)
        for j = 1:nv(h)
            add_edge!(z, id(i1, j), id(i2, j), weight(e))
        end
    end

    for e in edges(h)
        j1, j2 = Tuple(e)
        for i in vertices(g)
            add_edge!(z, id(i, j1), id(i, j2), weight(e))
        end
    end
    return z
end

# Connected Components on a Sparse Matrix

function _cc(g::SimpleWeightedGraph{T,U}) where T where U
    a = weights(g)
    comp = 0
    n = size(a, 1)
    marks = zeros(T, n)
    queue = Vector{T}()
    for i = 1:n
        if marks[i] == 0
            comp += 1
            push!(queue, i)
            while !isempty(queue)
                v = pop!(queue)
                marks[v] = comp
                for index in nzrange(a,v)
                    n = a.rowval[index]
                    if marks[n] == 0
                        push!(queue, n)
                    end
                end
            end
        end
    end
    marks, comp
end

function connected_components(g::SimpleWeightedGraph{T,U}) where T where U
    marks, num_cc = _cc(g)
    cc = [Vector{T}() for i = 1:num_cc]
    for (i,v) in enumerate(marks)
        push!(cc[v], i)
    end
    cc
end

function induced_subgraph(g::T, vlist::AbstractVector{U}) where T <: AbstractSimpleWeightedGraph where U <: Integer
    E = eltype(g)
    allunique(vlist) || throw(ArgumentError("Vertices in subgraph list must be unique"))
    new_weights = g.weights[E.(vlist), E.(vlist)]
    newg = zero(g)
    newg.weights = new_weights
    return newg, Vector{E}(vlist)
end


text:
function baselinefilter(Fₕ, Fₗ=zeros(size(Fₕ)); filtervariance=true, thresh=0.05)
    initialsize = size(Fₕ, 1);

    nanfilter = _nonanrows(Fₕ);
    Fₕ, Fₗ = nanfilter.((Fₕ, Fₗ))

    inffilter = _noinfrows(Fₕ);
    Fₕ, Fₗ = inffilter.((Fₕ, Fₗ))

    integerfilter = _nointegerfeatures(Fₕ, 10, 3) # * Threshold at 10 unique values, to 3 significant figures.
    Fₕ, Fₗ = integerfilter.((Fₕ, Fₗ))

    function lowdimfilter(Fₕ, Fₗ)
        # * Get features that have a low dim variance smaller than thresh% the high dim variance
        σₕ, σₗ = std.((Fₕ, Fₗ), dims=2)
        prop = Array(σₗ)./Array(σₕ)
        idxs = vec(prop .< thresh)
        return F -> ((@warn "Filtered by variances down to $(sum(idxs))/$(length(idxs)) features"),
                    F[idxs, :])[2]
    end
    variancefilter = filtervariance ? lowdimfilter(Fₕ, Fₗ) : identity

    @warn "Filtered down to $(size(Fₕ, 1))/$initialsize features"
    return F -> F |> nanfilter |> inffilter |> integerfilter |> variancefilter # The order of these matters
end
export baselinefilter

function baselinefilter(; kwargs...)
    return (x, y) -> baselinefilter(x, y; kwargs...)
end


text:
module ObjectFile
using Reexport

# Include base utilities
include("utils.jl")
include("string_utils.jl")

# Include our Abstract definitions
include("Abstract/Abstract.jl")

# Include ELF format
include("ELF/ELF.jl")
@reexport using .ELF

# Include MachO format
include("MachO/MachO.jl")
@reexport using .MachO

# Include COFF format
include("COFF/COFF.jl")
@reexport using .COFF

function __init__()
    global ObjTypes

    push!(ObjTypes, ELFHandle)
    push!(ObjTypes, MachOHandle)
    push!(ObjTypes, COFFHandle)
end

end #module ObjectFile


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AvailabilitySetsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Create or update an availability set.
Param: resourceGroupName::String (required)
Param: availabilitySetName::String (required)
Param: parameters::AvailabilitySet (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: AvailabilitySet
"""
function availabilitySetsCreateOrUpdate(_api::AvailabilitySetsApi, resourceGroupName::String, availabilitySetName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", AvailabilitySet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/availabilitySets/{availabilitySetName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "availabilitySetName", availabilitySetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Delete an availability set.
Param: resourceGroupName::String (required)
Param: availabilitySetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: OperationStatusResponse
"""
function availabilitySetsDelete(_api::AvailabilitySetsApi, resourceGroupName::String, availabilitySetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", OperationStatusResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/availabilitySets/{availabilitySetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "availabilitySetName", availabilitySetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Retrieves information about an availability set.
Param: resourceGroupName::String (required)
Param: availabilitySetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: AvailabilitySet
"""
function availabilitySetsGet(_api::AvailabilitySetsApi, resourceGroupName::String, availabilitySetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", AvailabilitySet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/availabilitySets/{availabilitySetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "availabilitySetName", availabilitySetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all availability sets in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: AvailabilitySetListResult
"""
function availabilitySetsList(_api::AvailabilitySetsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", AvailabilitySetListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/availabilitySets", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all available virtual machine sizes that can be used to create a new virtual machine in an existing availability set.
Param: resourceGroupName::String (required)
Param: availabilitySetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineSizeListResult
"""
function availabilitySetsListAvailableSizes(_api::AvailabilitySetsApi, resourceGroupName::String, availabilitySetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineSizeListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/availabilitySets/{availabilitySetName}/vmSizes", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "availabilitySetName", availabilitySetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the list of Microsoft.Compute SKUs available for your Subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ResourceSkusResult
"""
function resourceSkusList(_api::AvailabilitySetsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ResourceSkusResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/skus", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export availabilitySetsCreateOrUpdate, availabilitySetsDelete, availabilitySetsGet, availabilitySetsList, availabilitySetsListAvailableSizes, resourceSkusList


text:
# pyplot()
using Plots
default(grid=false) #src
plot(randn(10))
# closeall();gui()
using Parameters, LTVModelsBase, ValueHistories, DSP, ForwardDiff, AutoGrad #Zygote #Yota#, JLD, ReverseDiff#, JacProp
const Diff = ForwardDiff
# const RDiff = ReverseDiff
@with_kw struct LinearSys
    A
    B
    N     = 1000
    nx    = size(A,1)
    nu    = size(B,2)
    h     = 0.02
    σ0    = 0
    sind  = 1:nx
    uind  = nx+1:(nx+nu)
    s1ind = (nx+nu+1):(nx+nu+nx)
end

function LinearSys(seed; nx = 10, nu = nx, h=0.02, kwargs...)
    srand(seed)
    A = randn(nx,nx)
    A = A-A'        # skew-symmetric = pure imaginary eigenvalues
    A = A - h*I     # Make 'slightly' stable
    A = expm(h*A)   # discrete time
    B = h*randn(nx,nu)
    LinearSys(;A=A, B=B, nx=nx, nu=nu, h=h, kwargs...)
end

function generate_data(sys::LinearSys, seed, validation=false)
    Parameters.@unpack A,B,N, nx, nu, h, σ0 = sys
    srand(seed)
    u      = DSP.filt(ones(5),[5], 10randn(N+2,nu))'
    t      = h:h:N*h+h
    x0     = randn(nx)
    x      = zeros(nx,N+1)
    x[:,1] = x0
    for i = 1:N-1
        x[:,i+1] = A*x[:,i] + B*u[:,i]
    end

    validation || (x .+= σ0 * randn(size(x)))
    u = u[:,1:N]
    @assert all(isfinite, u)
    x,u
end

function true_jacobian(sys::LinearSys, x, u)
    [sys.A sys.B]
end

num_params = 30
wdecay     = 0
stepsize   = 0.02
const sys  = LinearSys(1, nx=5, N=100, h=0.02, σ0 = 0.01)
true_jacobian(x,u) = true_jacobian(sys,x,u)
nu         = sys.nu
nx         = sys.nx

# Generate validation data
function valdata()
    vx,vu,vy = Vector{Float64}[],Vector{Float64}[],Vector{Float64}[]
    for i = 20:60
        x,u = generate_data(sys,i, true)
        for j in 10:5:(sys.N-1)
            push!(vx, x[:,j])
            push!(vy, x[:,j+1])
            push!(vu, u[:,j])
        end
    end
    hcat(vx...),hcat(vu...),hcat(vy...)
end
vx,vu,vy = valdata()
const vt = Trajectory(vx,vu,vy)

# Generate training trajectories
const trajs = [Trajectory(generate_data(sys, i)...) for i = 1:3]
# TODO: take care of separate trajs
function todata(trajs::Vector{Trajectory})
    xu = hcat(getfield.(trajs,:xu)...)
    y = hcat(getfield.(trajs,:y)...)
    xu, y
end

function pred(w,x)
    state = copy(x)
    for i=1:2:length(w)-2
        state = tanh.(w[i]*state .+ w[i+1])
    end
    return w[end-1]*state .+ w[end] #.+ x[1:nx,:]
end

cost(w,x,y)  = sum(abs2.( pred(w,x) - y))/size(y,2)
# cost(w,data) = sum(cost(w, d...) for d in data)/length(data)
cost(w,data) = cost(w, data...)

function loss(w,x,y)
    chunk = Diff.Chunk(x[:,1])
    function lf(w)
        # println("Entering loss function, typeof(w):", typeof(w))
        model(x)    = pred(w,x)
        jcfg        = Diff.JacobianConfig(model, x[:,1], chunk)
        jacobian(x) = Diff.jacobian(model, x, jcfg)
        l = cost(w,x,y)
        J1 = jacobian(x[:,1])
        for t = 2:size(x,2)
            J2 = jacobian(x[:,t])
            l += 2sum(abs2.(J1.-J2))
            J1 = J2
        end
        l
    end
end

##

function runtrain(w, loss; epochs = 500)
    dtrn        = todata(trajs)
    dtst        = todata([vt])
    x,y         = dtrn
    lossfun     = loss(w,x,y)
    # gradient(w) = Diff.gradient(lossfun, w, gcfg)
    trace       = History(Float64)
    tracev      = History(Float64)
    push!(trace, 0, lossfun(w))
    push!(tracev, 0, cost(w,dtst))
    g = similar.(w)
    m = zeros.(size.(g))
    # plot(reuse=false)
    for epoch=1:epochs
        lossfun     = loss(w,x,y)
        γ = 0.85
        # g = ∇(lossfun)(w.-γ.*m) # Nabla
        # val,g = grad(lossfun, w) # Yota
        # @show g = Zygote.gradient(w->cost(w,x,y), w)
        g = AutoGrad.grad(lossfun)(w)
        for i in eachindex(g)
            @. m[i] = 0.002g[i] + γ*m[i]
            @. w[i] -= m[i]
        end
        if epoch % 10 == 0
            push!(trace, epoch, lossfun(w))
            push!(tracev, epoch, cost(w,dtst))
        #     plot(trace, reuse=true)
        #     plot!(tracev)
        #     gui()
        #     println(last(trace), last(tracev))
        end
    end
    trace,tracev
end

Base.promote_rule(::Type{Float64}, ::Type{Yota.TReal}) = Yota.TReal
Base.convert(::Type{Yota.TReal}, f::Float64) = Yota.TReal(f)
const sizes = ((num_params,nx+nu), (num_params,1), (nx,num_params), (nx,1))
srand(1)
w1 = [ 0.1randn(Float64,sizes[1]), zeros(Float64,sizes[2]),0.1randn(Float64,sizes[3]),  zeros(Float64,sizes[4]) ]
# w1 = tovec(w1) # Not for yota
w2 = deepcopy(w1)

# res = runtrain(w1, (w,x,y)-> (w-> cost(w,x,y)), epochs=300)
resj = runtrain(w2, loss, epochs=1000)
# @save "res1.jld" res w1
# @save "resj.jld" resj w2
##

# Eval jac
model(x)    = pred(w2,x,sizes)
jcfg        = Diff.JacobianConfig(model, vt.xu[:,1])
jacobian(x) = Diff.jacobian(model, x, jcfg)
jacs = vcat([eigvals(jacobian(vt.xu[:,i])[1:nx,1:nx]) for i = 1:length(vt)]...)
scatter(real.(jacs), imag.(jacs))
scatter!(real.(eigvals(sys.A)), imag.(eigvals(sys.A)))
phi = linspace(-π/5,π/5,200)
plot!(real.(exp.(phi.*im)), imag.(exp.(phi.*im)), l=(:dash, :black), yaxis=[-0.5,0.5], xaxis=[0,1.2])

@load "resj.jld"
model(x)    = pred(w2,x,sizes)
jfg         = ForwardDiff.JacobianConfig(model, zeros(nx+nu))
jacobian(x) = ForwardDiff.jacobian(model, x)
m = JacProp.ADDiffSystem(w2,sizes,nx,nu, 1, jacobian)
JacProp.eigvalplot(m,vt, true_jacobian)

# Simple test
using BenchmarkTools
using ForwardDiff
const Diff = ForwardDiff
w = randn(20,20)
# x = [randn(20) for i = 1:100]
x = randn(20,100)
f(w::AbstractMatrix,x) = w*x
f(w::AbstractVector,x) = f(reshape(w, 20,20),x)
# const jcfg  = Diff.JacobianConfig(model, x)
function get_loss(w,x)
    modelw(w)    = [f(w,x) for x = x]
    function loss(w)
        modelx(x::AbstractVector)    = f(w,x)
        jacobian(x) = Diff.jacobian(modelx, x)
        # @assert jacobian(x) == w
        # sum(sum(modelw(w))) + sum(sum(jacobian(xi) for xi in x))
        sum(sum(modelw(w))) + sum(sum(jacobian(x[:,i]) for i in 1:size(x,2)))
        # sum(modelw(w)) + sum(jacobian(x))
    end
    loss
end
loss = get_loss(w,x)
grad(w) = Diff.gradient(loss,vec(w))
grad(w)

@btime grad($w)

# f(w,x) = NeuralNetwork(w,x)
# function loss(w,x)
#     ∇xf(x) = jacobian(x->f(w,x), x)
#     w -> norm(∇xf(x))
# end
#
# l = loss(w,x)
# gradient(loss,w) # Good luch with this

function i2m(w,i)
    s = [1; cumsum([prod.(sizes)...]) .+ 1]
    reshape(w[s[i]:(s[i+1]-1)], sizes[i])
end

function pred(w,x)
    state = copy(x)
    for i=1:2:length(sizes)-2
        state = tanh.(i2m(w,i)*state .+ i2m(w,i+1))
    end
    return i2m(w,length(sizes)-1)*state .+ i2m(w,length(sizes))
end

cost(w,x,y)  = sum(abs2.( pred(w,x) .- y))/size(y,2)

function loss(w,x,y)
    chunk = Diff.Chunk(x[:,1])
    function lf(w)
        jacobian(x) = Diff.jacobian(x->pred(w,x), x)
        @show l = cost(w,x,y)
        @show J1 = jacobian(x[:,1])
        for t = 2:size(x,2)
            J2 = jacobian(x[:,t])
            l += 2sum(abs2.(J1.-J2))
            J1 = J2
        end
        l
    end
end


text:
isweekend(dt::Date) = Dates.dayofweek(dt) ∈ (Dates.Friday, Dates.Saturday, Dates.Sunday)

function hasfiveweekend(month::Integer, year::Integer)
    dmin = Date(year, month, 1)
    dmax = dmin + Dates.Day(Dates.daysinmonth(dmin) - 1)
    return count(isweekend, dmin:dmax) ≥ 15
end

months = collect((y, m) for y in 1900:2100, m in 1:12 if hasfiveweekend(m, y))

println("Number of months with 5 full-weekends: $(length(months))")
println("First five such months:")
for (y, m) in months[1:5] println(" - $y-$m") end
println("Last five such months:")
for (y, m) in months[end-4:end] println(" - $y-$m") end

# extra credit
yrs = getindex.(months, 1)
nyrs = 2100 - 1899 - length(unique(yrs))

println("Number of year with not one 5-full-weekend month: $nyrs")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export grad, jacobian, jvp, j′vp, to_vec
replace_arg(x, xs::Tuple, k::Int) = ntuple(p -> p == k ? x : xs[p], length(xs))

"""
    grad(fdm, f, xs...)

Approximate the gradient of `f` at `xs...` using `fdm`. Assumes that `f(xs...)` is scalar.
"""
function grad end

function grad(fdm, f, x::AbstractArray{T}) where T <: Number
    dx = similar(x)
    tmp = similar(x)
    for k in eachindex(x)
        dx[k] = fdm(zero(T)) do ϵ
            tmp .= x
            tmp[k] += ϵ
            return f(tmp)
        end
    end
    return (dx, )
end

grad(fdm, f, x::Real) = (fdm(f, x), )
grad(fdm, f, x::Tuple) = (grad(fdm, (xs...)->f(xs), x...), )

function grad(fdm, f, d::Dict{K, V}) where {K, V}
    dd = Dict{K, V}()
    for (k, v) in d
        function f′(x)
            tmp = copy(d)
            tmp[k] = x
            return f(tmp)
        end
        dd[k] = grad(fdm, f′, v)[1]
    end
    return (dd, )
end

function grad(fdm, f, x)
    v, back = to_vec(x)
    return (back(grad(fdm, x->f(back(v)), v)), )
end

function grad(fdm, f, xs...)
    return ntuple(length(xs)) do k
        grad(fdm, x->f(replace_arg(x, xs, k)...), xs[k])[1]
    end
end

"""
    jacobian(fdm, f, xs::Union{Real, AbstractArray{<:Real}}; len::Int=length(f(x)))

Approximate the Jacobian of `f` at `x` using `fdm`. `f(x)` must be a length `len` vector. If
`len` is not provided, then `f(x)` is computed once to determine the output size.
"""
function jacobian(fdm, f, x::Union{T, AbstractArray{T}}; len::Int=length(f(x))) where {T <: Number}
    J = Matrix{float(T)}(undef, len, length(x))
    for d in 1:len
        gs = grad(fdm, x->f(x)[d], x)[1]
        for k in 1:length(x)
            J[d, k] = gs[k]
        end
    end
    return (J, )
end

function jacobian(fdm, f, xs...; len::Int=length(f(xs...)))
    return ntuple(length(xs)) do k
        jacobian(fdm, x->f(replace_arg(x, xs, k)...), xs[k]; len=len)[1]
    end
end

"""
    _jvp(fdm, f, x::Vector{<:Number}, ẋ::AbstractVector{<:Number})

Convenience function to compute `jacobian(f, x) * ẋ`.
"""
_jvp(fdm, f, x::Vector{<:Number}, ẋ::AV{<:Number}) = fdm(ε -> f(x .+ ε .* ẋ), zero(eltype(x)))

"""
    _j′vp(fdm, f, ȳ::AbstractVector{<:Number}, x::Vector{<:Number})

Convenience function to compute `transpose(jacobian(f, x)) * ȳ`.
"""
_j′vp(fdm, f, ȳ::AV{<:Number}, x::Vector{<:Number}) = transpose(jacobian(fdm, f, x; len=length(ȳ))[1]) * ȳ

"""
    jvp(fdm, f, x, ẋ)

Compute a Jacobian-vector product with any types of arguments for which [`to_vec`](@ref)
is defined.
"""
function jvp(fdm, f, (x, ẋ)::Tuple{Any, Any})
    x_vec, vec_to_x = to_vec(x)
    _, vec_to_y = to_vec(f(x))
    return vec_to_y(_jvp(fdm, x_vec->to_vec(f(vec_to_x(x_vec)))[1], x_vec, to_vec(ẋ)[1]))
end
function jvp(fdm, f, xẋs::Tuple{Any, Any}...)
    x, ẋ = collect(zip(xẋs...))
    return jvp(fdm, xs->f(xs...)[1], (x, ẋ))
end

"""
    j′vp(fdm, f, ȳ, x...)

Compute an adjoint with any types of arguments for which [`to_vec`](@ref) is defined.
"""
function j′vp(fdm, f, ȳ, x)
    x_vec, vec_to_x = to_vec(x)
    ȳ_vec, _ = to_vec(ȳ)
    return (vec_to_x(_j′vp(fdm, x_vec->to_vec(f(vec_to_x(x_vec)))[1], ȳ_vec, x_vec)), )
end
j′vp(fdm, f, ȳ, xs...) = j′vp(fdm, xs->f(xs...), ȳ, xs)[1]

"""
    to_vec(x) -> Tuple{<:AbstractVector, <:Function}

Transform `x` into a `Vector`, and return a closure which inverts the transformation.
"""
to_vec(x::Number) = ([x], first)

# Vectors
to_vec(x::Vector{<:Number}) = (x, identity)
function to_vec(x::Vector)
    x_vecs_and_backs = map(to_vec, x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    return vcat(x_vecs...), function(x_vec)
        sz = cumsum([map(length, x_vecs)...])
        return [backs[n](x_vec[sz[n]-length(x_vecs[n])+1:sz[n]]) for n in eachindex(x)]
    end
end

# Arrays
to_vec(x::Array{<:Number}) = vec(x), x_vec->reshape(x_vec, size(x))
function to_vec(x::Array)
    x_vec, back = to_vec(reshape(x, :))
    return x_vec, x_vec->reshape(back(x_vec), size(x))
end

# AbstractArrays
function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    return x_vec, x_vec->T(reshape(back(x_vec), size(x)))
end
to_vec(x::Symmetric) = vec(Matrix(x)), x_vec->Symmetric(reshape(x_vec, size(x)))
to_vec(X::Diagonal) = vec(Matrix(X)), x_vec->Diagonal(reshape(x_vec, size(X)...))

function to_vec(X::Transpose)
    return vec(Matrix(X)), x_vec->Transpose(permutedims(reshape(x_vec, size(X))))
end
function to_vec(X::Adjoint)
    return vec(Matrix(X)), x_vec->Adjoint(conj!(permutedims(reshape(x_vec, size(X)))))
end

# Non-array data structures

function to_vec(x::Tuple)
    x_vecs, x_backs = zip(map(to_vec, x)...)
    sz = cumsum([map(length, x_vecs)...])
    return vcat(x_vecs...), function(v)
        return ntuple(n->x_backs[n](v[sz[n]-length(x_vecs[n])+1:sz[n]]), length(x))
    end
end

# Convert to a vector-of-vectors to make use of existing functionality.
function to_vec(d::Dict)
    d_vec_vec = [val for val in values(d)]
    d_vec, back = to_vec(d_vec_vec)
    return d_vec, function(v)
        v_vec_vec = back(v)
        return Dict([(key, v_vec_vec[n]) for (n, key) in enumerate(keys(d))])
    end
end


text:
"""
Abstract type for geodetic datums

A datum is a set of reference objects and assigned coordinates, relative to
which other objects may be positioned.  We model these in code with subtypes of
`Datum`.  Each geodetic datum has an associated ellipsoid model of the Earth
which is required when transforming between coordinate systems.  The ellipsoid
can be accessed with the `ellipsoid()` function.
"""
abstract Datum

#------------------------------------------------------------------------------
# Worldwide geodetic datums

const valid_WGS84_frame_weeks = [0, 730,  873, 1150, 1674, 1762]
const WGS84_frames_last_updated = 1907

"""
    WGS84()

Construct an object representing the World Geodetic System of 1984, as a dynamic
datum - see below for gory details about specific WGS84 frames).

If you're getting positions from a consumer GPS device, you're probably going to
have WGS84 by default because it's the datum in which GPS satellites
broadcast their position ("broadcast ephemerides").  Note however that many
devices can also provide position in a national datum, so you should check your
device settings to be sure.

As a special case for low accuracy work (worse than a meter or so), `WGS84` will
assume that coordinates supplied without a capture time are in the *latest*
frame realization known to Geodesy.jl, WGS84 (G$(valid_WGS84_frame_weeks[end])).
Note that this may not be correct if you're processing historical data, or
Geodesy.jl itself is out of date.  For higher accuracy, you should supply a date
of capture - either with each coordinate, or explicitly using the `GpsWeek`
parameter to the WGS84 type:

    WGS84{GpsWeek}()

Construct an object representing the static WGS84 datum computed using data
gathered prior to the given `GpsWeek`.  WGS84 is maintained and updated by the
US National Geospatial-Intelligence-Agency (NGA) at irregular intervals to align
with the ITRF to within 0.1m (see [1]); if you care about accuracy at that
level, you probably want to be solving for position in a different datum, for
example, ITRF.  As of 2016, `GpsWeek` should be one out of
$valid_WGS84_frame_weeks.

Note that the dates of implementation of these frames as broadcast by the
satellites are not the same as the associated GPS week - see Ref. [1], table
2.1.  (TODO: Perhaps Geodesy should have a table to figure out which frame to
use at a given date?  Does anybody care?)

1. "World Geodetic System 1984", NGA standard NGA.STND.0036_1.0.0_WGS84, 2014-07-08,
   http://earth-info.nga.mil/GandG/publications/NGA_STND_0036_1_0_0_WGS84/NGA.STND.0036_1.0.0_WGS84.pdf
"""
immutable WGS84{GpsWeek} <: Datum
    WGS84() = check_wgs84_params(new())
end

WGS84() = WGS84{Void}()

check_wgs84_params(wgs84::WGS84{Void}) = wgs84

@generated function check_wgs84_params{GpsWeek}(wgs84::WGS84{GpsWeek})
    if GpsWeek <= WGS84_frames_last_updated && GpsWeek ∉ valid_WGS84_frame_weeks
        :(throw(ErrorException("No WGS84 realization exists for week $GpsWeek")))
    else
        :(wgs84)
    end
end

Base.show(io::IO, ::WGS84{Void}) = print(io,"WGS84")
Base.show{W}(io::IO, ::WGS84{W}) = print(io,"WGS84 (G$W)")

ellipsoid(::Type{WGS84}) = wgs84_ellipsoid
ellipsoid(::WGS84)       = wgs84_ellipsoid
ellipsoid{GpsWeek}(::Type{WGS84{GpsWeek}}) = wgs84_ellipsoid

"""
    ITRF{Year}([epoch])

Construct an object representing the International Terrestrial Reference Frame
for the given `Year` of realization.  ITRF is the standard high accuracy
terrestrial reference frame for worldwide use.  An optional `epoch` parameter
defines the time of interest - typically a date at which coordinates were
measured, using, eg a GPS device.  Without the `epoch` parameter, the resulting
`ITRF{Year}` object represents the full dynamic datum.

A *realization* is created every few years by computing the position of a large
set of ground control stations from satellite and celestial measurements.  The
`Year` parameter represents the last year from which data was used in the
frame processing regression problem.  A list of recent realizations is available
at http://itrf.ensg.ign.fr/ITRF_solutions; as of 2016-07 valid `Year`s were
2014 2008 2005 2000 1997 1996 1994 1993 1992 1991 1990 1989 1988.

See http://itrf.ensg.ign.fr/general.php for a technical overview.  Useful
technical papers:

* "IERS Conventions (2010)", Petit and Luzum (eds.), IERS Technical note No.36,
   Chapter 4, https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html
* "ITRF2008: an improved solution of the international terrestrial reference
   frame", Altamimi et al., J. Geodesy (2011) 85: 457,
   http://dx.doi.org/10.1007/s00190-011-0444-4
* "IGS08: the IGS realization of ITRF2008", Rebischung et al., GPS Solutions
   (2012) 16: 483, http://dx.doi.org/10.1007/s10291-011-0248-2,
   ftp://igs.org/pub/resource/pubs/IGS08_The_IGS_Realization_of_ITRF2008.pdf
* "The IGS contribution to ITRF2014", Rebischung et al., J. Geodesy (2016) 90:
   611, http://dx.doi.org/10.1007/s00190-016-0897-6
"""
immutable ITRF{Year, EpochT} <: Datum
    epoch::EpochT

    function ITRF(epoch::EpochT)
        check_itrf_year(new(epoch))
    end
end

@generated function check_itrf_year{Y}(itrf::ITRF{Y})
    if Y < 100
        :(throw(ErrorException("Two digit year $Y for ITRF found - this library expects a full four-digit year.")))
    elseif Y <= 2016 && Y ∉ [2014, 2008, 2005, 2000, 1997, 1996,
                             1994, 1993, 1992, 1991, 1990, 1989, 1988]
        :(throw(ErrorException("No ITRF realization exists for year $Y")))
    else
        :(itrf)
    end
end

(::Type{ITRF{Year}}){Year}() = ITRF{Year,Void}(nothing)
(::Type{ITRF{Year}}){Year}(epoch) = ITRF{Year,typeof(epoch)}(epoch)

Base.show{Y}(io::IO, ::ITRF{Y,Void}) = print(io,"ITRF{$Y}")
Base.show{Y}(io::IO, itrf::ITRF{Y}) = print(io,"ITRF{$Y}($(itrf.epoch))")

ellipsoid(::ITRF)             = grs80
ellipsoid{D<:ITRF}(::Type{D}) = grs80

#------------------------------------------------------------------------------
# National geodetic datums
"""
`OSGB36` - Datum for Ordinance Survey of Great Britain, 1936
"""
immutable OSGB36 <: Datum; end
Base.show(io::IO, ::OSGB36) = print(io,"osbg84")
ellipsoid(::Union{OSGB36,Type{OSGB36}}) = airy1830

"""
`NAD27` - North American Datum of 1927
"""
immutable NAD27 <: Datum; end
Base.show(io::IO, ::NAD27) = print(io,"nad27")
ellipsoid(::Union{NAD27,Type{NAD27}}) = clarke1866

"""
`NAD83` - North American Datum of 1983

For technical details, see "NAD83 (NSRS2007) National Readjustment Final Report"
http://www.ngs.noaa.gov/PUBS_LIB/NSRS2007/NOAATRNOSNGS60.pdf
"""
immutable NAD83 <: Datum; end
Base.show(io::IO, ::NAD83) = print(io,"nad83")
ellipsoid(::Union{NAD83,Type{NAD83}}) = grs80

"""
`GDA94` - Geocentric Datum of Australia, 1994
"""
immutable GDA94 <: Datum; end
ellipsoid(::Union{GDA94,Type{GDA94}})   = grs80

#-------------------------------------------------------------------------------
# Datum instances
const wgs84 = WGS84()
const osgb36 = OSGB36()
const nad27 = NAD27()
const nad83 = NAD83()
const gda94 = GDA94()


text:
### Exports ###

export SymmetricRuleLearn, PowerLawLearn, ExpWeightDepLearn, ConvSynapses

### Types ###

# TODO: Use Value type to select additional modes?
mutable struct ConvSynapses{L} <: AbstractSynapses
  inputSize::Tuple{Int,Int,Int}
  filterSize::Int
  numFilters::Int
  stride::Int

  learn::L

  W::Array{Float32,4}
end
function ConvSynapses(inputSize::Tuple{Int,Int,Int}, filterSize::Int, numFilters::Int, stride::Int=1; learn=SymmetricRuleLearn())
  return ConvSynapses(
    inputSize,
    filterSize,
    numFilters,
    stride,
    learn,
    rand(Float32, inputSize[1], filterSize, filterSize, numFilters) .* 0.3f0
  )
end
ConvSynapses(size::Tuple{Tuple{Int,Int,Int}, Int, Int, Int}; kwargs...) =
  ConvSynapses(size...; kwargs...)

### Methods ###

function clear!(synapses::ConvSynapses)
  rand!(synapses.W, 0f0:0.01f0:0.3f0)
end
# TODO: Make one-liner
function Base.show(io::IO, synapses::ConvSynapses)
  outputSize = output_size(synapses)
  println(io, "ConvSynapses ($(synapses.inputSize) => $(outputSize))")
  println(io, "  Filter Size: $(synapses.filterSize) x $(synapses.filterSize)")
  println(io, "  Stride: $(synapses.stride)")
  println(io, "  Learning Algorithm: $(synapses.learn)")
end
Base.size(cs::ConvSynapses) = (cs.inputSize, cs.filterSize, cs.numFilters, cs.stride)
function output_size(cs::ConvSynapses)
  # FIXME: Stride
  return (cs.inputSize[2]-cs.filterSize+1, cs.inputSize[3]-cs.filterSize+1, cs.numFilters)
end

function update!(f, cs, stride, W, I, G, F, O)
  @inbounds for fj = axes(W, 3)
    @inbounds for fi = axes(W, 2)
      @inbounds for fc = axes(W, 1)
        @inbounds for oj = axes(O, 2)
          @inbounds for oi = axes(O, 1)
            I_ = Float32(I[fc, oi+fi-1, oj+fj-1]) #Float32(I[fc, (stride*(oi-1))+fi, (stride*(oj-1))+fj])

            # Convolve weights with inputs
            O[oi,oj,f] = muladd(W[fc,fi,fj,f], I_, O[oi,oj,f])
          end
        end
      end
    end
  end
end
function learn!(f, cs, stride, W, I, G, F, O)
  lrn = cs.learn
  α_pre, α_post, xtar, Wmax, μ = lrn.α_pre, lrn.α_post, lrn.xtar, lrn.Wmax, lrn.μ
  @inbounds for oj = axes(O, 2)
    @inbounds for oi = axes(O, 1)
      @inbounds for fj = axes(W, 3)
        @inbounds for fi = axes(W, 2)
          @inbounds @simd for fc = axes(W, 1)
            I_ = Float32(I[fc, oi+fi-1, oj+fj-1]) # TODO: Enable stride: Float32(I[fc, (stride*(oi-1))+fi, (stride*(oj-1))+fj])

            # Learn
            W[fc,fi,fj,f] += learn!(cs.learn, I_, G[oi,oj,f], F[oi,oj,f], W[fc,fi,fj,f])
          end
        end
      end
    end
  end
end
# TODO: Specify input and output types
function _eforward!(scont::CPUContainer{S}, input, output) where S<:ConvSynapses
  cs = root(scont)
  stride, W = cs.stride, cs.W

  I = transient(input).data     #@param input[O]
  G = transient(output).state.T #@param output[T]
  F = transient(output).state.F #@param output[F]
  O = transient(output).state.I #@param output[I]

  # Clear outputs
  fill!(O, 0f0)

  # Convolve
  @inbounds for f = axes(W, 4)
    update!(f, cs, stride, W, I, G, F, O)
    learn!(f, cs, stride, W, I, G, F, O)
  end

  # Clamp weights
  # TODO: Dispatch on type (or just pull into learn!)
  clamp!(cs.W, 0f0, cs.learn.Wmax)
end


text:
########################################
###all hamiltonian-building functions###
########################################
"""
    buildLocalHamiltonian(hx, hz, maxLen)

Return an array of local hamiltonian :X_jX_{j+1} + 0.5*hx*(X[j]+X[j+1])+ 0.5*hz*(Z[j]+Z[j+1])
"""
function buildLocalHamiltonian(hx, hz, maxLen)
    dicHam=Dict()
    for site=1:maxLen-1
        dicHam[site]=spzeros(2^maxLen, 2^maxLen)
        for x=1:2^maxLen
            dig=2*digits(x-1, 2, maxLen)-1
            ###Kinetic terms and sz
            dicHam[site][x, x-dig[site]*2^(site-1)-dig[site+1]*2^(site)]+=0.25
            dicHam[site][x, x-dig[site]*2^(site-1)]+=hx/2/2
            dicHam[site][x, x]+=hz/2*dig[site]/2
        end
    end
    return dicHam
end

"""
    buildInsideHamiltonian(hx, hz, maxLen)

Return the Hamiltonian acting only on site  1:maxLen by dumping extremity links
"""
function buildInsideHamiltonian(hx, hz, maxLen)
    dicHam=spzeros(2^maxLen, 2^maxLen)
    for x=1:2^maxLen
        dig=2*digits(x-1, 2, maxLen)-1
        for site=1:maxLen-1
            dicHam[x, x-dig[site]*2^(site-1)-dig[site+1]*2^(site)]+=0.25
        end
        for site=1:maxLen
            dicHam[x, x]+=hz/2*dig[site]
            dicHam[x, x-dig[site]*2^(site-1)]+=hx/2
        end
    end
    return dicHam
end

"""
    buildLinkHamiltonian(hx, hz, maxLen)

Return the array (left-right) of links: h_{1:2} and h_{maxLen-1:maxLen} for states defined on 1:maxLen
"""
function buildLinkHamiltonian(hx, hz, maxLen)
    dicHam=Array{typeof(spzeros(2,2))}(2)
    dicHam[1]=spzeros(2^maxLen, 2^maxLen);  dicHam[2]=spzeros(2^maxLen, 2^maxLen)
    for x=1:2^maxLen
        dig=2*digits(x-1, 2, maxLen)-1
        dicHam[1][x, x-dig[1]-dig[2]*2]+=0.25
        dicHam[2][x, x-dig[end-1]*2^(maxLen-2)-dig[end]*2^(maxLen-1)]+=0.25
    end
    return dicHam
end

"""
    buildSingleSiteHamiltonian(hx, hz, maxLen)

Build the single site part of the Hamiltonian for sites 1:maxLen
Here, special case as only Z field is taken into account
"""
function buildSingleSiteHamiltonian(hx, hz, maxLen)
    dicHam=spzeros(2^maxLen, 2^maxLen)
    for x=1:2^maxLen
        dig=2*digits(x-1, 2, maxLen)-1
        for site=1:maxLen
            dicHam[x, x]+=hz/2*dig[site]
        end
    end
    return dicHam
end

"""
    buildMultiSiteHamiltonian(hx, hz, maxLen)

Return the multisite part of the Hamiltonian for length 1:maxLen
Special case here as the X field is taken into account here directly
"""
function buildMultiSiteHamiltonian(hx, hz, maxLen)
    dicHam=spzeros(2^maxLen, 2^maxLen)
    for x=1:2^maxLen
        dig=2*digits(x-1, 2, maxLen)-1
        for site=1:maxLen-1
            dicHam[x, x-dig[site]*2^(site-1)-dig[site+1]*2^(site)]+=0.25
        end
        for site=1:maxLen
            dicHam[x, x-dig[site]*2^(site-1)]+=hx/2
        end
    end
    return dicHam
end

"""
    generateInitialState_z(maxLen)

Project the spin 0 on +_z and return all DM of length maxLen containing it
"""
function generateInitialState_z(maxLen; ini=1.0)
    listRho=Array{Array{Complex128}}(maxLen)
    localRho=[1.0-ini+0.0im 0.0; 0.0 ini]
    for l in 1:maxLen
        listRho[l]=kron(kron(eye(Complex128, 2^(l-1)),  localRho), eye(Complex128, 2^(maxLen-l)))/2^(maxLen-1)
    end
    return listRho
end

"""
    generateInitialState_x(maxLen)

Project the spin 0 on +_x and return all DM of length maxLen containing it
"""
function generateInitialState_x(maxLen)
    listRho=Array{Array{Complex128}}(maxLen)
    localRho=[0.5+0.0im 0.5; 0.5 0.5]
    for l in 1:maxLen
        listRho[l]=kron(kron(eye(Complex128, 2^(l-1)),  localRho), eye(Complex128, 2^(maxLen-l)))/2^(maxLen-1)
    end
    return listRho
end

"""
    generateInitialState_z_center(maxLen)

Project the spin 0 on +_z and return the DM of length maxLen where it is centered at maxLen÷2 +1
"""
function generateInitialState_z_center(maxLen)
    localRho=[0.0+0.0im 0.0; 0.0 1.0]
    len=(maxLen÷2)
    listRho=kron(kron(eye(Complex128, 2^(len)),  localRho), eye(Complex128, 2^(maxLen-len-1)))/2^(maxLen-1)
    return listRho
end

max_stars_repo_name: LoicHerviou/DensitySimulation

max_issues_repo_name: LoicHerviou/DensitySimulation

max_forks_repo_name: LoicHerviou/DensitySimulation


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie
import VideoIO
import Images

using Porta

startframe = 1
FPS = 60
resolution = (3840, 2160)
segments = 36
speed = 1
factor = 0.0025
scale = 1.0
radius = 0.005
basemapradius = 1.015
basemapsegments = 2segments
basepointradius = 0.015
basepointsegments = 10
basepointtransparency = false
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "hopfdance"

# Loping Sting by Kevin MacLeod is licensed under a Creative Commons Attribution
# 4.0 license. https://creativecommons.org/licenses/by/4.0/

# Source: http://incompetech.com/music/royalty-free/index.html?isrc=USUAN1200014

# Artist: http://incompetech.com/

audioname = "eric"
extension = ".wav"
audiopath = joinpath("data", audioname * extension)
signal = Signal(audiopath)
chunkspersecond = FPS
#totalsamples = Integer(getframerate(signal) ÷ chunkspersecond) - 1
totalsamples = length(getdata(getfftchunk(signal,1, chunkspersecond, 3)))
number = totalsamples
#indices = convert(Array{Int64},
#                  floor.(range(2, stop = totalsamples-1, length = number)))
frames = countchunks(signal, chunkspersecond)

videopath = "data/eric.mp4"
videoframes = VideoIO.openvideo(videopath)

# The scene object that contains other visual objects
#Makie.reasonable_resolution() = (800, 800)
scene = Makie.Scene(backgroundcolor = :black,
                               show_axis = false,
                               resolution = resolution)

function composeframes(image1, image2)
    height1, width1 = size(image1)
    height2, width2 = size(image2)
    scalefactor = height2 / height1
    newsize = (Int(floor(height1 * scalefactor)), Int(floor(width1 * scalefactor)))
    image3 = Images.imresize(image1, newsize)
    height3, width3 = size(image3)
    pad = Int((width2 - width3) ÷ 2)
    image4 = deepcopy(image2)
    image4[:, pad:pad+width3-1] = image2[begin:end, pad:pad+width3-1] .* 0.5 + image3 .* 0.5
    image4
end

fmap(b::S²) = b

"""
λ⁻¹map(p)
Sends a point on the plane back to a point on a unit sphere with the given
point. This is the inverse stereographic projection of a 3-sphere.
"""
function λ⁻¹map(p::S²)
    r3 = ℝ³(Cartesian(p))
    p₁, p₂, p₃ = vec(r3)
    mgnitude² = norm(r3)^2
    x₁ = 2p₁ / (1 + mgnitude²)
    x₂ = 2p₂ / (1 + mgnitude²)
    x₃ = 2p₃ / (1 + mgnitude²)
    x₄ = (-1 + mgnitude²) / (1 + mgnitude²)
    Quaternion(x₁, x₂, x₃, x₄)
end

"""
    getpointonpath(t)

Get a point on a path on S² with the given parameter `t`. t ∈ [0, 1].
"""
getpointonpath(t::Real) = begin
    t = t * 12π
    d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
    x, y = sin(t) * d, cos(t) * d
    ComplexLine(x + im * y)
end

"""
    getbutterflycurve(points)

Get butterfly curve by Temple H. Fay (1989) with the given number of `points`.
"""
function getbutterflycurve(points::Int)
    array = Array{ComplexLine,1}(undef, points)
    # 0 ≤ t ≤ 12π
    for i in 1:points
        t = (i - 1) / points * 12π
        d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
        x, y = sin(t) * d, cos(t) * d
        array[i] = ComplexLine(x + im * y)
    end
    array
end

quaternionicrgb(p::Quaternion, q::Quaternion, s::Real, v::Real) = begin
    hue = acos(dot(p, q)) / π
    hsvtorgb([hue * 360; s; v])
end

"""
    getpoints(z [, segments])

Get a circle of point in the base space with the given center `z` in the Riemann
sphere. The optional argument `segments` determines how many points should the
circle have.
"""
function getpoints(z::S²; segments::Int = 30)
    g = Geographic(z)
    r, ϕ, θ = vec(g)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    array = Array{Geographic,1}(undef, segments)
    for i in 1:segments
        ϕ′, θ′ = g.ϕ + factor * cos(lspace[i]), g.θ + factor * sin(lspace[i])
        array[i] = Geographic(r, ϕ′, θ′)
    end
    array
end

basemapconfig = Biquaternion(ℝ³(0, 1.1, 0))
# basemapconfig = Biquaternion(Quaternion(1, 0, 0, 0))
#basemapcolor = Makie.RGBAf0(0.75, 0.75, 0.75, 0.25)
#basemapcolor = FileIO.load("gallery/5672_mars_4k_color.jpg")
mars = Images.load("data/mars.jpg")
basemapcolor = mars
basemaptransparency = false
basemap = RGBSphere(basemapconfig,
                    scene,
                    basemapcolor,
                    radius = basemapradius,
                    segments = basemapsegments,
                    transparency = basemaptransparency)

curvepoints = getbutterflycurve(number)
solidcolors = Array{Makie.RGBAf,1}(undef, number)
ghostcolors = similar(solidcolors)
basepoints = Array{Sphere,1}(undef, number)
for i in 1:number
    config = Biquaternion(ℝ³(Cartesian(curvepoints[i])) * basemapradius)
    hue = (i - 1) / number * 360
    rgb = hsvtorgb([hue, 1, 1])
    solidcolors[i] = Makie.RGBAf(rgb..., 0.9)
    ghostcolors[i] = Makie.RGBAf(rgb..., 0.1)
    sphere = Sphere(basemapconfig * config,
                    scene,
                    radius = basepointradius,
                    segments = basepointsegments,
                    color = solidcolors[i],
                    transparency = basepointtransparency)
    basepoints[i] = sphere
end

#bundleconfig = Biquaternion(ℝ³(0, 1.25, 0))
bundleconfig = Biquaternion(ℝ³(0, -1.1, 0))
s3rotation = Quaternion(1, 0, 0, 0)
solidfibers = Array{Fiber,1}(undef, number)
ghostfibers = Array{Fiber,1}(undef, number)
α = 40 / 180 * pi
solidtop = U1(pi - 2α)
solidbottom = U1(-pi)
ghosttop = U1(pi)
ghostbottom = solidtop
for i in 1:number
    point = curvepoints[i]
    transparency = false
    solidfiber = Fiber(scene,
                       point,
                       λ⁻¹map,
                       fmap,
                       radius = radius,
                       top = solidtop,
                       bottom = solidbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = solidcolors[i],
                       transparency = transparency,
                       scale = scale)
    solidfibers[i] = solidfiber
    transparency = true
    ghostfiber = Fiber(scene,
                       point,
                       λ⁻¹map,
                       fmap,
                       radius = radius,
                       top = ghosttop,
                       bottom = ghostbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = ghostcolors[i],
                       transparency = transparency,
                       scale = scale)
    ghostfibers[i] = ghostfiber
end

basespacepoints = Array{ComplexLine,1}(undef, number)
phases = Array{U1,2}(undef, number, 4)
s3rotations = Array{Quaternion,1}(undef, number)
basespacecircles = Array{Any,1}(undef, number)

len = 0.1
width = 0.01
transparency = false
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrows = []
for i in 1:number
    hue = (i - 1) / number * 360
    rgb = hsvtorgb([hue, 1, 1])
    color = Makie.RGBAf(rgb..., 0.9)
    push!(arrows, Arrow(tail, head, scene, width = width, color = color))
end

a, b = ComplexLine(1 + im), ComplexLine(0.5 - 1.5 * im)
c, d = ComplexLine(-0.5 + 1 * im), ComplexLine(-1 + 0.3 * im)
f(z, a, b, c, d) = (a * z + b) / (c * z + d)

videoframe = read(videoframes)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    frequencies = getdata(getfftchunk(signal, i, chunkspersecond, 3))#[indices]
    if i % 2 == 0
        try
            read!(videoframes, videoframe)
            println("Read the new video frame.")
        catch
            println("Reached the end of video stream.")
        end
    end
    basemapcolor = composeframes(videoframe, mars)
    println(basemapcolor == mars)
    update(basemap, basemapcolor)
    step = (i - 1) / frames
    step = 1 - step
    u = getpointonpath(step)
    if isnan(vec(u)[1]) || isnan(vec(u)[2])
        u = Geographic(1, 0, π/2)
    end
    τ = step * speed * 2pi
    u = ℝ³(Cartesian(u))
    q = Quaternion(τ, u)
    basemapconfig = Biquaternion(q, ℝ³(0, 1.1, 0))
    update(basemap, basemapconfig)
    a′, b′, c′, d′ = [ComplexLine(Cartesian(rotate(ℝ³(Cartesian(point)), q))).z for point in (a, b, c, d)]
    s3rotation = Quaternion(1, 0, 0, 0)
    for (index, item) in enumerate(frequencies)
        z = item
        z = ComplexLine(f(2π * z, a′, b′, c′, d′))
        if isnan(vec(z)[1]) || isnan(vec(z)[2])
            z = ComplexLine(Geographic(1, 0, -π / 2))
        end
        r, ϕ, θ = vec(Geographic(z))
        if isnan(r)
            z = ComplexLine(Geographic(1, ϕ, θ * 0.99))
        end

        basespacepoints[index] = z
        α = angle(z) + π
        fraction = α / 2π
        complement = 1 - fraction
        g = Quaternion(τ, ℝ³(0, 0, 1))
        q = λ⁻¹map(z)
        rgb = quaternionicrgb(g, q, 0.2 * fraction + 0.2 * rand() + 0.6, 0.2 * complement + 0.2 * rand() + 0.6)
        solidcolors[index] = Makie.RGBAf(rgb..., fraction)
        ghostcolors[index] = Makie.RGBAf(rgb..., complement)
        solidtop = U1(0)
        solidbottom = U1(α)
        ghosttop = solidbottom
        ghostbottom = U1(2π)
        phases[index, 1] = solidtop
        phases[index, 2] = solidbottom
        phases[index, 3] = ghosttop
        phases[index, 4] = ghostbottom
        #s3rotations[index] = g
    end
    for (index, item) in enumerate(basepoints)
        z = basespacepoints[index]
        p = ℝ³(Cartesian(z)) * basemapradius
        config = Biquaternion(p)
        update(item, basemapconfig * config)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(solidfibers)
        solidtop = phases[index, 1]
        solidbottom = phases[index, 2]
        update(item, basespacepoints[index], s3rotation, radius, solidtop,
               solidbottom)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(ghostfibers)
        ghosttop = phases[index, 3]
        ghostbottom = phases[index, 4]
        update(item, basespacepoints[index], s3rotation, radius, ghosttop,
               ghostbottom)
        update(item, ghostcolors[index])
    end

    # update tangent vectors
    for (index, item) in enumerate(arrows)
        w = basespacepoints[index]
        tail = ℝ³(Cartesian(w)) * basemapradius
        initial = ℝ³(0, 0, 1) * basemapradius
        q = getrotation(normalize(tail), normalize(initial))
        α = angle(phases[index, 2])
        factor = max(0.01, α / 2π) / 5
        head = initial + ℝ³(1, 0, 0) * factor
        head = rotate(head, q)
        q = Quaternion(α / 2, normalize(tail))
        head = rotate(head, q)
        tail, head = applyconfig([tail; head], basemapconfig)
        update(item, tail, head - tail)
        update(item, solidcolors[index])
    end

#=     z = Complex(0 + 0im)
    z = ComplexLine(f(π * z, a′, b′, c′, d′))
    z = ℝ³(Cartesian(z))

    z₀ = Complex(0 + im)
    z₀ = ComplexLine(f(π * z₀, a′, b′, c′, d′))
    z₀ = ℝ³(Cartesian(z₀)) =#

    #n = normalize(z₀ - z)
    #v = ℝ³(-1, 0, 1) * 2.5
    #v = ℝ³(-1, 0, 1) * 2.1
    #distance = norm(v)
    #v = normalize(z) * distance
    # update eye position
    # scene.camera.eyeposition.val
    #upvector = GeometryBasics.Vec3f0(vec(n)...)
    #eyeposition = GeometryBasics.Vec3f0(vec(v)...)
    #lookat = GeometryBasics.Vec3f0(0, 0, 0)
    #Makie.update_cam!(scene, eyeposition, lookat, upvector)
    #scene.center = false # prevent scene from recentering on display
end

n = ℝ³(0, 0, 1)
v = ℝ³(1, 0, 1) * 2.5
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
camera = Makie.Camera3D(scene)
Makie.update_cam!(scene, camera, eyeposition, lookat, upvector)
#scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    directory = joinpath("gallery", modelname, audioname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(audioname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        FileIO.save(filepath,
                    scene;
                    resolution = resolution,
                    pt_per_unit = 400.0,
                    px_per_unit = 400.0)
        elapsed = time() - start
        println("Saving file $filepath took $elapsed (s).")

        step = (i - 1) / frames
        println("Completed step $(100step).\n")
        sleep(1)

        stitch() = begin
            part = i ÷ (frames / 10)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            #WxH = "$(resolution[1])x$(resolution[2])"
            WxH = "1920x1080"
            commonpart = "$(audioname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            outputwithaudioname = "$(commonpart)_withaudio_$(part).mp4"
            outputwithaudiopath = joinpath(exportdir, outputwithaudioname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
            sleep(1)
            command2 = `ffmpeg -y -i $outputpath -i $audiopath -c:v copy -map 0:v:0 -map 1:a:0 -c:a aac -b:a 192k -shortest $(outputwithaudiopath)`
            run(command2)
        end

        if isapprox(i % (frames / 10), 0)
            stitch()
            sleep(1)
        end
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
"""Testing the SmuMolecule.
The 1 based array indexing in Julia makes this somewhat discordant with the phython
implementation.
"""

module TestSmuMolecule

using ProtoBuf
using Test

using SmuMolecules
using SmuUtilities

using dataset_pb2

export all_tests

"""Make a BondTopology consisting of `n` disconnected Carbon atoms.
"""
function carbon_atoms(n::T)::BondTopology where {T<:Integer}
  result = BondTopology()
  for i in 1:n
    add_atom_by_atomic_number!(result, 6)
  end
  return result
end

"""The simplest multi-atom molecule, CC"""
function test_ethane()
  matching_parameters = MatchingParameters()
  matching_parameters.must_match_all_bonds = false

  ethane = carbon_atoms(2)
  add_bond!(0, 1, 1, ethane)
  scores = [0.1, 1.1, 2.1, 3.1]
  bonds_to_scores = Dict{Tuple{Int32, Int32}, Vector{Float32}}((1, 2) => scores)
  mol = SmuMolecule(ethane, bonds_to_scores, matching_parameters)
  state = generate_search_state(mol)
  @test length(state) == 1
  @test state == [[1, 2, 3, 4]]

  for s in Iterators.product(state...)
    res = place_bonds!(s, mol)
    @test res !== nothing
    # self.assertAlmostEqual(res.score, scores[i])
  end
  true
end

function test_ethane_all_btypes()
  matching_parameters = MatchingParameters()
  matching_parameters.must_match_all_bonds = false

  test_cases = [
    [0, dataset_pb2.BondTopology_BondType.BOND_UNDEFINED],
    [1, dataset_pb2.BondTopology_BondType.BOND_SINGLE],
    [2, dataset_pb2.BondTopology_BondType.BOND_DOUBLE],
    [3, dataset_pb2.BondTopology_BondType.BOND_TRIPLE],
  ]

  for t in test_cases
    numeric_btype = t[1]
    smu_btype = t[2]
    cc = carbon_atoms(2)
    bonds_to_scores = Dict{Tuple{Int32,Int32}, Vector{Float32}}((1, 2) => zeros(Float32, 4))
    bonds_to_scores[(1, 2)][numeric_btype + 1] = 1.0
    mol = SmuMolecule(cc, bonds_to_scores, matching_parameters)
    state = generate_search_state(mol)
    for s in Iterators.product(state...)
      res = place_bonds!(s, mol)
      @test res !== nothing
      if numeric_btype == 0
        @test !hasproperty(res, :bonds)
      else
        @test length(res.bonds) == 1
        @test res.bonds[1].bond_type == smu_btype
      end
    end
  end
  true
end

"""Checking to see if a BondTopology has the right number of bonds is
complicated by the fact that if there are no bonds present, the :bonds
attribute will be missing.
Args:
  bond_topology:
  expected_bonds: the number of bonds expected to be in `bond_topology`.
Returns:
  true if the number of bonds present matches `expected_bonds`.
"""
function bond_count_matches(bond_topology::BondTopology, expected_bonds)::Bool
  if !hasproperty(bond_topology, :bonds)
    return expected_bonds == 0
  end

  return length(bond_topology.bonds) == expected_bonds
end

struct A1A2BtypeScore
  btype1::Int32
  btype2::Int32
  expected_bonds::Int32
  expected_score::Union{Float32, Nothing}
end

function test_propane_all()
  matching_parameters = MatchingParameters()
  matching_parameters.must_match_all_bonds = false

  test_cases = [
    A1A2BtypeScore(0, 0, 0, 2.0),
    A1A2BtypeScore(0, 1, 1, 2.0),
    A1A2BtypeScore(0, 2, 1, 2.0),
    A1A2BtypeScore(0, 3, 1, 2.0),
    A1A2BtypeScore(1, 1, 2, 2.0),
    A1A2BtypeScore(1, 2, 2, 2.0),
    A1A2BtypeScore(1, 3, 2, 2.0),
    A1A2BtypeScore(2, 2, 2, 2.0),
    A1A2BtypeScore(2, 3, 0, nothing),
    A1A2BtypeScore(3, 3, 0, nothing)
  ]
  for t in test_cases
    btype1 = t.btype1
    btype2 = t.btype2
    expected_bonds = t.expected_bonds
    expected_score = t.expected_score
    cc = carbon_atoms(3)
    bonds_to_scores = Dict{Tuple{Int32,Int32}, Vector{Float32}}((1, 2) => zeros(Float32, 4),
                       (2, 3) => zeros(Float32, 4))
    bonds_to_scores[(1, 2)][btype1 + 1] = 1.0
    bonds_to_scores[(2, 3)][btype2 + 1] = 1.0
    mol = SmuMolecule(cc, bonds_to_scores, matching_parameters)
    state = generate_search_state(mol)
    for s in Iterators.product(state...)
      res = place_bonds!(s, mol)
      if expected_score !== nothing
        @test res !== nothing
        @test bond_count_matches(res, expected_bonds)
        @test isapprox(res.score, expected_score)
        if btype1 == 0
          if btype2 > 0
            @test res.bonds[1].bond_type == btype2
          end
        else
          @test Set((btype1, btype2)) == Set((res.bonds[1].bond_type, res.bonds[2].bond_type))
        end
      else
        @test res === nothing
      end
    end
  end
  true
end

function test_operators()
  matching_parameters = MatchingParameters()
  matching_parameters.must_match_all_bonds = false

  cc = carbon_atoms(3)
  add_bond!(0, 1, 1, cc)
  add_bond!(1, 2, 1, cc)
  bonds_to_scores = Dict{Tuple{Int32,Int32}, Vector{Float32}}((1, 2) => zeros(Float32, 4),
                                            (2, 3) => zeros(Float32, 4))
  scores = [1.0, 3.0]
  bonds_to_scores[(1, 2)][1] = scores[1]
  bonds_to_scores[(2, 3)][1] = scores[2]
  mol = SmuMolecule(cc, bonds_to_scores, matching_parameters)
  set_initial_score_and_incrementer!(1.0, Base.:*, mol)
  state = generate_search_state(mol)
  for s in Iterators.product(state...)
    res = place_bonds!(s, mol)
    @test isapprox(res.score, prod(scores))
  end
  true
end

function all_tests()
  @test test_ethane()
  @test test_ethane_all_btypes()
  @test test_propane_all()
  @test test_operators()
end

end  # module TestSmuMolecule

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu