text:
"""
`erdos_renyi_undirected`
========================

Generate an undirected Erdős-Rényi graph.

**The current implementation uses sprand, this may change in the future.**
**Do not depend on this routine for reliable output between versions.**

Input
-----
- `n`: the number of nodes
- `p`: the probability of an edge, or the average degree.
   if \$`p` >= 1\$, then \$`p`\$ is interpreted as an average degree
   instead of a probability. (There is no point in generating
   an Erdős-Rényi graph with probability \$`p`=1\$)

Output
------
- A matrix network type for the Erdős-Rényi graph.

Example
-------
~~~~
# show the connected phase transition
n = 100
avgdegs = linspace(1.,2*log(n),100)
compsizes = map( (dbar) ->
        maximum(scomponents(erdos_renyi_undirected(n,dbar)).sizes),
    avgdegs )
using Plots
unicodeplots()
plot(avgdegs,compsizes,xaxis=("average degree"),yaxis=("largest component size"))
~~~~
"""
function erdos_renyi_undirected(n::Int, p::Float64)
    if p < 0 || p > n throw(DomainError()) end

    if p >= 1. # interpret as average degree
        p = p/n # convert to probability
    end
    A = sprand(n,n,p)
    Aup = triu(A,1)
    Asym = max(Aup,Aup')
    return _matrix_network_direct(Asym,1)
end

erdős_rényi_undirected(n,p) = erdos_renyi_undirected(n,p)

"""
`erdos_renyi_directed`
========================

Generate an directed Erdős-Rényi graph.

**The current implementation uses sprand, this may change in the future.**
**Do not depend on this routine for reliable output between versions.**

Input
-----
- `n`: the number of nodes
- `p`: the probability of an edge.

Output
------
- A matrix network type for the Erdős-Rényi graph.
"""
function erdos_renyi_directed(n::Int, p::Float64)
    if p < 0 || p > n throw(DomainError()) end
    if p >= 1. # interpret as average degree
        p = p/n # convert to probability
    end

    A = sprand(n,n,p)

    return _matrix_network_direct(A-spdiagm(diag(A),0)) # directions don't matter
end

erdős_rényi_directed(n,p) = erdos_renyi_directed(n,p)

"""
`chung_lu_undirected`
========================

Generate an approximate undirected Chung-Lu graph. The approximation is because we draw
exactly |E| edges where each edge is sampled from the Chung-Lu model. But then
we discard duplicate edges and self-loops. So the new graph will always have fewer
edges than the input degree sequence.

**This will likely change in future versions and provide an exact Chung-Lu model.**

If the graph is

Usage
-----
- `chung_lu_undirected(d)`
- `chung_lu_undirected(d,nedges)`

Input
-----
- `d`: the degree sequence vector. This vector is non-negative and has the expected
degree for each vertex.

Output
------
- A MatrixNetwork for the undirected graph that results from the Chung-Lu sample.

Example
-------
~~~~
A = load_matrix_network("tapir")
d = vec(sum(A,1))
B = sparse(chung_lu_undirected(d))
nnz(A)
nnz(B)
~~~~

"""
function chung_lu_undirected end

function chung_lu_undirected(d::Vector{Int})
    chung_lu_undirected(d, floor(Int,sum(d)/2))
end

function chung_lu_undirected(d::Vector{Int}, nedges::Int)
    n = length(d)

    # TODO find some standard function for this
    for v in d
        if v < 0
            throw(DomainError())
        end
    end

    if nedges < 0 || nedges > div(n*(n-1),2)
        throw(ArgumentError("nedges $nedges is too large for $n node undirected graph"))
    end

    nodevec = zeros(Int,sum(d))
    curedge = 1
    for i=1:n
        for j=1:d[i]
            nodevec[curedge] = i
            curedge += 1
        end
    end

    ei, ej = unique_edge_sample_undirected(nodevec, nedges)
    A = sparse(ei,ej,1.,n,n)

    return _matrix_network_direct(A) # avoid the transpose
end

function unique_edge_sample_undirected(nodevec, nedges::Int)
    edges = Set{Tuple{Int,Int}}()
    sizehint!(edges, nedges)

    curedges = 0
    while curedges < nedges
        src = rand(nodevec)
        dst = rand(nodevec)
        if src == dst
            continue
        else
            # fix the order
            if src < dst
                dst, src = src, dst
            end
            if !((src,dst) in edges)
                push!(edges, (src,dst))
                curedges += 1
            end
        end
    end

    ei = zeros(Int,nedges*2)
    ej = zeros(Int,nedges*2)

    k = 1
    for edge in edges
        ei[k] = edge[1]
        ej[k] = edge[2]
        k+=1
        ei[k] = edge[2]
        ej[k] = edge[1]
        k+=1
    end

    return ei, ej
end

"""
Not public right now
"""
function _chung_lu_dense_undirected(d::Vector{Int})
    M = zeros(Int,n,n)
    idenom = 1/sum(d)
    for i=1:n
        for j=i+1:n
            if rand(Float64) < d[i]*d[j]*idenom
                M[i,j] = 1
            end
        end
    end
    M = M + M'
    return _matrix_network_direct(sparse(M))
end

"""
The internal Havel-Hakimi function has an optional store
behavior that saves the edges as they come out of the algorithm.
This enables us to generate a Havel Hakimi graph, which can
be useful.
"""
function _havel_hakimi(degs::Vector{Int}, store::Bool, ei::Vector{Int}, ej::Vector{Int})
    q = Collections.PriorityQueue(Int,Int,Base.Order.Reverse)
    n = length(degs)
    effective_n = n
    degsum = 0
    dmax = 0
    for (i,d) in enumerate(degs)
        q[i] = d
        degsum += d
        dmax = max(d,dmax)
        if d < 0
            throw(ArgumentError("the degree sequence must be non-negative"))
        end
    end

    if mod(degsum,2) != 0; return false; end
    if n > 0 && dmax >= n; return false; end # n > 0 checks for the empty graph

    if store
        resize!(ei,degsum)
        resize!(ej,degsum)
    end

    dlist = Vector{Pair{Int,Int}}(dmax)
    enum = 1

    while !isempty(q)
        vi,d = Collections.peek(q) # vi is the cur vertex, d is the cur deg
        Collections.dequeue!(q)    # remove it
        for n=1:d                  # make a list of each neighbor
            if isempty(q); return false; end
            dlist[n] = Collections.peek(q)
            Collections.dequeue!(q)
        end
        # now "add" an edge from vi->neighbor, and thus, decrease it's
        # degree when we re-add it
        for n=1:d
            if store
                ei[enum] = vi
                ej[enum] = dlist[n][1]
                enum += 1
                ei[enum] = dlist[n][1]
                ej[enum] = vi
                enum += 1
            end
            if dlist[n][2] < 1
                return false
            elseif dlist[n][2] == 1
                # don't both re-adding the vertex
            else
                q[dlist[n][1]] = dlist[n][2] - 1
            end
        end
    end
    return true
end

"""
`is_graphical_sequence`
=======================

Check whether or not a degree sequence is graphical,
which means that it is a valid degree sequence for
an undirected graph.

Note that this does not mean it is a valid degree
sequence for a connected undirected graph. So,
for instance,
`[1,1,1,1]`
is a valid degree sequence for two
disconnected edges

Usage
-----
`is_graphical_sequence(d)` returns true or false

Input
-----
- `d::Vector{Int}`:  a vector of integer valued degrees

Output
------
- a boolean that is true if the sequence is graphical
"""
function is_graphical_sequence(d::Vector{Int})
    return _havel_hakimi(d, false, Int[],Int[])
end

"""
`havel_hakimi_graph`
====================

Create a graph with a given degree sequence

Usage
-----
`A = havel_hakimi_graph(d)` returns an instance of the
a graph with degree sequence d or throws ArgumentError
if the degree sequence is not graphical.

Input
-----
- `d::Vector{Int}`:  a vector of integer valued degrees

Output
------
-`A`: a matrix network for the undirected graph that
results from the Havel-Hakimi procedure.
"""
function havel_hakimi_graph(d::Vector{Int})
    ei = Int[]
    ej = Int[]
    if _havel_hakimi(d, true, ei, ej) == false
        throw(ArgumentError("the degree sequence is not graphical"))
    end
    return MatrixNetwork(ei,ej,length(d))
end

# TODO Add chung-lu for general floating point weights
# via union of ER graphs add


text:
"""
generic estimators for the CPIMC formalism
"""
module Estimators

using ..CPIMC

import ..CPIMC: right_type_1_count, longest_type_1_chain_length, energy, w, signum

export E, Ekin, W, W_off_diag, W_diag, K, occupations, signum, longest_type_1_chain_length, right_type_1_count

"""
    Ekin(m::Model, e::Ensemble, c::Configuration)

estimator for the kinetic energy
    """
function Ekin(m::Model, e::Ensemble, c::Configuration)
    if isempty(c.kinks)
        return sum(energy(m, n) for n in c.occupations)
    end

    occs = copy(c.occupations)
    E_kin = 0
    # betrachte das komplette Intervall zwischen letztem und erstem Kink als
    # erstes intervall
    old_τ = first(last(c.kinks)) - 1
    for (τ,kink) in c.kinks
        E_kin += sum(energy(m, n) for n in occs) * float(τ - old_τ)
        old_τ = τ
        excite!(occs, kink)
    end
    return E_kin
end

"""
    W_off_diag(m::Model, e::Ensemble, c::Configuration)

estimator for the offdiagonal contribution to the interaction energy
This estimator is redundant because we can calculate the offdiagonal interaction energy from K_fermion.
"""
function W_off_diag(m::Model, e::Ensemble, c::Configuration)
    return -(length(c.kinks)/e.β)
end

"""
    W_diag(m::Model, e::Ensemble, c::Configuration)

estimator for the diagonal contribution to the interaction energy
"""
function W_diag(m::Model, e::Ensemble, c::Configuration)
    W_diag = 0

    if isempty(c.kinks)
        for occ1 in c.occupations
            redundant = true
            for occ2 in c.occupations
                if !redundant & (occ1.spin == occ2.spin)
                    W_diag += e.λ * w(m, occ1, occ2, occ2, occ1)
                end
                if occ1 == occ2
                    redundant = false
                end
            end
        end
    else
        occs = copy(c.occupations)
        old_τ = first(last(c.kinks)) - 1
        for (τ,kink) in c.kinks
            for occ1 in occs
                redundant = true
                for occ2 in occs
                    if !redundant & (occ1.spin == occ2.spin)
                        W_diag += e.λ * w(m, occ1, occ2, occ2, occ1) * (τ-old_τ)
                    end
                    if occ1 == occ2
                        redundant = false
                    end
                end
            end
            old_τ = τ
            excite!(occs, kink)
        end
    end
    W_diag *= -1
    return W_diag
end

"""
    W(m::Model, e::Ensemble, c::Configuration)

estimator for the interaction energy
This estimator is redundant because we can calculate the interaction energy from K_fermion and W_diag.
"""
function W(m::Model, e::Ensemble, c::Configuration)
    return W_diag(m, e, c) + W_off_diag(m, e,c)
end

"""
    E(m::Model, e::Ensemble, c::Configuration)

estimator for the interaction energy
This estimator is redundant because we can calculate the full energy as the sum of the kinetic and the interaction part.
"""
function E(m::Model, e::Ensemble, c::Configuration)
    return Ekin(m, e, c) + Epot(m, e, c)
end

"""
    K(m::Model, e::Ensemble, c::Configuration)

estimator for the number of kinks
"""
function K(m::Model, e::Ensemble, c::Configuration)
    return length(c.kinks)
end

"""
    occupations(m::Model, e::Ensemble, c::Configuration, emax::Int=100) :: Array{Float64,1}

estimator for the occupations of the emax:Int lowest single particle energy eigenvalues
"""
function occupations(m::Model, e::Ensemble, c::Configuration, emax::Int=100) :: Array{Float64,1}
    nk = zeros(Float64, emax)
    if isempty(c.kinks)
        for en in [ energy(m,n) for n in c.occupations ]
            if en < emax
                nk[en+1] = nk[en+1] + 1.0
            end
        end
    else
        occs = copy(c.occupations)
        old_τ = first(last(c.kinks)) - 1
        for (tau,k) in c.kinks
            for en in [ energy(m,n) for n in c.occupations ]
                if en < emax
                    nk[en+1] = nk[en+1] + float(tau - old_τ)
                end
            end
            old_τ = tau
            excite!(occs,k)
        end
    end
    nk
end

"""
    signum(m::Model, e::Ensemble, c::Configuration)

estimator for the sign of the weight function
"""
signum(m::Model, e::Ensemble, c::Configuration) = signum(m, c)

longest_type_1_chain_length(m::Model, e::Ensemble, c::Configuration) = longest_type_1_chain_length(c.kinks)

right_type_1_count(m::Model, e::Ensemble, c::Configuration) = right_type_1_count(c.kinks)

end

ext: jl


text:
########### Defines the parameters & plot for Fig 1 b,d,f,h,j ###########
bi = 1.0; # set the birth rate of infected individuals as a constant 1.0
# import required variables/modules common to all simulations
include("../Model_scripts/fig1params.jl")
include("../Model_scripts/Functions.jl")
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N';
########### Figure 1b ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "by=1.0\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_1b.png"
)
########### Figure 1h ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# plot the weighted average mortality rate through time
plot_mortality(1, size(Np)[1], "$figure_directory/graph_1h.png")
########### Figure 1f ###########
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculate weighted H0
H0_w = calculate_H0()
# calculate weighted V0
V0_w = calculate_V0()
# plot the average weighted V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_1f.png")
########### Figure 1d ###########
# plot the average weighted R0
plot_reproductive_rate(6, size(Np)[1], "$figure_directory/graph_1d.png")
########### Figure 1j ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_1j.png")


text:
module Function
import ...PatternStructure.Trees: bindings
import ...Analyzer.Function: analyze
import ...Matching.Function: matchtree, newstate
import ...Matching.Environment: Variables
import ...Helper: clean_code
import Base.Meta: quot
export destructure

function destructure(pattern, ex, body)
  ptree = analyze(pattern).child
  @gensym vars
  quote
    $vars = $Variables()
    if(!($matchtree($ptree, $ex, $newstate($vars))))
      throw(ArgumentError("$($(ptree)) cannot be matched with $($(clean_code(ex)))"))
    end

    $(code(vars, ptree, body))
  end
end

code(vars, ptree, body) =
  :(let $(declarations(vars, ptree)...); $body end)

declarations(vars, ptree) =
  map(x->declare(vars, x), bindings(ptree))

declare(vars, name) =
  :($name = $(vars)[$(quot(name))])

end


text:
using Distributions, LinearAlgebra
using Random: seed!
using AugmentedGaussianProcesses
using CSV, DataFrames
seed!(42)
N = 3000
D = 20
K = 4
k = RBFKernel(1.0)
X = rand(N,D)
y = rand(MvNormal(kernelmatrix(X,k)+1e-3I))
df = DataFrame(X)
df.y_reg = y
df.y_class = sign.(y)
width = maximum(y)-minimum(y);normy = (y.-minimum(y))/width*K
df.y_multi = floor.(Int64,normy)
CSV.write(joinpath(@__DIR__,"benchmarkdata.csv"),df)


text:
@testset "RepeatPartitionIterator" begin

    @testset "RepeatPartitionIterator basic" begin

        bitr = RepeatPartitionIterator(1:20, 5)

        for (itr, exp) in zip(bitr, [1:5, 6:10, 11:15, 16:20])
            @test collect(itr) == exp
            @test collect(itr) == exp
            @test collect(itr) == exp
        end
    end

    @testset "RepeatPartitionIterator partition" begin

        bitr = RepeatPartitionIterator(Iterators.partition(1:20, 5), 2)

        for (itr, exp) in zip(bitr, [[1:5, 6:10], [11:15, 16:20]])
            @test collect(itr) == exp
            @test collect(itr) == exp
            @test collect(itr) == exp
        end
    end

    @testset "RepeatPartitionIterator repeated partition" begin
        import IterTools: ncycle
        bitr = RepeatPartitionIterator(ncycle(Iterators.partition(1:20, 5), 3), 2)

        cnt = 0;
        for (itr, exp) in zip(bitr, [[1:5, 6:10], [11:15, 16:20],[1:5, 6:10], [11:15, 16:20],[1:5, 6:10], [11:15, 16:20]])
            @test collect(itr) == exp
            @test collect(itr) == exp
            @test collect(itr) == exp
            cnt += 1
        end
    end
end

@testset "SeedIterator" begin
    rng = MersenneTwister(123)
    testitr = SeedIterator(Iterators.map(x -> x * rand(rng, Int), ones(10)); rng=rng, seed=12)
    @test collect(testitr) == collect(testitr)

    rng = MersenneTwister(1234)
    nesteditr = SeedIterator(Iterators.map(x -> x * rand(rng, Int), testitr); rng=rng, seed=1)
    @test collect(nesteditr) == collect(nesteditr)
end

@testset "GpuIterator" begin
    dorg = 1:100;
    itr = GpuIterator(zip([view(dorg,1:10)], [dorg]))
    d1,d2 = first(itr) |> cpu
    @test !isa(d1, SubArray)
    @test d1 == dorg[1:10]
    @test d2 == dorg
end

@testset "BatchIterator" begin

    @testset "Single array" begin
        itr = BatchIterator(collect(reshape(1:2*3*4*5,2,3,4,5)), 2)
        for (i, batch) in enumerate(itr)
            @test size(batch) == (2,3,4,i==3 ? 1 : 2)
        end

        @test "biter: $itr" == "biter: BatchIterator(size=(2, 3, 4, 5), batchsize=2, shuffle=false)"
    end

    @testset "Tuple data shuffle=$shuffle" for shuffle in (true, false)
        itr = BatchIterator((collect([1:10 21:30]'), 110:10:200), 3; shuffle)
        for (i, (x, y)) in enumerate(itr)
            expsize = i == 4 ? 1 : 3
            @test size(x) == (2, expsize)
            @test size(y) == (expsize,)
        end
    end

    @testset "BatchIterator singleton" begin
        itr = BatchIterator(Singleton([1,3,5,7,9,11]), 2)
        for (i, b) in enumerate(itr)
            @test b == [1,3] .+ 4(i-1)
        end
    end

    @testset "BatchIterator shuffle basic" begin
        @test reduce(vcat, BatchIterator(1:20, 3; shuffle=true)) |> sort == 1:20

        itr = BatchIterator(ones(2,3,4), 4; shuffle=MersenneTwister(2))
        @test "siter: $itr" == "siter: BatchIterator(size=(2, 3, 4), batchsize=4, shuffle=true)"
    end

    @testset "BatchIterator shuffle ndims $(length(dims))" for dims in ((5), (3,4), (2,3,4), (2,3,4,5), (2,3,4,5,6), (2,3,4,5,6,7))
        sitr = BatchIterator(collect(reshape(1:prod(dims),dims...)), 2;shuffle=MersenneTwister(12))
        bitr = BatchIterator(collect(reshape(1:prod(dims),dims...)), 2)
        sall, nall = Set{Int}(), Set{Int}()
        for (sb, nb) in zip(sitr, bitr)
            @test sb != nb
            @test size(sb) == size(nb)
            push!(sall, sb...)
            push!(nall, nb...)
        end
        @test sall == nall
    end
end

@testset "RepeatPartitionIterator and ShuffleIterator" begin
    import IterTools: ncycle

    @testset "Single epoch small" begin
        ritr = RepeatPartitionIterator(BatchIterator(1:20, 3; shuffle=MersenneTwister(123)), 4)
        for itr in ritr
            @test collect(itr) == collect(itr)
        end
    end

    @testset "Multi epoch small" begin
        sitr = BatchIterator(1:20, 3;shuffle=MersenneTwister(123))
        citr = ncycle(sitr, 2)
        ritr = RepeatPartitionIterator(SeedIterator(citr; rng=sitr.rng), 4)
        for itr in ritr
            @test collect(itr) == collect(itr)
        end
    end

    @testset "Multi epoch big" begin
        sitr = BatchIterator(1:20, 3;shuffle= MersenneTwister(123))
        citr = ncycle(sitr, 4)
        ritr = RepeatPartitionIterator(SeedIterator(citr; rng=sitr.rng), 10)
        for (i, itr) in enumerate(ritr)
            @test collect(itr) == collect(itr)
        end
    end
end

@testset "StatefulGenerationIter" begin
    import NaiveGAflux: itergeneration, StatefulGenerationIter
    ritr = RepeatPartitionIterator(BatchIterator(1:20, 3), 4)
    sitr = RepeatPartitionIterator(BatchIterator(1:20, 3), 4) |> StatefulGenerationIter
    for (i, itr) in enumerate(ritr)
        @test collect(itr) == collect(itergeneration(sitr, i))
    end
end

@testset "TimedIterator" begin

    @testset "No stopping accumulate = $acc" for (acc, exp) in (
        (true, 7),
        (false, 0)
    )
        timeoutcnt = 0

        titer = TimedIterator(;timelimit=0.1, patience=2, timeoutaction = () -> timeoutcnt += 1, accumulate_timeouts=acc, base=1:10)

        @test collect(titer) == 1:10
        @test timeoutcnt === 0 # Or else we'll have flakey tests...

        for i in titer
            if iseven(i)
                sleep(0.11) # Does not matter here if overloaded CI VM takes longer than this to get back to us
            end
        end
        # When accumulating timeouts: after 1,2,3,4 our patience is up, call timeoutaction for 4,5,6,7,8,9,10
        # When not accumulating: We never reach patience level
        @test timeoutcnt == exp
    end

    @testset "Stop iteration at timeout" begin
        # also test that we really timeout when not accumulating here
        titer = TimedIterator(;timelimit=0.1, patience=4, timeoutaction = () -> TimedIteratorStop, accumulate_timeouts=false, base=1:10)

        last = 0
        for i in titer
            last = i
            if i > 2
                sleep(0.11) # Does not matter here if overloaded CI VM takes longer than this to get back to us
            end
        end
        @test last === 6 # Sleep after 2, then 4 patience
    end

end


text:
export train

"""
Main training driver function
"""
function train(con_path, sty_path, n_epoch = 140;
              lr = 2., noise = 0.6f0)

    info("Loading data ....")
    content = load(con_path)
    orig_size = size(content)
    content = preprocess_img(content, (224, 224)) |> gpu

    style = load(sty_path)
    style = preprocess_img(style, (224, 224)) |> gpu
    info("Loaded data")

    # Random starting point
    # generated = rand(eltype(content), size(content)) * 256
    noise = Float32(noise)
    generated = (1-noise)*content + gpu(noise*rand(-20.0f0:0.01f0:20.0f0, size(content)))
    generated = generated |> gpu
    generated = generated |> param

    m = VGG19() |> gpu

    act_content = forward_pass(m, content)
    act_style = forward_pass(m, style)

    STYLE_LAYERS =
    [(1,  0.2),
     (4,  0.2),
     (7,  0.2),
     (12, 0.2),
     (17, 0.2)]

    #optimize(x -> compute_cost(m, content, style, x, STYLE_LAYERS),
    #         generated, GradientDescent(), Optim.Options(show_trace = true, iterations = 2))
    loss() = compute_cost(m, act_content, act_style, generated, STYLE_LAYERS)

    @epochs n_epoch Flux.train!(loss,
                                Iterators.repeated((), 1), ADAM([generated], lr))

    output = generated.data |> collect |> postprocess_img
    output = imresize(output, orig_size)

    save("$(con_path)_output.jpg", output)

    output
end

function compute_cost(m, act_content, act_style, generated, STYLE_LAYERS)

    act_generated = forward_pass(m, generated)
    J_content = compute_content_cost(act_content, act_generated)
    println("Content cost = $J_content")
    J_style = compute_style_cost(act_style, act_generated, STYLE_LAYERS)
    println("Style cost = $J_style")
    J = compute_total_cost(J_content, J_style)
    println("Total cost = $J")
    println("---------------")
    J
end

function compute_content_cost(act_content, act_generated)
    _content = act_content[13]
    _generated = act_generated[13]
    cost = 0.
    #=for i = 1:endof(_content)
        cost += (_content[i] - _generated[i]) ^ 2
    end=#
    cost = sum(abs2, _content - _generated)
    cost /= (prod(size(_content)) * 4)
end

function compute_style_cost(act_style, act_generated, style_layers)
    J = 0.
    for (layer, coeff) in style_layers
        _style = act_style[layer]
        _generated = act_generated[layer]
        cost = compute_layer_style_cost(_style, _generated)
        J += coeff*cost
    end
    J
end

function compute_layer_style_cost(style, generated)

    m, n, c, _ = size(style)
    style_unrolled = reshape(style, m * n, c)
    generated_unrolled = reshape(generated, m * n, c)

    GS = gram_matrix(style_unrolled)
    GG = gram_matrix(generated_unrolled)

    cost = 0.
    #=for i = 1:endof(GS)
        cost += (GS[i] - GG[i]) ^ 2
    end=#
    cost = sum(abs2, GS - GG)

    cost /= (4 * c^2 * (m*n)^2)
    # sum(abs2, GS - GG) / (4 * c * (m*n)^2)
end

function compute_total_cost(content_cost, style_cost, α = 10, β = 40)
    α*content_cost + β*style_cost
end

function gram_matrix(A)
    #A * transpose(A)
    #A_mul_Bt(A, A)
    A'A
end

function forward_pass(m, arr)

    n = m.layers.layers |> length
    act = []

    input = arr
    for i = 1:n
        output = m.layers[i](input)
        push!(act, output)
        input = output
    end

    act
end

end


text:
using MPSKit,TensorKit,Test

th = nonsym_xxz_ham()
ts = InfiniteMPS([ℂ^3],[ℂ^48]);
(ts,pars,_) = find_groundstate(ts,th,Vumps(maxiter=400));
(energies,Bs) = quasiparticle_excitation(th,Float64(pi),ts,pars);
@test energies[1] ≈ 0.41047925 atol=1e-4


text:
include("../src/read/read_datasets.jl")
include("../src/utils.jl")
include("../src/opq/OPQ.jl");
include("../src/chainq/chainq.jl");
include("../src/lsq/LSQ.jl");
include("../src/encodings/encode_icm_cuda.jl");
include("../src/linscan/Linscan.jl");

function demo_lsq(
  dataset_name="SIFT1M",
  nread::Integer=Int(1e4)) # Increase this to 1e5 to use the full dataset

  # === Hyperparams ===
  m       = 7 # In LSQ we use m-1 codebooks
  h       = 256
  verbose = true
  nquery  = Int(1e4)
  knn     = Int(1e3) # Compute recall up to
  b       = Int( log2(h) * m )
  niter   = 10

  # === OPQ initialization ===
  x_train              = read_dataset(dataset_name, nread )
  d, _                 = size( x_train )
  C, B, R, train_error = train_opq(x_train, m, h, niter, "natural", verbose)
  @printf("Error after OPQ is %e\n", train_error[end])

  # === ChainQ initialization ===
  B                    = convert( Matrix{Int16}, B )
  C, B, R, train_error = train_chainq( x_train, m, h, R, B, C, niter )
  @printf("Error after ChainQ is %e\n", train_error[end])

  # === LSQ train ===
  ilsiter = 8
  icmiter = 4
  randord = true
  npert   = 4

  C, B, cbnorms, B_norms, obj = train_lsq( x_train, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert )
  cbnorms = vec( cbnorms[:] )

  # === Encode the base set ===
  nread_base   = Int(1e6)
  x_base       = read_dataset(dataset_name * "_base", nread_base )
  B_base       = randinit(nread_base, m, h) # initialize B at random

  ilsiter_base = 16 # LSQ-16 in the paper
  nsplits = 2 # Split the base set in two so it fits in memory. Split more if you run out of GPU memory.
  B_base, _ = encode_icm_cuda( x_base, B_base, C, [ilsiter_base], icmiter, npert, randord, nsplits, verbose )
  B_base    = B_base[end]

  base_error = qerror( x_base, B_base, C )
  @printf("Error in base is %e\n", base_error)

  # Compute and quantize the database norms
  dbnormsB = quantize_norms(B_base, C, cbnorms)
  db_norms = vec( cbnorms[dbnormsB] )

  # === Compute recall ===
  x_query = read_dataset( dataset_name * "_query", nquery, verbose )
  gt      = read_dataset( dataset_name * "_groundtruth", nquery, verbose )
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt + 1;
  end
  gt           = convert( Vector{UInt32}, gt[1,1:nquery] )
  B_base       = convert( Matrix{UInt8}, B_base-1 )

  print("Querying m=$m ... ")
  @time dists, idx = linscan_lsq( B_base, x_query, C, db_norms, eye(Float32, d), knn )
  println("done")

  idx = convert( Matrix{UInt32}, idx );
  rec = eval_recall( gt, idx, knn )

end

# train
demo_lsq()

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
module Simulation

using CSV, DataFrames, Distributions, IntervalCensoredEventTimes,
Ipopt, NLopt, Random, SparseArrays, Statistics, Test

Random.seed!(123)
# number of observations
n = 1_000
# 2 covariates in base model: one binary and one continuous
Z  = [rand(Bernoulli(0.5), n) randn(n)]
βz = [1.0, 0.8]
# auto-correlated covariates for variable selection
p  = 100
ρ  = 0.1
Σx = [ρ^(abs(i - j)) for i in 1:p, j in 1:p]
X  = Matrix(transpose(rand!(MvNormal(Σx), Matrix{Float64}(undef, p, n))))
# first 5 coefficients are signal
βx = zeros(p)
βx[1:5] .= 0.25
# systematic component
η  = Z * βz + X * βx
expη = exp.(-η)
# actual event times
T = rand.(Exponential.(expη))
# examination time points
s = quantile(T, 0.0:0.1:1.0)
s = [s[1]/2; s; s[end]+1]
display(s); println()
# censoring
C = rand([:left_censored, :right_censored, :interval_censored], n)
L = Vector{Float64}(undef, n)
R = Vector{Float64}(undef, n)
for i in 1:n
    if C[i] == :left_censored
        L[i] = 0
        R[i] = s[findfirst(sj -> sj ≥ T[i], s)]
    elseif C[i] == :right_censored
        L[i] = s[findlast(sj -> sj < T[i], s)]
        R[i] = Inf
    elseif C[i] == :interval_censored
        L[i] = s[findlast(sj -> sj < T[i], s)]
        R[i] = s[findfirst(sj -> sj ≥ T[i], s)]
    end
end
@info "Simulated interval-censored data:"
show(DataFrame(T=T, C=C, L=L, R=R)); println()
@info "Base model:"
icm = IntervalCensoredModel(Z, L, R)
solver = Ipopt.IpoptSolver(print_level=3)
@time fit!(icm, solver, init=initialize_uniform!(icm))
display(icm); println()
@info "Sparse model by IHT:"
icsm = IntervalCensoredSparseModel(icm, X, 5)
@time iht!(icsm, verbose=true, debiasing=true, printfreq=1)
βx = sparsevec(icsm.βx)
display(DataFrame(j=βx.nzind, βj=βx.nzval)); println()
end


text:
"""
```
tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; verbose = :high,
    n_particles = 1000, fixed_sched = [], r_star = 2, findroot = bisection,
    xtol = 1e-3, resampling_method = :multionial, n_mh_steps = 1, c_init = 0.3,
    target_accept_rate = 0.4, n_presample_periods = 0, allout = true,
    parallel = false, verbose = :low,
    dynamic_measurement = false, poolmodel = false)
```

### Inputs

- `data::Matrix{S}`: `Ny` x `Nt` matrix of historical data
- `Φ::Function`: state transition equation: s_t = Φ(s_{t-1}, ϵ_t)
- `Ψ::Function`: measurement equation: y_t = Ψ(s_t) + u_t
- `F_ϵ::Distribution`: shock distribution: ϵ_t ~ F_ϵ
- `F_u::Distribution`: measurement error distribution: u_t ~ F_u
- `s_init::Matrix{S}`: `Ns` x `n_particles` matrix of initial state vectors s_0

where `S<:AbstractFloat` and

- `Ns` is the number of states
- `Ny` is the number of observables
- `Nt` is the number of data periods

### Keyword Arguments

- `n_particles::Int`: number of particles used to approximate the
  log-likelihood. Using more particles yields a more accurate estimate, at the
  cost of being more computationally intensive

**Correction:**

- `fixed_sched::Vector{S}`: array of monotone increasing values in (0,1] which
  specify the fixed tempering schedule. Set to `[]` if you wish to adaptively
  choose the tempering schedule (below)
- `r_star::S`: target inefficiency ratio used to adaptively choose φ, i.e.
  Ineff(φ_n) := 1/M Σ_{j=1}^M (Wtil_t^{j,n}(φ_n))^2 = r_star
- `findroot::Function`: root-finding algorithm, one of `bisection` or `fzero`
- `xtol::S`: root-finding tolerance on x bracket size

**Selection:**

- `resampling_method::Symbol`: either `:multinomial` or `:systematic`

**Mutation:**

- `n_mh_steps::Int`: number of Metropolis-Hastings steps to take in each
  tempering stage
- `c_init::S`: initial scaling of proposal covariance matrix
- `target_accept_rate::S`: target Metropolis-Hastings acceptance rate

**Other:**

- `n_presample_periods::Int`: number of initial periods to omit from the
  log-likelihood calculation
- `allout::Bool`: whether to return all outputs or just `sum(loglh)`
- `parallel::Bool`: whether to use `SharedArray`s
- `verbose::Symbol`: amount to print to STDOUT. One of `:none`, `:low`, or
  `:high`

### Outputs

- `sum(loglh)::S`: log-likelihood p(y_{1:T}|Φ,Ψ,F_ϵ,F_u) approximation
- `loglh::Vector{S}`: vector of conditional log-likelihoods p(y_t|y_{1:t-1},Φ,Ψ,F_ϵ,F_u)
- `times::Vector{S}`: vector of runtimes per period t
"""
function tempered_particle_filter(data::AbstractArray, Φ::Function, Ψ::Function,
                                  F_ϵ::Distribution, F_u::Distribution,
                                  s_init::AbstractArray{S}; n_particles::Int = 1000,
                                  fixed_sched::Vector{S} = zeros(0), r_star::S = 2.0,
                                  findroot::Function = bisection, xtol::S = 1e-3,
                                  resampling_method = :multinomial, n_mh_steps::Int = 1,
                                  c_init::S = 0.3, target_accept_rate::S = 0.4,
                                  n_presample_periods::Int = 0, allout::Bool = true,
                                  parallel::Bool = false, get_t_particle_dist::Bool = false,
                                  verbose::Symbol = :high,
                                  dynamic_measurement::Bool = false,
                                  poolmodel::Bool = false) where S<:AbstractFloat
    #--------------------------------------------------------------
    # Setup
    #--------------------------------------------------------------

    # If using fixed φ schedule, check well-formed
    adaptive_φ = isempty(fixed_sched)
    if !adaptive_φ
        @assert fixed_sched[end] == 1 "φ schedule must be a range from [a,1] s.t. a > 0."
    end

    # Initialize constants
    n_obs, T  = size(data)
    n_shocks  = length(F_ϵ)
    n_states  = size(s_init, 1)
    QQerr = false
    HHerr = false
    try
        QQ = cov(F_ϵ)
    catch
        QQerr = true
    end
    try
        HH = cov(F_u)
    catch
        HHerr = true
    end
    if QQerr
        QQ = F_ϵ.σ * ones(1,1)
    end
    if HHerr
        HH = zeros(1,1)
    end
    @assert @isdefined HH

    # Initialize output vectors
    loglh = zeros(T)
    times = zeros(T)

    # Initialize matrix of normalized weight per particle by time period
    # and a matrix of the particle locations, if desired by user
    if get_t_particle_dist
        t_norm_weights = Matrix{Float64}(undef,n_particles,T)
        t_particle_dist = Dict{Int64,Matrix{Float64}}()
    end

    # Initialize working variables
    s_t1_temp     = parallel ? SharedMatrix{Float64}(copy(s_init)) :
        Matrix{Float64}(copy(s_init))
    s_t_nontemp   = parallel ? SharedMatrix{Float64}(n_states, n_particles) :
        Matrix{Float64}(undef, n_states, n_particles)
    ϵ_t           = parallel ? SharedMatrix{Float64}(n_shocks, n_particles) :
        Matrix{Float64}(undef, n_shocks, n_particles)
    coeff_terms   = parallel ? SharedVector{Float64}(n_particles) :
        Vector{Float64}(undef, n_particles)
    log_e_1_terms = parallel ? SharedVector{Float64}(n_particles) :
        Vector{Float64}(undef, n_particles)
    log_e_2_terms = parallel ? SharedVector{Float64}(n_particles) :
        Vector{Float64}(undef, n_particles)

    inc_weights   = Vector{Float64}(undef, n_particles)
    norm_weights  = Vector{Float64}(undef, n_particles)

    c = c_init
    accept_rate = target_accept_rate

    # If not using a dynamic measurement equation, then define measurement equation
    # applying to all states, even if they are missing (but not time variables)
    if !dynamic_measurement
        Ψ_allstates = Ψ
    end

    #--------------------------------------------------------------
    # Main Algorithm: Tempered Particle Filter
    #--------------------------------------------------------------

    for t = 1:T
        begin_time = time_ns()
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("============================================================")
            @show t
        end

        #--------------------------------------------------------------
        # Initialization
        #--------------------------------------------------------------

        # Remove rows/columns of series with NaN values
        y_t = data[:, t]

        # Handle measurement equation
        if !(poolmodel || dynamic_measurement)
            Ψ_t  = x -> Ψ(x)[nonmissing]
        elseif poolmodel && dynamic_measurement
            Ψ_t = x -> Ψ(x,y_t,t)[nonmissing]
            Ψ_allstates = x -> Ψ(x,y_t,t)
        elseif poolmodel
            Ψ_t = x -> Ψ(x,y_t)[nonmissing]
            Ψ_allstates = x -> Ψ(x,y_t)
        else
            Ψ_t  = x -> Ψ(x,t)[nonmissing]
            Ψ_allstates = x -> Ψ(x,t)
        end

        # Adjust other values to remove rows/columns with NaN values
        nonmissing = isfinite.(y_t)
        y_t        = y_t[nonmissing]
        n_obs_t    = length(y_t)
        HH_t     = poolmodel ? HH : HH[nonmissing, nonmissing] # poolmodel -> keep missing is ok
        inv_HH_t = poolmodel ? zeros(1,1) : inv(HH_t) # poolmodel -> don't need inv_HH
        det_HH_t = poolmodel ? 0. : det(HH_t) # poolmodel -> don't need det_HH

        # Initialize s_t_nontemp and ϵ_t for this period
        if parallel
            @sync @distributed for i in 1:n_particles
                ϵ_t[:, i] .= rand(F_ϵ)
                s_t_nontemp[:, i] = Φ(s_t1_temp[:, i], ϵ_t[:, i])
            end
        else
             for i in 1:n_particles
                ϵ_t[:, i] .= rand(F_ϵ)
                s_t_nontemp[:, i] = Φ(s_t1_temp[:, i], ϵ_t[:, i])
            end
        end

        # Tempering initialization
        φ_old = 1e-30
        stage = 0

        #--------------------------------------------------------------
        # Main Algorithm
        #--------------------------------------------------------------
        while φ_old < 1
            stage += 1

            ### 1. Correction
            # Modifies coeff_terms, log_e_1_terms, log_e_2_terms
            weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old,
                           Ψ_allstates, y_t, s_t_nontemp, det_HH_t, inv_HH_t;
                           initialize = stage == 1, parallel = parallel,
                           poolmodel = poolmodel)

            φ_new = next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, n_obs_t,
                           r_star, stage; fixed_sched = fixed_sched, findroot = findroot,
                           xtol = xtol)

            if VERBOSITY[verbose] >= VERBOSITY[:high]
                @show φ_new
            end

            # Modifies inc_weights, norm_weights
            correction!(inc_weights, norm_weights, φ_new, coeff_terms,
                        log_e_1_terms, log_e_2_terms, n_obs_t)

            ### 2. Selection
            # Modifies s_t1_temp, s_t_nontemp, ϵ_t
            selection!(norm_weights, s_t1_temp, s_t_nontemp, ϵ_t;
                       resampling_method = resampling_method)

            loglh[t] += log(mean(inc_weights))

            c = update_c(c, accept_rate, target_accept_rate)
            if VERBOSITY[verbose] >= VERBOSITY[:high]
                @show c
                println("------------------------------")
            end

            ### 3. Mutation
            # Modifies s_t_nontemp, ϵ_t
            if stage != 1
                accept_rate = mutation!(Φ, Ψ_t, QQ, det_HH_t, inv_HH_t, φ_new, y_t,
                                        s_t_nontemp, s_t1_temp, ϵ_t, c, n_mh_steps;
                                        parallel = parallel,
                                        poolmodel = poolmodel)
            end

            φ_old = φ_new
        end # of loop over stages

        if get_t_particle_dist
            # save the normalized weights in the column for period t
            t_norm_weights[:,t] = norm_weights
            t_particle_dist[t] = copy(s_t_nontemp)
        end

        times[t] = time_ns() - begin_time
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            print("\n")
            @show loglh[t]
            print("Completion of one period $times[t]")
        end
        s_t1_temp .= s_t_nontemp
    end # of loop over periods

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("=============================================")
    end

    if get_t_particle_dist && allout
        return sum(loglh[n_presample_periods + 1:end]), loglh[n_presample_periods + 1:end], times, t_particle_dist, t_norm_weights
    elseif get_t_particle_dist
        return sum(loglh[n_presample_periods + 1:end]), t_particle_dist, t_norm_weights
    elseif allout
        return sum(loglh[n_presample_periods + 1:end]), loglh[n_presample_periods + 1:end], times
    else
        return sum(loglh[n_presample_periods + 1:end])
    end
end


text:
zero(a::GroupElem) = zero(parent(a))
one(a::GroupElem) = one(parent(a))

# http://sagenb.org/src/groups/generic.py
multiple(a::GroupElem, n::Integer, op) = begin
    if op === +
        return multiple(a, n, +, x -> -x, zero(a))
    elseif op === *
        return multiple(a, n, *, inv, one(a))
    else
        throw(MethodError())
    end
end

multiple(a::T, n::Integer, op, inverse, identity::T) where {T<:GroupElem} = begin
    n == 0 && return identity

    if n < 0
        n = -n
        a = inverse(a)
    end

    n == 1 && return a

    # check for idempotence, and store the result otherwise
    aa = op(a, a)
    aa == a && return a

    n == 2 && return aa
    n == 3 && return op(aa, a)
    n == 4 && return op(aa, aa)

    # since we've computed a^2, let's start squaring there
    # so, let's keep the least-significant bit around, just
    # in case.
    m = n & 1
    n = n >> 1

    # One multiplication can be saved by starting with
    # the second-smallest power needed rather than with 1
    # we've already squared a, so let's start there.
    apow = aa
    while n & 1 == 0
        apow = op(apow, apow)
        n = n >> 1
    end
    power = apow
    n = n >> 1

    # now multiply that least-significant bit in...
    m == 1 && (power = op(power, a))

    # and this is straight from the book.
    while n != 0
        apow = op(apow, apow)
        n & 1 != 0 && (power = op(power, apow))
        n = n >> 1
    end

    power
end

# TODO: improve! => http://sagenb.org/src/groups/generic.py
log(a::T, base::T, ord::fmpz) where {T<:GroupElem} = log(a, base, BigInt(ord))
log(a::T, base::T, ord::Integer) where {T<:GroupElem} = begin
    n = 1
    b = base
    while a != b
        n += 1
        b *= base
    end
    n
end


text:
module Pomodoro

using Uhlarm, ProgressMeter

include("functions/pomodoro_stuff.jl")
export pomodoro_stuff

end

max_stars_repo_name: marcelovmaciel/Pomodoro

max_issues_repo_name: marcelovmaciel/Pomodoro

max_forks_repo_name: marcelovmaciel/Pomodoro


text:
using CSV, Feather
using JDF, FileIO, Blosc, StatsPlots, RCall
using DataFrames, WeakRefStrings # required for JLD2, JDF
using Fread:fread
using PyCall
Blosc.set_num_threads(6)

dirpath = "C:/data/"
largest_file = "feature_matrix_cleaned.csv"
outpath = "c:/data/jdf-bench/feature_matrix_cleaned.csv"
data_label =  "Mortgage Risk with featuretools"
delim = ','
header = true

gen_benchmark3(dirpath, largest_file, outpath, data_label; delim = ',', header=true) = begin
    if !isdir(outpath)
        mkpath(outpath)
    end

    println("CSV.read")
    csv_read1 = @elapsed df = CSV.read(joinpath(dirpath, largest_file), delim = delim, header = header);
    csv_read2 = @elapsed df = CSV.read(joinpath(dirpath, largest_file), delim = delim, header = header);

    csv_write1 = 0
    csv_write2 = 0
    try
        csv_write1 = @elapsed CSV.write(joinpath(outpath, largest_file*".csv"), df);
        csv_write2 = @elapsed CSV.write(joinpath(outpath, largest_file*".csv"), df);
    catch err
    end

    delims = string(delim)
    @rput delims
    @rput header
    R"""
    library(data.table)
    library(fst)
    # memory.limit(Inf)
    readr_read1 = system.time(a <- readr::read_delim($(joinpath(dirpath, largest_file)), delim = delims, progress=FALSE, col_names=header))[3]
    # readr_read2 = system.time(a <- readr::read_delim($(joinpath(dirpath, largest_file)), delim = delims, progress=FALSE, col_names=header))[3]
    # rbase_read1 = system.time(a <- readr::read_csv($(joinpath(dirpath, largest_file)), progress=FALSE))[3]
    # rbase_read2 = system.time(a <- readr::read_csv($(joinpath(dirpath, largest_file)), progress=FALSE))[3]
    data_table_read1 = system.time(a <- data.table::fread($(joinpath(dirpath, largest_file))))[3]
    # data_table_read2 = system.time(a <- data.table::fread($(joinpath(dirpath, largest_file))))[3]
    rm(a)
    gc()
    """

    py"""
    2+2
    """

    @rget data_table_read1
    @rget data_table_read2
    # @rget rbase_read1
    # @rget rbase_read2
    @rget readr_read1
    @rget readr_read2

    fread_read1 = 0
    fread_read2 = 0
    fread_read1 = @elapsed fread(joinpath(dirpath, largest_file))
    fread_read2 = @elapsed fread(joinpath(dirpath, largest_file))

    #read_perf  = [ csv_read1, csv_read2, fread_read1, fread_read2, data_table_read1, data_table_read2, readr_read1, readr_read2, rbase_read1, rbase_read2]
    read_perf  = [ csv_read1, csv_read2, fread_read1, fread_read2, data_table_read1, data_table_read2, readr_read1, readr_read2]

    read_pkgs = repeat(["CSV.jl", "Fread.jl", "data.table", "readr"], inner = 2)
    read_run_group = repeat(["1st", "2nd"], outer = 4)

    plot_read = groupedbar(
        read_pkgs,
        read_perf,
        group = read_run_group,
        ylab = "Seconds",
        title = "Disk-format Read performance comparison \n Data: $data_label data \n Size: $(size(df)) filesize:$(round(filesize(joinpath(dirpath, largest_file))/1024^3, digits=1))GB \n Julia $(VERSION)"
    )
    savefig(plot_read, joinpath(outpath, largest_file*"plot3_read.png"))

    (read_perf, dirpath, outpath, largest_file, df)
end

#@time res = gen_benchmark3("c:/data/AirOnTimeCSV/", "airOT199302.csv", "c:/data/jdf-bench/airOT199302.csv", "Air On Time 199302")

dirpath = "C:/data/"
largest_file = "feature_matrix_cleaned.csv"
outpath = "c:/data/jdf-bench/feature_matrix_cleaned.csv"
data_label =  "Mortgage Risk with featuretools"

@time res = gen_benchmark3(dirpath, largest_file, outpath, data_label);
#
# dirpath = "C:/data/perf/"
# largest_file = "Performance_2004Q3.txt"
# outpath = "c:/data/jdf-bench/Performance_2004Q3.txt"
# data_label =  "Fannie Mae Performance 2004Q3"
# #
# @time res = gen_benchmark3(dirpath, largest_file, outpath, data_label, delim = '|', header = false);


text:
@testset "Sparse (semi-)nonnegative Tucker decomposition" begin
Random.seed!(12345)

@testset "nonnegative decomposition" begin
    # An example of nonnegative and semi-nonnegative Tucker decomposition
    tucker_orig = rand_tucker((4, 5, 6), (40, 50, 60), factors_nonneg=true, core_nonneg=true)
    tnsr_orig = compose(tucker_orig)

    tnsr_max = maximum(tnsr_orig)
    tucker_orig.core ./= tnsr_max
    tnsr_orig ./= tnsr_max

    tnsr = add_noise(tnsr_orig, 0.6, true)

    # Solve the problem
    @time tucker_spnn = spnntucker(tnsr, size(tucker_orig.core), tol=1E-5, ini_decomp=:hosvd,
                            core_nonneg=true,
                            max_iter=1000, verbose=true, lambdas=fill(0.1, 4))

    # Reporting
    @test rel_residue(tucker_spnn) < 0.05
    @info "Relative error of decomposition : $(rel_residue(tucker_spnn))"
end

@testset "semi-nonnegative decomposition" begin
    tucker_orig = rand_tucker((4, 5, 6), (40, 50, 60), factors_nonneg=true, core_nonneg=false)
    tnsr_orig = compose(tucker_orig)

    tnsr_max = maximum(tnsr_orig)
    tucker_orig.core ./=  tnsr_max
    tnsr_orig ./= tnsr_max

    tnsr = add_noise(tnsr_orig, 0.6, false)

    # Solve the problem
    @time tucker_spnn = spnntucker(tnsr, size(tucker_orig.core), tol=1E-5, ini_decomp=:hosvd,
                             core_nonneg=false,
                             max_iter=1000, verbose=true, lambdas=fill(0.1, 4))

    # Reporting
    @test rel_residue(tucker_spnn) < 0.05
    @info("Relative error of decomposition : $(rel_residue(tucker_spnn))")
end

end


text:
mutable struct Variable
    varid::Int # 0 is target, 1 is One()
    name::Symbol
    parentfuncs::Vector{Int}#id of Func creating it.
    useids::Vector{Int}#ids of Funcs that use it.
    paireddeps::Vector{Int}
    tracked::Bool
    lowered_count::Int
    # initialized::Bool
    ref::Any
    # function Variable(name::Symbol, id::Int)
        # new(id, name, 0, Int[], false, true)
    # end
    function Variable(name::Symbol, id::Int)
        new(id, name, Int[], Int[], Int[], 0, true)
    end
end
# Base.ndims(d::Dimensions) = length(d.sizehints)
# Base.ndims(v::Variable) = length(v.dims)
# isscalar(v::Variable) = iszero(ndims(v))
istracked(v::Variable) = v.tracked
num_lowered(v::Variable) = v.lowered_count
isinitialized(v::Variable) = num_lowered(v) == length(parents(v))
LoopVectorization.name(v::Variable) = v.name
identifier(v::Variable) = v.varid

LoopVectorization.parents(v::Variable) = v.parentfuncs
hasparent(v::Variable) = !iszero(parent(v))

isref(v::Variable) = isdefined(v, :ref)

function Base.push!(x::Vector{Any}, v::Variable)
    isref(v) ? push!(x, v.ref) : push!(x, v.name)
end
# function Base.pushfirst!(x::Vector{Any}, v::Variable)
#     isref(v) ? pushfirst!(x, v.ref) : pushfirst!(x, v.name)
# end

# function Base.hash(v::Variable, u::UInt)

# end
# function Base.isequal(v1::Variable, v2::Variable)
#     v1.varid == v2.varid && return true
#     v1.name == v2.name && return true

# end

# const var"#TARGET#" = Variable(Symbol("##TARGET##"), Dimensions(Int[],Bool[]),true,0)
# const var"#TARGET#" = Variable(Symbol("##TARGET##"), 0, true)
# const var"#ONE#" = Variable(Symbol("##ONE##"), 0, true)


text:
module NewsLookout

using ArgParse
using Dates
using Logging

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--conf", "-c"
            help = "Configuration file name"
			arg_type = String
			default = "conf/file1.conf"
        "--rundate", "-d"
            help = "Run-date to extract data"
            arg_type = String
            default = Dates.format(Dates.now(), "yyyy-mm-dd")
    end

    return parse_args(s)
end

function set_logger(io::IOStream)
	logger = Logging.SimpleLogger(io)
	with_logger(logger) do
           @info("Started logging.")
    end
	flush(io)
	Logging.global_logger(logger)
end

function greet()
	return("-- The NewsLookout Web Scraping App ---")
end

function print_usage()
	println("Usage: newslookout.jl --conf <configfile> -rundate <run date>")
end

function julia_main()
	println(greet())
	# start logger
	io = open("logfile.txt", "w+")
	set_logger(io)

    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("\t$arg  =>  $val")
    end

	# close logger
	close(io)
end

function __init__()

end

if !isdefined(Base, :active_repl)
    julia_main()
end

end # module


text:
"""
    CauchyEndgame(x::AbstractVector; samples_per_loop = 8, max_winding_number = 12)

The Cauchy endgame tries to to predict the value a solution path `x(t)` at `0` , i.e. `x(0)`,
by using a generalization of [Cauchy's integral formula]. This is derived in [^MSW90].
The general idea is that we make loops around the origin until we have again the start point.
The number of loops is referred to as the *winding* or *cycle number*.

[^MSW90]: Morgan, A.P., Sommese, A.J. & Wampler, C.W. Numer. Math. (1990) 58: 669. https://doi.org/10.1007/BF01385648
[Cauchy's integral formula]: https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula
"""
struct CauchyEndgame{AV<:AbstractVector}
    base_point::AV
    samples_per_loop::Int
    max_winding_number::Int

    function CauchyEndgame(
        x::AV;
        samples_per_loop::Int = 8,
        max_winding_number::Int = 12,
    ) where {AV<:AbstractVector}
        new{AV}(copy(x), samples_per_loop, max_winding_number)
    end
end

"""
    CauchyEndgameResult

An enum indicating the result of the [`predict!`](@ref) computation.

# Cases
* `CAUCHY_SUCCESS`: The endgame was successfull.
* `CAUCHY_TERMINATED_MAX_WINDING_NUMBER`: The endgame was terminated since the winding
  number is larger than the provided threshold.
* `CAUCHY_TERMINATED_ACCURACY_LIMIT`: The endgame was terminated since during the loop
  tracking the achievable accuracy was no more sufficient.
* `CAUCHY_TERMINATED`: The endgame was terminated due to some other error in the path
  tracking.
"""
@enum CauchyEndgameResult begin
    CAUCHY_SUCCESS
    CAUCHY_TERMINATED_MAX_WINDING_NUMBER
    CAUCHY_TERMINATED_ACCURACY_LIMIT
    CAUCHY_TERMINATED_MAX_ITERS
    CAUCHY_TERMINATED
end

"""
    predict!(p, tracker::CoreTracker, ::CauchyEndgame)

Try to predict the value of `x(0)` using the [`CauchyEndgame`](@ref).
For this we track the polygon defined by ``te^{i2πk/n}`` until we end again at ``x``.
Here ``n`` is the number of samples we take per loop, `samples_per_loop`.
The computation gives up if we have a winding number larger than `max_winding_number`.
It returns a tuple denoting the success ([`CauchyEndgameResult`](@ref)) the computed
winding number `m::Int` and th expected accuracy of the solution.

[Cauchy's integral formula]: https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula
"""
function predict!(p, tracker::CoreTracker, cauchy::CauchyEndgame)
    @unpack base_point, samples_per_loop, max_winding_number = cauchy

    tracker.options.logarithmic_time_scale || throw(ArgumentError("The `CoreTracker` needs to use a logarithmic time scale for the cauchy endgame."))

    s = current_t(tracker)
    s_target = tracker.state.segment.target

    # store the current value in p
    base_point .= current_x(tracker)
    p .= 0.0
    x = current_x(tracker)

    # during the loop we fix the affine patch
    fix_patch!(tracker)

    m = k = 1
    Δθ = 2π / samples_per_loop
    result = CAUCHY_TERMINATED_MAX_WINDING_NUMBER
    # The accuracy of the Cauchy endgame depends on the minimal accuracy
    # of a sample point. If all sample points are computed with an accuracy τ
    # we can expect τ as an error of the solution
    max_acc = tracker.state.accuracy
    while m ≤ max_winding_number
        θⱼ = 0.0
        for j = 1:samples_per_loop
            θⱼ₋₁ = θⱼ
            θⱼ += Δθ

            retcode = track!(tracker, x, s + im * θⱼ₋₁, s + im * θⱼ; loop = true)
            max_acc = max(max_acc, tracker.state.accuracy)

            if !is_success(retcode)
                init!(tracker, base_point, s, s_target; keep_steps = true)
                if retcode == CoreTrackerStatus.terminated_accuracy_limit
                    return CAUCHY_TERMINATED_ACCURACY_LIMIT, m, max_acc
                elseif retcode == CoreTrackerStatus.terminated_maximal_iterations
                    return CAUCHY_TERMINATED_MAX_ITERS, m, max_acc
                else
                    return CAUCHY_TERMINATED, m, max_acc
                end
            end

            p .+= x
        end
        # Check that loop is closed
        if norm(tracker)(base_point, x) < 4 * tracker.options.accuracy
            p ./= m .* samples_per_loop
            result = CAUCHY_SUCCESS
            break
        end

        m += 1
    end

    # we have to undo the fixing of the patch
    unfix_patch!(tracker)
    init!(tracker, base_point, s, s_target; loop = true, keep_steps = true)

    result, m, max_acc
end

fix_patch!(tracker::CoreTracker) = tracker.options.update_patch = false
unfix_patch!(tracker::CoreTracker) = tracker.options.update_patch = true


text:
using Logging
using JuMP, Gurobi
using DecisionProgramming

@info("Creating the influence diagram.")
diagram = InfluenceDiagram()

add_node!(diagram, ChanceNode("O", [], ["lemon", "peach"]))
add_node!(diagram, DecisionNode("T", [], ["no test", "test"]))
add_node!(diagram, ChanceNode("R", ["O", "T"], ["no test", "lemon", "peach"]))
add_node!(diagram, DecisionNode("A", ["R"], ["buy without guarantee", "buy with guarantee", "don't buy"]))

add_node!(diagram, ValueNode("V1", ["T"]))
add_node!(diagram, ValueNode("V2", ["A"]))
add_node!(diagram, ValueNode("V3", ["O", "A"]))

generate_arcs!(diagram)

X_O = ProbabilityMatrix(diagram, "O")
X_O["peach"] = 0.8
X_O["lemon"] = 0.2
add_probabilities!(diagram, "O", X_O)

X_R = ProbabilityMatrix(diagram, "R")
X_R["lemon", "no test", :] = [1,0,0]
X_R["lemon", "test", :] = [0,1,0]
X_R["peach", "no test", :] = [1,0,0]
X_R["peach", "test", :] = [0,0,1]
add_probabilities!(diagram, "R", X_R)

Y_V1 = UtilityMatrix(diagram, "V1")
Y_V1["test"] = -25
Y_V1["no test"] = 0
add_utilities!(diagram, "V1", Y_V1)

Y_V2 = UtilityMatrix(diagram, "V2")
Y_V2["buy without guarantee"] = 100
Y_V2["buy with guarantee"] = 40
Y_V2["don't buy"] = 0
add_utilities!(diagram, "V2", Y_V2)

Y_V3 = UtilityMatrix(diagram, "V3")
Y_V3["lemon", "buy without guarantee"] = -200
Y_V3["lemon", "buy with guarantee"] = 0
Y_V3["lemon", "don't buy"] = 0
Y_V3["peach", :] = [-40, -20, 0]
add_utilities!(diagram, "V3", Y_V3)

generate_diagram!(diagram)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, diagram)
x_s = PathCompatibilityVariables(model, diagram, z)
EV = expected_value(model, diagram, x_s)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)
S_probabilities = StateProbabilities(diagram, Z)
U_distribution = UtilityDistribution(diagram, Z)

@info("Printing decision strategy:")
print_decision_strategy(diagram, Z, S_probabilities)

@info("Printing utility distribution.")
print_utility_distribution(U_distribution)

@info("Printing expected utility.")
print_statistics(U_distribution)


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Gherkin
using Behavior.Gherkin: ScenarioStep, Background
using Behavior
using Behavior: StepDefinitionContext, StepDefinition, StepDefinitionLocation, StepDefinitionMatch
using Behavior: Executor, StepExecutionResult, QuietRealTimePresenter, executefeature
import Behavior: present

successful_step_definition(::StepDefinitionContext, args) = Behavior.SuccessfulStepExecution()
failed_step_definition(::StepDefinitionContext, args) = Behavior.StepFailed("")
error_step_definition(::StepDefinitionContext, args) = error("Some error")

struct FakeStepDefinitionMatcher <: Behavior.StepDefinitionMatcher
    steps::Dict{Behavior.Gherkin.ScenarioStep, Function}
end

function Behavior.findstepdefinition(s::FakeStepDefinitionMatcher, step::Behavior.Gherkin.ScenarioStep)
    if step in keys(s.steps)
        StepDefinitionMatch(StepDefinition("some text", s.steps[step], StepDefinitionLocation("", 0)))
    else
        throw(Behavior.NoMatchingStepDefinition())
    end
end

struct ThrowingStepDefinitionMatcher <: Behavior.StepDefinitionMatcher
    ex::Exception
end

Behavior.findstepdefinition(matcher::ThrowingStepDefinitionMatcher, ::Behavior.Gherkin.ScenarioStep) = throw(matcher.ex)

@testset "Executor             " begin
    @testset "Execute a one-step scenario; No matching step found; Result is NoStepDefinitionFound" begin
        stepdefmatcher = ThrowingStepDefinitionMatcher(Behavior.NoMatchingStepDefinition())
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[Given("some precondition")])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.NoStepDefinitionFound)
    end

    @testset "Execute a one-step scenario; The matching step is successful; Result is Successful" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
    end

    @testset "Execute a one-step scenario; The matching step fails; Result is Failed" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.StepFailed)
    end

    @testset "Execute a one-step scenario; The matching step throws an error; Result is Error" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => error_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.UnexpectedStepError)
    end

    @testset "Execute a two-step scenario; First step throws an error; Second step is Skipped" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => error_step_definition,
                                                        when => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given, when])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[2], Behavior.SkippedStep)
    end

    @testset "Execute a two-step scenario; First step fails; Second step is Skipped" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                        when => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given, when])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[2], Behavior.SkippedStep)
    end

    @testset "Execute a two-step scenario; Both steps succeed; All results are Success" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition,
                                                        when => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given, when])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[2], Behavior.SuccessfulStepExecution)
    end

    @testset "Execute a three-step scenario; All steps succeeed; All results are Success" begin
        given = Given("Some precondition")
        when = When("some action")
        then = Then("some postcondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition,
                                                        when => successful_step_definition,
                                                        then => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[given, when, then])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[2], Behavior.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[3], Behavior.SuccessfulStepExecution)
    end

    @testset "Execute a scenario; Scenario is provided; Scenario is returned with the result" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("This is a scenario", String[], ScenarioStep[given])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test scenarioresult.scenario == scenario
    end

    @testset "Execute a scenario; No unique step definition found; Result is NonUniqueMatch" begin
        stepdefmatcher = ThrowingStepDefinitionMatcher(Behavior.NonUniqueStepDefinition([]))
        executor = Behavior.Executor(stepdefmatcher)
        scenario = Scenario("Description", String[], ScenarioStep[Given("some precondition")])

        scenarioresult = Behavior.executescenario(executor, Background(), scenario)

        @test isa(scenarioresult.steps[1], Behavior.NonUniqueMatch)
    end

    @testset "Execute a ScenarioOutline; Outline has two examples; Two scenarios are returned" begin
        outline = ScenarioOutline("", String[], ScenarioStep[Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)

        outlineresult = Behavior.executescenario(executor, Background(), outline)

        @test length(outlineresult) == 2
    end

    @testset "Execute a ScenarioOutline; Outline has a successful and a failing example; First is success, second is fail" begin
        outline = ScenarioOutline("", String[], ScenarioStep[Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => failed_step_definition))
        executor = Behavior.Executor(stepdefmatcher)

        outlineresult = Behavior.executescenario(executor, Background(), outline)

        @test outlineresult[1].steps[1] isa Behavior.SuccessfulStepExecution
        @test outlineresult[2].steps[1] isa Behavior.StepFailed
    end

    @testset "Execute a ScenarioOutline; Outline has three examples; Three scenarios are returned" begin
        outline = ScenarioOutline("", String[], ScenarioStep[Given("step <stepnumber>")], ["stepnumber"], ["1" "2" "3"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        step3 = Given("step 3")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => successful_step_definition,
                                                        step3 => successful_step_definition))
        executor = Behavior.Executor(stepdefmatcher)

        outlineresult = Behavior.executescenario(executor, Background(), outline)

        @test length(outlineresult) == 3
    end

    @testset "Block text" begin
        @testset "Scenario step has a block text; Context contains the block text" begin
            given = Given("Some precondition", block_text="Some block text")
            function check_block_text_step_definition(context::StepDefinitionContext, _args)
                if context[:block_text] == "Some block text"
                    Behavior.SuccessfulStepExecution()
                else
                    Behavior.StepFailed("")
                end
            end
            stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => check_block_text_step_definition))
            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], ScenarioStep[given])

            scenarioresult = Behavior.executescenario(executor, Background(), scenario)

            @test isa(scenarioresult.steps[1], Behavior.SuccessfulStepExecution)
        end

        @testset "First step has block text, but second doesn't; The block text is cleared in the second step" begin
            given = Given("Some precondition", block_text="Some block text")
            when = When("some action")
            function check_block_text_step_definition(context::StepDefinitionContext, _args)
                if context[:block_text] == ""
                    Behavior.SuccessfulStepExecution()
                else
                    Behavior.StepFailed("")
                end
            end
            stepdefmatcher = FakeStepDefinitionMatcher(Dict(
                given => successful_step_definition,
                when => check_block_text_step_definition))
            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], [given, when])

            scenarioresult = Behavior.executescenario(executor, Background(), scenario)

            @test isa(scenarioresult.steps[2], Behavior.SuccessfulStepExecution)
        end
    end

    @testset "Backgrounds" begin
        @testset "Execute a one-step Background; No matching step found; Result is NoStepDefinitionFound" begin
            given = Given("some precondition")
            stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))

            background = Background("A background", ScenarioStep[Given("some background precondition")])

            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], ScenarioStep[Given("some precondition")])

            scenarioresult = Behavior.executescenario(executor, background, scenario)

            @test isa(scenarioresult.backgroundresult[1], Behavior.NoStepDefinitionFound)
        end

        @testset "Execute a one-step Background; A successful match found; Background result is Success" begin
            given = Given("some precondition")
            bgiven = Given("some background precondition")
            stepdefmatcher = FakeStepDefinitionMatcher(
                Dict(
                    given => successful_step_definition,
                    bgiven => successful_step_definition,
                ))

            background = Background("A background", ScenarioStep[Given("some background precondition")])

            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], ScenarioStep[Given("some precondition")])

            scenarioresult = Behavior.executescenario(executor, background, scenario)

            @test isa(scenarioresult.backgroundresult[1], Behavior.SuccessfulStepExecution)
        end

        @testset "Execute a one-step background; The matching step fails; Result is Failed" begin
            given = Given("Some precondition")
            bgiven = Given("some background precondition")
            stepdefmatcher = FakeStepDefinitionMatcher(
                Dict(
                    given => successful_step_definition,
                    bgiven => failed_step_definition,
                ))
            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], ScenarioStep[given])
            background = Background("Background description", ScenarioStep[bgiven])

            scenarioresult = Behavior.executescenario(executor, background, scenario)

            @test isa(scenarioresult.backgroundresult[1], Behavior.StepFailed)
        end

        @testset "Execute a one-step background; The background step fails; The Scenario step is skipped" begin
            given = Given("Some precondition")
            bgiven = Given("some background precondition")
            stepdefmatcher = FakeStepDefinitionMatcher(
                Dict(
                    given => successful_step_definition,
                    bgiven => failed_step_definition,
                ))
            executor = Behavior.Executor(stepdefmatcher)
            scenario = Scenario("Description", String[], ScenarioStep[given])
            background = Background("Background description", ScenarioStep[bgiven])

            scenarioresult = Behavior.executescenario(executor, background, scenario)

            @test isa(scenarioresult.steps[1], Behavior.SkippedStep)
        end

        # TODO Additional tests to ensure that Backgrounds work as any Scenario section does
    end
end


text:
@testset "approximate_inference" begin
    @testset "gaussian_reconstruction_term" begin
        _, _, Σ, y = generate_synthetic_problem(MersenneTwister(123456))
        mq = randn(size(y))
        Σq = Diagonal(rand(length(y)) .+ 1)
        gaussian_reconstruction_term(y, Σ, mq, Σq)
        Zygote.gradient(gaussian_reconstruction_term, y, Σ, mq, Σq)
    end
    @testset "approx_posterior" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))

        # Compute natural pseudo-observations and approximate posterior.
        η1, η2 = natural_from_canonical(y, Σ)
        f_post_approx = approx_posterior(f, x, η1, η2)

        # Compute the exact posterior.
        f_post_exact = posterior(f(x, Σ), y)

        # The exact and approximate posteriors should be identical in this case.
        # Check that this is the case by comparing the statistics at some test points.
        x_test = range(-7.5, 7.5; length=101)
        ms_exact = marginals(f_post_exact(x_test))
        ms_approx = marginals(f_post_approx(x_test))
        @test mean.(ms_exact) ≈ mean.(ms_approx) rtol=1e-6
        @test std.(ms_exact) ≈ std.(ms_approx) rtol=1e-6
    end
    @testset "update_approx_posterior" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))

        # Reconstruction term and its gradient for Gaussian likelihood.
        r(m̃, σ̃²) = gaussian_reconstruction_term(y, Σ, m̃, σ̃²)

        @testset "optimal parameters don't change when updated" begin

            # Run a step of the update procedure.
            η1, η2 = natural_from_canonical(y, Σ)
            η1_new, η2_new, delta_norm = update_approx_posterior(f, x, η1, η2, r, 1.0)

            # Verify that the new parameters are equal to the old parameters.
            @test η1 ≈ η1_new
            @test η2 ≈ η2_new
            @test delta_norm < 1e-12
        end

        # Verify the optimal approximate posterior parameters are found after only a single
        # iteration of the algorithm in the Gaussian case.
        @testset "optimal in 1 step" begin

            # Perform a single update step with unit step size.
            N = length(y)
            η1, η2 = natural_from_canonical(y .+ randn(N), Σ + Diagonal(rand(N)))
            η1_opt, η2_opt, _ = update_approx_posterior(f, x, η1, η2, r, 1.0)

            # Ensure that we've reached the optimum (the exact posterior).
            y_opt, Σ_opt = canonical_from_natural(η1_opt, η2_opt)
            @test y ≈ y_opt
            @test Σ ≈ Σ_opt
        end
    end
    @testset "optimise_approx_posterior" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))

        # Reconstruction term and its gradient for Gaussian likelihood.
        r(m̃, σ̃²) = gaussian_reconstruction_term(y, Σ, m̃, σ̃²)

        @testset "converges quickly for step size $ρ" for ρ in [0.1, 0.5, 0.9, 1.0]
            η1, η2 = natural_from_canonical(y, Σ)
            η1_0 = η1 + randn(length(η1))
            η2_0 = η2 - Diagonal(rand(length(y)))
            η1_opt, η2_opt, _ = optimise_approx_posterior(f, x, η1_0, η2_0, r, ρ)
            @test η1_opt ≈ η1
            @test η2_opt ≈ η2
        end
    end
    @testset "elbo" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))

        # Set initial parameters to be sub-optimal.
        η1, η2 = natural_from_canonical(y, Σ)
        η1_0 = η1 + randn(length(η1))
        η2_0 = η2 - Diagonal(rand(length(y)))

        # Specify the reconstruction term.
        r(m̃, σ̃²) = gaussian_reconstruction_term(y, Σ, m̃, σ̃²)

        # Verify that the ELBO is less than the log marginal likelihood for a Gaussian.
        @test elbo(f, x, η1_0, η2_0, r) < logpdf(f(x, Σ), y)

        # Verify that the ELBO equals the log marginal likelihood at the optimum.
        @test elbo(f, x, η1, η2, r) ≈ logpdf(f(x, Σ), y)

        # Verify that the gradient w.r.t. everything can be computed using Zygote.
        Zygote.gradient(elbo, f, x, η1, η2, r)
    end
    @testset "batch_quadrature" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))
        function make_integrand(y, σ²)
            return f -> logpdf(Normal(f, sqrt(σ²)), y)
        end
        integrands = map(make_integrand, y, diag(Σ))

        # Specify the reconstruction term.
        r(m̃, σ̃²) = sum(batch_quadrature(integrands, m̃, sqrt.(diag(σ̃²)), 25))

        # Ensure that it's close to the exact reconstruction term in the Gaussian case.
        mq = randn(length(y))
        Σq = Diagonal(rand(length(y)) .+ 1)
        @test r(mq, Σq) ≈ gaussian_reconstruction_term(y, Σ, mq, Σq)

        # Check that the gradients agree.
        Δmq_exact, ΔΣq_exact = Zygote.gradient(
            (m̃, σ̃²) -> gaussian_reconstruction_term(y, Σ, m̃, σ̃²), mq, Σq,
        )
        Δmq_quad, ΔΣq_quad = Zygote.gradient(r, mq, Σq)
        @test Δmq_exact ≈ Δmq_quad
        @test ΔΣq_exact ≈ ΔΣq_quad
    end
    @testset "Gaussian likelihood via Quadature" begin
        f, x, Σ, y = generate_synthetic_problem(MersenneTwister(123456))
        function make_integrand(y, σ²)
            return f -> logpdf(Normal(f, sqrt(σ²)), y)
        end
        integrands = map(make_integrand, y, diag(Σ))

        # Set initial parameters to be sub-optimal.
        η1, η2 = natural_from_canonical(y, Σ)
        η1_0 = η1 + randn(length(η1))
        η2_0 = η2 - Diagonal(rand(length(y)))

        # Specify the reconstruction term.
        r(m̃, σ̃²) = sum(batch_quadrature(integrands, m̃, sqrt.(diag(σ̃²)), 25))

        # Ensure that it's close to the exact reconstruction term in the Gaussian case.
        mq = randn(length(y))
        Σq = Diagonal(rand(length(y)) .+ 1)
        @test r(mq, Σq) ≈ gaussian_reconstruction_term(y, Σ, mq, Σq)

        @testset "converges quickly for step size $ρ" for ρ in [0.1, 0.5, 0.9, 1.0]
            η1, η2 = natural_from_canonical(y, Σ)
            η1_0 = η1 + randn(length(η1))
            η2_0 = η2 - Diagonal(rand(length(y)))
            η1_opt, η2_opt, _, _ = optimise_approx_posterior(f, x, η1_0, η2_0, r, ρ)
            @test η1_opt ≈ η1
            @test η2_opt ≈ η2
        end
        @testset "elbo" begin
            # Verify that the ELBO is less than the log marginal likelihood for a Gaussian.
            @test elbo(f, x, η1_0, η2_0, r) < logpdf(f(x, Σ), y)

            # Verify that the ELBO equals the log marginal likelihood at the optimum.
            @test elbo(f, x, η1, η2, r) ≈ logpdf(f(x, Σ), y)

            # Verify that the gradient w.r.t. everything can be computed using Zygote.
            Zygote.gradient(elbo, f, x, η1, η2, r)
        end
    end
    @testset "Bernoulli likelihood via Quadrature" begin
        # Really all that we can do here is ensure that it converges and that the ELBO
        # is _something_.
        f, x, _, _ = generate_synthetic_problem(MersenneTwister(123456))
        y = rand.(Bernoulli.(logistic.(rand(f(x, 1e-4)))))
        function make_integrand(y)
            return f -> logpdf(Bernoulli(logistic(f)), y)
        end
        integrands = map(make_integrand, y)

        # Specify the reconstruction term.
        r(m̃, σ̃²) = sum(batch_quadrature(integrands, m̃, sqrt.(diag(σ̃²)), 25))

        # Set initial parameters to be... something.
        η1_0 = randn(length(y))
        η2_0 = Diagonal(- rand(length(y)))

        @testset "solutions converge to the same thing" begin
            solutions = map([0.1, 0.5, 0.9, 0.95, 0.99]) do ρ
                η1_opt, η2_opt, iteration, delta_norm = optimise_approx_posterior(
                    f, x, η1_0, η2_0, r, ρ,
                )
                return η1_opt, η2_opt
            end
            ref_solutions = first(solutions)
            @test all(isapprox.(first.(solutions), Ref(ref_solutions[1]); rtol=1e-6))
            @test all(isapprox.(last.(solutions), Ref(ref_solutions[2]); rtol=1e-6))
        end
        @testset "elbo" begin

            # Check that we can compute the ELBO. There's nothing obvious to compare it to.
            elbo(f, x, η1_0, η2_0, r)

            # Verify that the gradient w.r.t. everything can be computed using Zygote.
            Zygote.gradient(elbo, f, x, η1_0, η2_0, r)
        end
    end
end


text:
"""
A `Quantity` represents a time series of some quantity
measured in a tidal channel. You might measure water
levels (`Stage`) or flow (`Discharge`) or whatever
water quality parameters you might think up (not implemented).

Each sample of a quantity is assumed to be drawn at a
real-world time, so a `Quantity` is a vector of values of the
quantity in some space `T` attached to a `DateTime`.
"""
abstract type Quantity{T} <: AbstractVector{Tuple{DateTime,T}} end

"""Return the times at which the quantity is measured"""
times(q::Quantity) = q.ts

"""Return the values of the quantity itself"""
quantity(q::Quantity) = q.q

"""
Return both of the parts of the quantity
"""
unzip(q::Quantity) = (times(q),quantity(q))

Base.IndexStyle{T<:Quantity}(::Type{T}) = IndexLinear()

Base.size(q::Quantity) = size(times(q))
Base.getindex(q::Quantity,i::Int) = (times(q)[i],quantity(q)[i])

function Base.setindex!{T}(q::Quantity, v::Tuple{DateTime,T}, i::Int)
    q.ts[i] = v[1]
    q.q[i] = v[2]
    v
end

Base.similar{T}(q::Quantity,::Type{T},dims::Dims) = typeof(q)(similar(times(q),dims),similar(quantity(q),dims))

"""
    @quantity name T

This macro simplifies the process of defining
a quantity that conforms to the abstract
specification.

# Examples

```julia
@quantity Stage Float64
```

generates a new type Stage <: Quantity{Float64}
which is a Real-valued time series
"""
macro quantity(name,T)
    quote
        Base.@__doc__ struct $name{T<:$T} <: Quantity{T}
            ts::Vector{DateTime}
            q::Vector{T}
        end
    end
end

"""
    @quantity_op Quantity = Quantity (:op) Quantity

This macro defines a Quantity as the output of an
operation on two other Quantities. The operation is
applied by a map over the two input quantities, so
it is possible to define any operation which makes
sense pointwise in a time series.

# Examples

```julia
@quantity_op Discharge = CrossSectionalArea * AlongChannelVelocity
```

"""
macro quantity_op(ex)
    ex.head == :(=) || error("Expression is not an equality")
    output = ex.args[1]
    expr = ex.args[2]
    expr.head == :call || error("Expression does not call a function")
    op = expr.args[1]
    length(expr.args[2:end]) == 2 || error("Expressions of length > 2 quantities not yet supported")
    args = (expr.args[2:end]...)
    quote
        function Base.$op{A<:$(esc(args[1])),B<:$(esc(args[2]))}(a::A,b::B)
            times(a) == times(b) || error("Time series are not sampled at the same times. Consider interpolating")
            $(esc(output))(times(a),map($op, quantity(a), quantity(b)))
        end
    end
end

"""
An ``R``-valued time series for water levels
"""
@quantity Stage Real

"""
An ``R``-valued time series for the flooded
cross-sectional area of a channel
"""
@quantity CrossSectionalArea Real

"""
An ``R^3`` valued time series for velocity
"""
@quantity Velocity Tuple{Real,Real,Real}

"""
An ``R^3`` valued time series for the velocity
along a channel
"""
@quantity AlongChannelVelocity Real

"""
An ``R``-valued time series for discharge
"""
@quantity Discharge Real

@quantity_op Discharge = CrossSectionalArea * AlongChannelVelocity

"""
An ``R``-valued time series for temperature
"""
@quantity Temperature Real

"""
An ``R``-valued time series for total suspended solids
"""
@quantity TSS Real

"""
A ``Bool`-valued time series used to indicate valid (true) and invalid (false) data point
"""
@quantity Mask Bool


text:
using PyPlot
using MRIReco

#### trajectories ####

trajectories = ["Spiral", "Radial", "Cartesian", "EPI", "SpiralVarDens"]

figure(1)
clf()

tr = trajectory("Spiral", 1, 600, windings=10)
subplot(2,2,1)
nodes = kspaceNodes(tr)
plot(nodes[1,:], nodes[2,:],"b.",lw=2)
title("Spiral")

tr = trajectory("Cartesian", 13, 50, EPI_factor=1)
subplot(2,2,2)
nodes = kspaceNodes(tr)
plot(nodes[1,:], nodes[2,:],"b.",lw=2)
title("Cartesian")

tr = trajectory("Radial", 13, 50)
subplot(2,2,3)
nodes = kspaceNodes(tr)
plot(nodes[1,:], nodes[2,:],"b.",lw=2)
title("Radial")

tr = trajectory("SpiralVarDens", 3, 200)
subplot(2,2,4)
nodes = kspaceNodes(tr)
plot(nodes[1,:], nodes[2,:],"b.",lw=2)
title("SpiralVarDens")

subplots_adjust(top=0.95,
       bottom=0.05,
       left=0.09,
       right=0.95,
       hspace=0.3,
       wspace=0.2)

filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/trajectories.png")
savefig(filename, dpi=200)


text:
function my_function()
    sleep(2.0)
    println("Hello world! Current time is $(time())s")
end

my_function()


text:
using HDF5, JLD
using PyPlot
using FITSIO

###############################
########## Rec Image ##########
file = string("/Users/deguignet/Documents/Julia/imagerec_400ite_woshar.jld")
x = load(file, "admmst.x");
x = x[1024-512:1024+512-1,1024-512:1024+512-1,:];

file2 = string("/Users/deguignet/Documents/Julia/imagerec_700ite.jld")
x2 = load(file2, "admmst.x");
x2 = x2[1024-512:1024+512-1,1024-512:1024+512-1,:];

###############################
####### Original object #######
skyfile = string("/Users/deguignet/Documents/Julia/I_HALO_CL_RS_SKY.FITS")
f = FITS(skyfile)
sky = read(f[1]);
close(f)

###############################
############ Dirty ############
dirtyfile = string("/Users/deguignet/Documents/Julia/example_sim_dirty.fits")
f = FITS(dirtyfile)
mydata = read(f[1]);
close(f)

###############################
############# PSF #############
psffile = string("/Users/deguignet/Documents/Julia/example_sim_psf.fits")
f = FITS(psffile)
psf = read(f[1]);
close(f)

###############################
nxy  = size(x)[1]
lastiter = load(file, "algost.lastiter")
nfreq = load(file, "algost.nfreq")
tol1 = load(file, "toolst.tol1")
tol2 = load(file, "toolst.tol2")
tol3 = load(file, "toolst.tol3")
tol4 = load(file, "toolst.tol4")
tol5 = load(file, "toolst.tol5")
nu = load(file, "psfst.nu")
nu0 = load(file, "psfst.nu0")

sumsky2 = Array(Float64,nfreq)
for z in 1:nfreq
    sumsky2[z] = sum(sky[:,:,z].^2)
end
###############################

###############################
############ RMSE #############
# rmse = Array(Float64,nfreq)
# for z in 1:nfreq
#     rmse[z] = sqrt(sum( (x[:,:,z] - sky[:,:,z]).^2)/sumsky2[z])
# end
#
#
# figure(1)
# clf()
# plot(nu,rmse)
##########################
#
# err = Array(Float64,nxy,nxy,nfreq)
# for z in 1:nfreq
#     err[:,:,z] = (sky[:,:,z] - x[:,:,z]).^2./(sky[:,:,z].^2)
# end

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

##########################
#
# figure(2)
# clf()
# for z in 1:nfreq
#     subplot(6,6,z)
#     a = nu[z]/1e9
#     imshow(x[:,:,z],cmap=ColorMap("ocean"))
#     axis("off")
#     title("v = $a GHz")
#
# end

#############################################################################
#############################################################################
######################################  sky / dirty / rec_obj / error
# figure(3)
# clf()
# title("multiplot")
# e = 1
# for z in [8 16 24 32]
#     subplot(4,4,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow((abs(sky[:,:,z]).^0.25),cmap=ColorMap("spectral"))
#     title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e += 1
# end
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.3, 0.025, 0.6])
# colorbar(cax=cax)
# e = 5
# for z in [8 16 24 32]
#     subplot(4,4,e)
#     axis("off")
#     imshow((mydata[:,:,z]),cmap=ColorMap("spectral"))
#     e += 1
# end
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.3, 0.025, 0.6])
# colorbar(cax=cax)
#
# e = 9
# for z in [8 16 24 32]
#     subplot(4,4,e)
#     axis("off")
#     imshow((abs(x[:,:,z])/maximum(x[:,:,z])*maximum(sky[:,:,z])).^0.5,cmap=ColorMap("spectral"))
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.3, 0.025, 0.6])
# colorbar(cax=cax)
#
#
# e = 13
# for z in [8 16 24 32]
#     subplot(4,4,e)
#     axis("off")
#     imshow((abs(sky[:,:,z]-x[:,:,z])/(sum(sky[:,:,z]))).^0.25,cmap=ColorMap("spectral"))
#     e += 1
# end
#
# subplots_adjust(bottom=0.1, right=0.8, top=0.9)
# cax = axes([0.85, 0.1, 0.025, 0.17])
# colorbar(cax=cax)
#############################################################################
#############################################################################

# ################# PSF
# figure(4)
# clf()
# e=1
# for z in [8 16 24 32]
#
#     subplot(2,2,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow((abs(psf[:,:,z]).^0.25),cmap=ColorMap("spectral"))
#         title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e+=1
# end
#
# ################# Object
# figure(5)
# clf()
# e=1
# for z in [8 16 24 32]
#
#     subplot(2,2,e)
#     axis("off")
#     a = nu[z]/1e9
#     imshow((abs(sky[:,:,z]).^0.25),cmap=ColorMap("spectral"))
#         title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
#     e+=1
# end

##################################
##################################
figure(6)
clf()
subplot(2,3,1)
    axis("off")
    colorbar(imshow((abs(x[:,:,1]).^0.25),cmap=ColorMap("spectral")))
    title("400 iterations w/o shaar")
    ylabel("puissance 1/4")
subplot(2,3,2)
    axis("off")
    colorbar(imshow((abs(x2[:,:,1]).^0.25),cmap=ColorMap("spectral")))
    title("700 iterations")
subplot(2,3,3)
    axis("off")
    colorbar(imshow((abs(sky[:,:,1]).^0.25),cmap=ColorMap("spectral")))
    title("original object")

    subplot(2,3,4)
        axis("off")
        colorbar(imshow((abs(x[:,:,1])),cmap=ColorMap("spectral")))
        title("400 iterations w/o shaar")
        ylabel("puissance 1")

    subplot(2,3,5)
        axis("off")
        colorbar(imshow((abs(x2[:,:,1])),cmap=ColorMap("spectral")))
        title("700 iterations")
    subplot(2,3,6)
        axis("off")
        colorbar(imshow((abs(sky[:,:,1])),cmap=ColorMap("spectral")))
        title("original object")

# subplot(1,3,1)
#     axis("off")
#     imshow(((x[:,:,1]).^1),cmap=ColorMap("spectral"))
# subplot(1,3,2)
#     axis("off")
#     imshow(((x2[:,:,1]).^1),cmap=ColorMap("spectral"))
# subplot(1,3,3)
#     axis("off")
#     imshow(((sky[:,:,1]).^1),cmap=ColorMap("spectral"))

##################################
##################################
# figure(7)
# clf()
# xt = linspace(1,lastiter,lastiter)
# plot(xt,log(tol2))
# plot(xt,log(tol3))
# plot(xt,log(tol4))
# plot(xt,log(tol5))

##################################
##################################
figure(6)
clf()
subplot(1,2,1)
    axis("off")
    (imshow((abs(x[:,:,1]).^0.25),cmap=ColorMap("spectral"),vmax=0.1))

subplot(1,2,2)
    axis("off")
    (imshow((abs(sky[:,:,1]).^0.25),cmap=ColorMap("spectral"),vmax=0.12))

    # subplots_adjust(bottom=0.1, right=0.8, top=0.9)
    subplots_adjust(bottom=0.1, right=0.8, top=1.1)
    cax = axes([0.85, 0.3, 0.025, 0.6])
    colorbar(cax=cax)


text:
include("test_immutable.jl")
include("test_dispatch.jl")
include("test_ops.jl")


text:
# high-level wrappers

# create
gco_create(width, height, labelNum) = GCoptimizationGridGraph(width, height, labelNum)
gco_create(siteNum, labelNum) = GCoptimizationGeneralGraph(siteNum, labelNum)

# set neighbors
gco_setneighbors(handle, site1, site2) = setNeighbors(handle, site1-1, site2-1)
gco_setneighbors(handle, site1, site2, weight) = setNeighbors(handle, site1-1, site2-1, weight)
function gco_setneighbors(handle, weight::Matrix)
    r, c = size(weight)
    for j in 1:c, i in 1:j-1
        gco_setneighbors(handle, i, j, weight[i,j])
    end
end

# datacost
gco_setdatacost(handle, datacost::Vector) = setDataCost(handle, datacost)
gco_setdatacost(handle, datacost::Matrix) = gco_setdatacost(handle, [datacost...;])
gco_setdatacost(handle, site, label, energy) = setDataCost(handle, site-1, label-1, energy)

# smoothcost
gco_setsmoothcost(handle, smoothcost::Vector) = setSmoothCost(handle, smoothcost)
gco_setsmoothcost(handle, smoothcost::Matrix) = gco_setsmoothcost(handle, [smoothcost...;])
gco_setsmoothcost(handle, label1, label2, energy) = setSmoothCost(handle, label1-1, label2-1, energy)

# labelcost
gco_setlabelcost(handle, labelcost::Integer) = setLabelCost(handle, labelcost)
gco_setlabelcost(handle, labelcost::Vector) = setLabelCost(handle, labelcost)
gco_setlabelcost(handle, cost, labelIndices) = setLabelSubsetCost(handle, labelIndices-1, length(labelIndices), cost)

# expansion & swap
gco_expansion(handle, maxIterationNum=-1) = expansion(handle, maxIterationNum)
gco_αexpansion(handle, αlabel) = alpha_expansion(handle, αlabel-1)
gco_swap(handle, maxIterationNum=-1) = swap(handle, maxIterationNum)
gco_αβswap(handle, αlabel, βlabel) = alpha_beta_swap(handle, αlabel-1, βlabel-1)

# compute_energy
gco_energy(handle) = return compute_energy(handle), giveDataEnergy(handle), giveSmoothEnergy(handle), giveLabelEnergy(handle)

# get labeling
gco_getlabeling(handle) = [whatLabel(handle, i-1)+1 for i in 1:numSites(handle)]

function gco_getlabeling(handle, site)
    @assert 0 < site <= numSites(handle) "invalid site: $site"
    return whatLabel(handle, site-1)+1
end

function gco_getlabeling(handle, sites::Vector)
    result = Vector{Int}(length(sites))
    for i in eachindex(sites)
        @assert 0 < sites[i] <= numSites(handle) "invalid site: $(sites[i])"
        result[i] = whatLabel(handle, sites[i]-1)+1
    end
    return result
end

# set labeling
function gco_setlabeling(handle, labeling::Vector)
    @assert 0 < length(labeling) <= numSites(handle) "no label or too many labels... $(length(labeling))"
    for i in eachindex(labeling)
        @assert 0 < labeling[i] <= numLabels(handle) "invalid label: $(labeling[i])"
        setLabel(handle, i-1, labeling[i]-1)
    end
end
gco_setlabeling(handle, labeling::Matrix) = gco_setlabeling(handle, [labeling...;])

function gco_setlabeling(handle, sites::Vector, labeling::Vector)
    @assert 0 < length(labeling) <= numSites(handle) "no label or too many labels... $(length(labeling))"
    for i in eachindex(sites)
        @assert 0 < labeling[i] <= numLabels(handle) "invalid label: $(labeling[i])"
        setLabel(handle, sites[i]-1, labeling[i]-1)
    end
end
gco_setlabeling(handle, sites::Matrix, labeling::Matrix) = gco_setlabeling(handle, [sites...;], [labeling...;])

# set label order
gco_setlabelorder(handle, random::Bool) = setLabelOrder(handle, random)
gco_setlabelorder(handle, order, size) = setLabelOrder(handle, order, size)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UnionProcess(p₁, p₂)

Union (or superposition) of spatial point processes `p₁` and `p₂`.
"""
struct UnionProcess{P₁<:PointProcess,P₂<:PointProcess} <: PointProcess
  p₁::P₁
  p₂::P₂
end

"""
    p₁ ∪ p₂

Return the union of point processes `p₁` and `p₂`.
"""
Base.union(p₁::PointProcess, p₂::PointProcess) = UnionProcess(p₁, p₂)

ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)

default_sampling_algorithm(::UnionProcess, ::Geometry) = UnionSampling()

function rand_single(rng::Random.AbstractRNG,
                     p::UnionProcess, g::Geometry,
                     ::UnionSampling)
  pp₁ = rand(rng, p.p₁, g)
  pp₂ = rand(rng, p.p₂, g)

  PointSet([coordinates.(pp₁); coordinates.(pp₂)])
end