text:
@testset "find_common_clusters" begin
    ref_tree = ParseNewick("((A,(B,(C,(D,E)))),(F,(G,H)));")
    tree2 = ParseNewick("((G,(C,(A,(F,E)))),(B,(D,H)));")

    A, B, C, D, E, F, G, H = find_by_name(tree2, "A"), find_by_name(tree2, "B"),
                             find_by_name(tree2, "C"), find_by_name(tree2, "D"),
                             find_by_name(tree2, "E"), find_by_name(tree2, "F"),
                             find_by_name(tree2, "G"), find_by_name(tree2, "H")
     expected_dict = Dict{Int64,Tuple{Bool,Union{Missing, Float64}}}([(D.mother.num, (false, missing)), (C.mother.num, (false, missing)), (B.mother.num, (false, missing)),
                           (A.mother.num, (false, missing)), (G.mother.num, (false, missing)), (F.mother.num, (false, missing)),
                           (G.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                           (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhyloTree.find_common_clusters(ref_tree, tree2) , expected_dict)

    tree3 = ParseNewick("((A,(C,(D,(B,E)))),(G,(F,H)));")

    A, B, C, D, E, F, G, H = find_by_name(tree3, "A"), find_by_name(tree3, "B"),
                             find_by_name(tree3, "C"), find_by_name(tree3, "D"),
                             find_by_name(tree3, "E"), find_by_name(tree3, "F"),
                             find_by_name(tree3, "G"), find_by_name(tree3, "H")
     expected_dict = Dict([(D.mother.num, (false, missing)), (C.mother.num, (true, 1.0)), (B.mother.num, (false, missing)),
                           (A.mother.num, (true, 1.0)), (G.mother.num, (true, 1.0)), (F.mother.num, (false, missing)),
                           (A.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                           (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhyloTree.find_common_clusters(ref_tree, tree3), expected_dict)

    tree4 = ParseNewick("((A,(B,(C,(D,E)))),(F,(G,H)));")

    A, B, C, D, E, F, G, H = find_by_name(tree4, "A"), find_by_name(tree4, "B"),
                             find_by_name(tree4, "C"), find_by_name(tree4, "D"),
                             find_by_name(tree4, "E"), find_by_name(tree4, "F"),
                             find_by_name(tree4, "G"), find_by_name(tree4, "H")
    expected_dict = Dict([(D.mother.num, (true, 1.0)), (C.mother.num, (true, 1.0)), (B.mother.num, (true, 1.0)),
                          (A.mother.num, (true, 1.0)), (G.mother.num, (true, 1.0)), (F.mother.num, (true, 1.0)),
                          (A.mother.mother.num, (true, 1.0)), (A.num, (true, 1.0)), (B.num, (true, 1.0)), (C.num, (true, 1.0)),
                          (D.num, (true, 1.0)), (E.num, (true, 1.0)), (F.num, (true, 1.0)), (G.num, (true, 1.0)), (H.num, (true, 1.0))])
    @test isequal(MCPhyloTree.find_common_clusters(ref_tree, tree4), expected_dict)

    tree5 = ParseNewick("((G,(X,(A,(F,E)))),(B,(D,H)));")

    @test_throws ArgumentError MCPhyloTree.find_common_clusters(ref_tree, tree5)

    tree6 = ParseNewick("(X,(G,(C,(A,(F,E)))),(B,(D,H)));")

    @test_throws ArgumentError MCPhyloTree.find_common_clusters(ref_tree, tree6)
end

@testset "one_way_compatible" begin
    tree = ParseNewick("((A,C,E),(B,D));")
    tree2 = ParseNewick("((A,C),(B,D,E));")
    expected_tree = ParseNewick("(A,C,E,(B,D));")

    @test newick(MCPhyloTree.one_way_compatible(tree, tree2)) == newick(expected_tree)
end

@testset "get_leaf_ranks" begin
    ref_tree = ParseNewick("((A,(B,(C,(D,E)))),(F,(G,H)));")
    nodes = post_order(ref_tree)
    @test MCPhyloTree.get_leaf_ranks(nodes) == Dict([("A", 1), ("B", 2), ("C", 3),
                                                 ("D", 4), ("E", 5), ("F", 6),
                                                 ("G", 7), ("H", 8)])
end

@testset "order_tree!" begin
    tree = ParseNewick("(A,B,(C,(D,E)F)G)H;")

    A, B, C, D, E, F, G, H = find_by_name(tree, "A"), find_by_name(tree, "B"),
                             find_by_name(tree, "C"), find_by_name(tree, "D"),
                             find_by_name(tree, "E"), find_by_name(tree, "F"),
                             find_by_name(tree, "G"), find_by_name(tree, "H")
    cluster_start_indeces = Dict([(A, 3), (B, 7), (C, 2), (D, 8),
                                  (E, 5), (F, 1), (G, 4), (H, 6)])
    ordered_tree = ParseNewick("(A,((E,D)F,C)G,B)H;")
    @test MCPhyloTree.order_tree!(tree, cluster_start_indeces) == [A, E, D, C, B]
    @test newick(tree) == newick(ordered_tree)
end

@testset "max/min_leaf_rank" begin
    tree = ParseNewick("(A,B,(C,(D,E)F)G)H;")
    F, G, H, A = find_by_name(tree, "F"), find_by_name(tree, "G"),
                 find_by_name(tree, "H"), find_by_name(tree, "A")
    leaf_ranks = Dict([("A", 1), ("B", 5), ("C", 2), ("D", 4), ("E", 3)])

    @testset "min_leaf_rank" begin
        @test MCPhyloTree.min_leaf_rank(leaf_ranks, F) == 3
        @test MCPhyloTree.min_leaf_rank(leaf_ranks, G) == 2
        @test MCPhyloTree.min_leaf_rank(leaf_ranks, H) == 1
        @test MCPhyloTree.min_leaf_rank(leaf_ranks, A) == 1
    end

    @testset "max_leaf_rank" begin
        @test MCPhyloTree.max_leaf_rank(leaf_ranks, F) == 4
        @test MCPhyloTree.max_leaf_rank(leaf_ranks, G) == 4
        @test MCPhyloTree.max_leaf_rank(leaf_ranks, H) == 5
        @test MCPhyloTree.max_leaf_rank(leaf_ranks, A) == 1
    end
end

@testset "x_left_right" begin
    tree = ParseNewick("(A,B,(C,(D,E)F)G)H;")

    A, B, C, D, E, F, G, H = find_by_name(tree, "A"), find_by_name(tree, "B"),
                             find_by_name(tree, "C"), find_by_name(tree, "D"),
                             find_by_name(tree, "E"), find_by_name(tree, "F"),
                             find_by_name(tree, "G"), find_by_name(tree, "H")

    @testset "x_left" begin
        @test MCPhyloTree.x_left(A) == (H, [A,H])
        @test MCPhyloTree.x_left(B) == (B, [B,H])
        @test MCPhyloTree.x_left(C) == (G, [C,G,H])
        @test MCPhyloTree.x_left(D) == (F, [D,F,G])
        @test MCPhyloTree.x_left(E) == (E, [E,F])
    end

    @testset "x_right" begin
        @test MCPhyloTree.x_right(A) == (A, [A,H])
        @test MCPhyloTree.x_right(B) == (B, [B,H])
        @test MCPhyloTree.x_right(C) == (C, [C,G])
        @test MCPhyloTree.x_right(D) == (D, [D,F])
        @test MCPhyloTree.x_right(E) == (H, [E,F,G,H])
    end
end

@testset "majority_consensus_tree" begin
    tree1 = ParseNewick("(((A,B),C),(D,E));")
    tree2 = ParseNewick("((A,C),(B,D,E));")
    tree3 = ParseNewick("(((B,C),A),D,E);")
    trees = [tree1, tree2, tree3]

    result = newick(ParseNewick("((A,B,C),D,E);"))
    @test newick(majority_consensus_tree(trees)) == result

    tree4 = ParseNewick("(((B,C),A),D,F);")
    push!(trees, tree4)
    @test_throws ArgumentError majority_consensus_tree(trees)

    tree1 = ParseNewick("(((A,B)AB,C)ABC,(D,E)DE)no_name;")
    tree2 = ParseNewick("((A,C)AC,(B,D,E)BDE)no_name;")
    tree3 = ParseNewick("(((B,C)BC,A)BCA,D,E)no_name;")
    trees = [tree1, tree2, tree3]

    result = newick(ParseNewick("((A,B,C),D,E);"))
    @test newick(majority_consensus_tree(trees)) == result
end

@testset "loose_consensus_tree" begin
    tree1 = ParseNewick("(((A,B),C),(D,E));")
    tree2 = ParseNewick("((A,C),(B,D,E));")
    tree3 = ParseNewick("(((B,C),A),D,E);")
    trees = [tree1, tree2, tree3]

    result = newick(ParseNewick("(A,B,C,(D,E));"))
    @test newick(loose_consensus_tree(trees)) == result

    tree4 = ParseNewick("(((B,C),A),D,F);")
    push!(trees, tree4)
    @test_throws ArgumentError loose_consensus_tree(trees)
    tree1 = ParseNewick("(((A,B)AB,C)ABC,(D,E)DE)no_name;")
    tree2 = ParseNewick("((A,C)AC,(B,D,E)BDE)no_name;")
    tree3 = ParseNewick("(((B,C)BC,A)BCA,D,E)no_name;")
    trees = [tree1, tree2, tree3]

    result = newick(ParseNewick("(A,B,C,(D,E)no_name);"))
    @test newick(loose_consensus_tree(trees)) == result
end

@testset "greedy_consensus_tree" begin
    tree1 = ParseNewick("(((A,B),C),(D,E));")
    tree2 = ParseNewick("((A,C),(B,D,E));")
    tree3 = ParseNewick("(((B,C),A),D,E);")
    trees = [tree1, tree2, tree3]

    result = newick(ParseNewick("(((A,B),C),(D,E));"))
    # result = newick(ParseNewick("(((A,C),B),(D,E));"))
    # result = newick(ParseNewick("(((B,C),A),(D,E));"))
    @test newick(greedy_consensus_tree(trees)) == result
end

@testset "count_cluster_occurences" begin
    bit_vec = BitArray.([[1,0,0,0,0], [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0],
              [0,0,0,0,1], [1,1,0,0,0], [1,1,1,0,0], [0,0,0,1,1], [1,0,0,0,0],
              [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0], [0,0,0,0,1], [1,0,1,0,0],
              [0,1,0,1,1], [0,1,1,0,0], [1,1,1,0,0], [0,0,0,1,1], [1,0,0,0,0],
              [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0], [0,0,0,0,1]])
    sort!(bit_vec, alg=QuickSort)
    queue = MCPhyloTree.PriorityQueue{BitVector, Int64}()
    queue[BitArray([0,0,0,0,1])] = -3
    queue[BitArray([0,0,0,1,0])] = -3
    queue[BitArray([0,1,0,0,0])] = -3
    queue[BitArray([0,0,1,0,0])] = -3
    queue[BitArray([1,0,0,0,0])] = -3
    queue[BitArray([0,0,0,1,1])] = -2
    queue[BitArray([1,1,0,0,0])] = -1
    queue[BitArray([1,0,1,0,0])] = -1
    queue[BitArray([0,1,0,1,1])] = -1
    queue[BitArray([0,1,1,0,0])] = -1
    queue[BitArray([1,1,1,0,0])] = -2

    @test MCPhyloTree.count_cluster_occurences(bit_vec) == queue
end


text:
import JSON
# mock a listener:
Core.eval(Atom, Meta.parse("sock = IOBuffer()"))
readmsg() = JSON.parse(String(take!(Atom.sock)))

# callback count
cb = 0

# check the different message handlers:
# pingpong
Atom.handlemsg(Dict("type" => "ping",
                    "callback" => (cb += 1)))
@test readmsg() == ["cb", cb, "pong"]

# echo
Atom.handlemsg(Dict("type" => "echo",
                    "callback" => (cb += 1)),
                    "echome!")
@test readmsg() == ["cb", cb, "echome!"]

# cd
old_path = pwd()
Atom.handlemsg(Dict("type" => "cd",
                    "callback" => (cb += 1)),
                    joinpath(old_path, ".."))
@test readmsg() == ["cb", cb, nothing]
@test pwd() == realpath(joinpath(old_path, ".."))
cd(old_path)

# evalsimple
Atom.handlemsg(Dict("type" => "evalsimple",
                    "callback" => (cb += 1)),
                    "1+1")
@test readmsg() == ["cb", cb, 2]

Atom.handlemsg(Dict("type" => "evalsimple",
                    "callback" => (cb += 1)),
                    "sin(pi)")
@test readmsg() == ["cb", cb, sin(pi)]

ext: jl


text:
using MacroTools

# Tail call operations

function lastcalls(f, ex)
  isexpr(ex, :block) && return :(begin $(lastcalls(f, ex.args)...) end)
  @match ex begin
    let __ end     => :(let $(lastcalls(f, ex.args)...) end)
    (c_ ? y_ : n_) => :($c ? $(lastcalls(f, y)) : $(lastcalls(f, n)))
    (a_ && b_)     => :($a && $(lastcalls(f, b)))
    (a_ || b_)     => :($a || $(lastcalls(f, b)))
    _              => f(ex)
  end
end

lastcalls(f, ex::Array) =
  isempty(ex) ? ex :
    [ex[1:end-1]..., lastcalls(f, ex[end])]

function retcalls(f, ex)
  MacroTools.postwalk(ex) do ex
    @capture(ex, return x_) ? :(return $(lastcalls(f, x))) : ex
  end
end

tailcalls(f, ex) = @>> ex lastcalls(f) retcalls(f)

# Tail recursion

export @rec, @bounce

"Generate an expression like `(a, b) = (c, d)`."
tupleassign(xs, ys) = Expr(:(=), Expr(:tuple, xs...), Expr(:tuple, ys...))

tco(ex, f, dummy, start) =
  @capture(ex, $f(args__)) ?
    :($(tupleassign(dummy, args)); @goto $start) :
    ex

"""
Enables efficient recursive functions, e.g.

    @rec reduce(f::Function, v, xs::List) =
      isempty(xs) ? v : reduce(f, f(v, first(xs)), tail(xs))

Without `@rec` this function would overflow the stack for
lists of 80,000 or more elements.

Caveats:

  • No support for trampolining, i.e. only calls to the
    given function are optimised away.
  • Ignores multiple dispatch – it is assumed that the function's
    name always refers to the given definition.
  • Don't rebind the function's name in a let (see above).
  • Don't use this with varargs functions.

Use the more flexible, but slower, `@bounce` to avoid these issues.
"""
macro rec(def)
  def = shortdef(macroexpand(def))
  @capture(def, f_(args__) = body_) || error("@rec: $def is not a function definition.")
  f = namify(f)
  dummy = @>> args map(namify) map(string) map(gensym)

  # Set up of variables
  @gensym start
  body = quote
           $(tupleassign(dummy, args))
           @label $start
           $(tupleassign(args, dummy))
           $body
         end

  def.args[2] = tailcalls(ex -> tco(ex, f, dummy, start), body)
  return esc(def)
end

# Trampolining

trampname(f) = Symbol(string("#__", f, "_tramp__"))

mutable struct Bounce
  f::Function
end

function trampoline(f, args...)
  val = f(args...)
  while isa(val, Bounce)
    val = val.f()
  end
  return val
end

function bounce(ex)
  @capture(ex, f_(args__)) || return ex
  f_tramp = trampname(f)
  :(Lazy.Bounce(() -> $f_tramp($(args...))))
end

function trampdef(f)
  f_tramp = trampname(f)
  isdefined(f_tramp) && return
  :($f_tramp(args...) = $f(args...))
end

"""
Tail recursion that doesn't blow the stack.

    @bounce even(n) = n == 0 ? true : odd(n-1)
    @bounce odd(n) = n == 0 ? false : even(n-1)

    even(1_000_000) # Blows up without `@bounce`.
    #> true

For simple cases you probably want the much faster `@rec`.
"""
macro bounce(def)
  def = macroexpand(def)
  @assert isdef(def)
  @assert isexpr(def.args[1].args[1], Symbol) # TODO: handle f{T}() = ...
  f = namify(def)
  f_tramp = trampname(f)
  args = def.args[1].args[2:end]
  def.args[1].args[1] = f_tramp

  calls = Symbol[]
  def.args[2] = tailcalls(ex -> (isexpr(ex, :call) && push!(calls, ex.args[1]);
                                 bounce(ex)),
                          def.args[2])
  quote
    $([trampdef(call) for call in calls]...)
    $def
    $f($(args...)) = Lazy.trampoline($f_tramp, $(args...))
  end |> esc
end


text:
d = TruncatedNormal(Normal(0, 1), -0.1, +0.1)

@assert pdf(d, 0.0) > pdf(Normal(0, 1), 0.0)
@assert pdf(d, -1.0) == 0.0
@assert pdf(d, +1.0) == 0.0

@assert logpdf(d, 0.0) > logpdf(Normal(0, 1), 0.0)
@assert isinf(logpdf(d, -1.0))
@assert isinf(logpdf(d, +1.0))

@assert cdf(d, -1.0) == 0.0
@assert cdf(d, -0.09) < cdf(Normal(0, 1), -0.09)
@assert cdf(d, 0.0) == 0.5
@assert cdf(d, +0.09) > cdf(Normal(0, 1), +0.09)
@assert cdf(d, +1.0) == 1.0

@assert quantile(d, 0.01) > -0.1
@assert abs(quantile(d, 0.5) - 0.0) < 1e-8
@assert quantile(d, 0.99) < +0.1

@assert abs(cdf(d, quantile(d, 0.01)) - 0.01) < 1e-8
@assert abs(cdf(d, quantile(d, 0.50)) - 0.50) < 1e-8
@assert abs(cdf(d, quantile(d, 0.99)) - 0.99) < 1e-8


text:
export SuperOpSpace

abstract type AbstractSuperOpBasis <: DoubledBasis end
@inline nsites(h::AbstractSuperOpBasis) = 2*nsites(physical(h))
@inline spacedimension(h::AbstractSuperOpBasis) = spacedimension(physical(h))^2

state(arrT::AbstractArray, T::Type{<:Number}, h::AbstractSuperOpBasis, dims::Dims) =
    (state(arrT, T, physical(h), dims),
     state(arrT, T, physical(h), dims))

mutable struct SuperOpSpace{H<:AbstractHilbert} <: AbstractSuperOpBasis
    basis::H
end

Base.show(io::IO, ::MIME"text/plain", h::SuperOpSpace) =
    print(io, "Hilbert Space for Superoperators with $(physical(h)))")

Base.show(io::IO, h::SuperOpSpace) =
    print(io, "SuperOpSpace($(physical(h))))")

##

function set!(σ::ADoubleState, h::SuperOpSpace, i::Integer)
    @boundscheck checkbounds_hilbert(h, i)

    hp = physical(h)
    np = nsites_physical(h)
    i -= 1

    i_r = div(i, spacedimension(hp))
    i_c = i - i_r*spacedimension(hp)

    @inbounds set!(col(σ), hp, i_r + 1) #i
    @inbounds set!(row(σ), hp, i_c + 1) #j
    return σ
end

Base.@propagate_inbounds function set!(σ::ADoubleState, h::SuperOpSpace, i_r::Integer, i_c::Integer)
    set!(row(σ), physical(h), i_r)
    set!(col(σ), physical(h), i_c)
    return σ
end

@inline function Random.rand!(rng::AbstractRNG, σ::AbstractDoubled, h::SuperOpSpace)
    rand!(rng, row(σ), physical(h))
    rand!(rng, col(σ), physical(h))
    return σ
end

function toint(σ::ADoubleState, h::SuperOpSpace)
    hp = physical(h)

    left = toint(col(σ), hp) - 1
    right = toint(row(σ), hp) - 1
    return left * spacedimension(hp) + right + 1
end


text:
abstract type AbstractGAMSCmdAttribute <: MOI.AbstractOptimizerAttribute end

struct ModelType <: MOI.AbstractOptimizerAttribute end
struct SysDir <: MOI.AbstractOptimizerAttribute end
struct WorkDir <: MOI.AbstractOptimizerAttribute end

struct HoldFixed <: AbstractGAMSCmdAttribute end
struct IterLim <: AbstractGAMSCmdAttribute end
struct License <: AbstractGAMSCmdAttribute end
struct LogOption <: AbstractGAMSCmdAttribute end
struct NodLim <: AbstractGAMSCmdAttribute end
struct OptCA <: AbstractGAMSCmdAttribute end
struct OptCR <: AbstractGAMSCmdAttribute end
struct ResLim <: AbstractGAMSCmdAttribute end
struct Solver <: AbstractGAMSCmdAttribute end
struct Threads <: AbstractGAMSCmdAttribute end
struct Trace <: AbstractGAMSCmdAttribute end
struct TraceOpt <: AbstractGAMSCmdAttribute end

struct CNS <: AbstractGAMSCmdAttribute end
struct DNLP <: AbstractGAMSCmdAttribute end
struct LP <: AbstractGAMSCmdAttribute end
struct MCP <: AbstractGAMSCmdAttribute end
struct MINLP <: AbstractGAMSCmdAttribute end
struct MIP <: AbstractGAMSCmdAttribute end
struct MIQCP <: AbstractGAMSCmdAttribute end
struct MPEC <: AbstractGAMSCmdAttribute end
struct NLP <: AbstractGAMSCmdAttribute end
struct QCP <: AbstractGAMSCmdAttribute end
struct RMINLP <: AbstractGAMSCmdAttribute end
struct RMIP <: AbstractGAMSCmdAttribute end
struct RMIQCP <: AbstractGAMSCmdAttribute end

function MOI.get(
   model::Optimizer,
   opt::AbstractGAMSCmdAttribute
)
   name = lowercase(replace(string(typeof(opt)), r"(GAMS.)" => ""))
   if haskey(model.gams_options, name)
      return model.gams_options[name]
   end
   error("GAMS option '$(name)' is not set.")
end

function MOI.set(
   model::Optimizer,
   opt::AbstractGAMSCmdAttribute,
   value
)
   name = lowercase(replace(string(typeof(opt)), r"(GAMS.)" => ""))
   model.gams_options[name] = value
   return
end

function MOI.get(
   model::Optimizer,
   ::ModelType
)
   return label(model.user_model_type)
end

function MOI.set(
   model::Optimizer,
   ::ModelType,
   value::String
)
   value = uppercase(value)
   model.user_model_type = model_type_from_label(value)
   if model.user_model_type == MODEL_TYPE_UNDEFINED
      error("Unsupported model type '$value'.")
   end
   return
end

function MOI.get(
   model::Optimizer,
   ::SysDir
)
   if isnothing(model.gamswork)
      if isnothing(model.sysdir)
         error("GAMS system directory has not been set.")
      end
      return model.sysdir
   else
      return model.gamswork.system_dir
   end
end

function MOI.set(
   model::Optimizer,
   ::SysDir,
   value::String
)
   if isnothing(model.gamswork)
      check_system_dir(value)
      model.sysdir = value
   else
      set_system_dir(model.gamswork, value)
   end
   return
end

function MOI.get(
   model::Optimizer,
   ::WorkDir
)
   if isnothing(model.gamswork)
      if isnothing(model.workdir)
         error("GAMS working directory has not been set.")
      end
      return model.workdir
   else
      return model.gamswork.working_dir
   end
end

function MOI.set(
   model::Optimizer,
   ::WorkDir,
   value::String
)
   if isnothing(model.gamswork)
      model.workdir = value
   else
      set_working_dir(model.gamswork, value)
   end
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.Silent
)
   return get(model.gams_options, "logoption", nothing) == 0
end

function MOI.set(
   model::Optimizer,
   ::MOI.Silent,
   silent::Bool
)
   if silent
      MOI.set(model, MOI.RawParameter("logoption"), 0)
   else
      MOI.set(model, MOI.RawParameter("logoption"), 1)
   end
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.TimeLimitSec
)
   return get(model.gams_options, "reslim", nothing)
end

function MOI.set(
   model::Optimizer,
   ::MOI.TimeLimitSec,
   value::Real
)
   MOI.set(model, MOI.RawParameter("reslim"), Float64(value))
   return
end

function MOI.set(
   model::Optimizer,
   ::MOI.TimeLimitSec,
   ::Nothing
)
   delete!(model.gams_options, "reslim")
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfThreads
)
   return get(model.gams_options, "threads", nothing)
end

function MOI.set(
   model::Optimizer,
   ::MOI.NumberOfThreads,
   value::Int
)
   MOI.set(model, MOI.RawParameter("threads"), value)
   return
end

function MOI.get(
   model::Optimizer,
   option::MOI.RawParameter
)
   name = lowercase(option.name)

   if name == "modeltype"
      return MOI.get(model, ModelType())
   elseif name == "sysdir"
      return MOI.get(model, SysDir())
   elseif name == "workdir"
      return MOI.get(model, WorkDir())
   end

   if name in GAMS.MOI_SUPPORTED_CLOPTIONS
      if haskey(model.gams_options, name)
         return model.gams_options[name]
      end
      error("GAMS option '$(option.name)' is not set.")
   else
      if haskey(model.solver_options, option.name)
         return model.solver_options[option.name]
      end
      error("GAMS solver option '$(option.name)' is not set.")
   end
end

function MOI.set(
   model::Optimizer,
   option::MOI.RawParameter,
   value::String
)
   name = lowercase(option.name)

   if name == "modeltype"
      MOI.set(model, ModelType(), value)
      return
   elseif name == "sysdir"
      MOI.set(model, SysDir(), value)
      return
   elseif name == "workdir"
      MOI.set(model, WorkDir(), value)
      return
   end

   if name in GAMS.MOI_SUPPORTED_CLOPTIONS
      model.gams_options[name] = lowercase(value)
   else
      model.solver_options[option.name] = value
   end
   return
end

function MOI.set(
   model::Optimizer,
   option::MOI.RawParameter,
   value::Number
)
   name = lowercase(option.name)
   if name in GAMS.MOI_SUPPORTED_CLOPTIONS
      model.gams_options[name] = value
   else
      model.solver_options[option.name] = value
   end
   return
end


text:
using Test
using KMLTracks
using TimeZones

@testset "parse_kml_string" begin
    s = """<?xml version="1.0" encoding="UTF-8"?>
    <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
        <ExtendedData>
            <Data name="GPSModelName">
                <value>MyGPS</value>
            </Data>
        </ExtendedData>
        <Placemark>
            <gx:Track>
                <altitudeMode>absolute</altitudeMode>
                    <when>2019-11-07T08:44:12Z</when><gx:coord>0.309042 46.586349 111.4</gx:coord>
                    <when>2019-11-07T08:44:16Z</when><gx:coord>0.309042 46.586349 111.1</gx:coord>
                    <when>2019-11-07T08:44:20Z</when><gx:coord>0.309042 46.586349 111.3</gx:coord>
            </gx:Track>
        </Placemark>
    </kml>"""
    kmldoc = parse_kml_string(s)
    @test length(kmldoc.placemark.track.points) == 3
    pt = kmldoc.placemark.track.points[1]
    @test pt.time == ZonedDateTime(2019, 11, 07, 8, 44, 12, tz"UTC")
    @test pt.longitude == 0.309042
    @test pt.latitude == 46.586349
    @test pt.altitude == 111.4
end

@testset "read_kml_file" begin
    fname = "sample.kml"
    kmldoc = read_kml_file(fname)
    @test length(kmldoc.placemark.track.points) == 45
    pt = kmldoc.placemark.track.points[1]
    @test pt.time == ZonedDateTime(2019, 11, 07, 8, 44, 12, tz"UTC")
    @test pt.longitude == 0.309042
    @test pt.latitude == 46.586349
    @test pt.altitude == 111.4

    pt = kmldoc.placemark.track.points[45]
    @test pt.time == ZonedDateTime(2019, 11, 07, 8, 47, 8, 1, tz"UTC")
    @test pt.longitude == 0.295266
    @test pt.latitude == 46.575615
    @test pt.altitude == 267.6
end

@testset "parse_kml_string_coma_issue" begin
    s = """<?xml version="1.0" encoding="UTF-8"?>
    <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
        <ExtendedData>
            <Data name="GPSModelName">
                <value>MyGPS</value>
            </Data>
        </ExtendedData>
        <Placemark>
            <gx:Track>
                <altitudeMode>absolute</altitudeMode>
                    <when>2019-11-07T08:44:12Z</when><gx:coord>0,309042 46,586349 111,4</gx:coord>
                    <when>2019-11-07T08:44:16Z</when><gx:coord>0,309042 46,586349 111,1</gx:coord>
                    <when>2019-11-07T08:44:20Z</when><gx:coord>0,309042 46,586349 111,3</gx:coord>
            </gx:Track>
        </Placemark>
    </kml>"""
    kmldoc = parse_kml_string(s)
    @test length(kmldoc.placemark.track.points) == 3
    pt = kmldoc.placemark.track.points[1]
    @test pt.time == ZonedDateTime(2019, 11, 07, 8, 44, 12, tz"UTC")
    @test pt.longitude == 0.309042
    @test pt.latitude == 46.586349
    @test pt.altitude == 111.4
end


text:
#=
Created on Tuesday 2 July 2019
Last update: Saturday 27 July 2019

@author: Michiel Stock
michielfmstock@gmail.com

Benchmarking Kroncker.jl compated to native functions.
=#

using Kronecker, Plots, LinearAlgebra
using BenchmarkTools

sizes = [5, 10, 25, 50, 100, 250, 500, 1000, 5000]

tmax = 10

# inverse

times_kron = []
times_naive = []

for s in sizes
    A = randn(s, s) / s^2
    B = rand(s, s) / s^2
    v = rand(s^2)
    K = A ⊗ B
    t = @belapsed inv($K)
    push!(times_kron, t)
    if s^2 < 5000
        K = collect(K)
        t = @belapsed inv($K)
        push!(times_naive, t)
        compute_naive = t < tmax
    end
end

plot(sizes.^2, times_kron, label="matrix inverse", color=:blue, legend=:bottomright, lw=2)
plot!(sizes[1:length(times_naive)].^2, times_naive, color=:blue, ls=:dash, label="", lw=2 )

# determinant

times_kron = []
times_naive = []

@belapsed det(Kwarmup)

for s in sizes
    A = randn(s, s)
    B = rand(s, s)
    v = rand(s^2)
    K = A ⊗ B
    t = @belapsed det($K)
    push!(times_kron, t)
    if s^2 < 5000
        K = collect(K)
        t = @belapsed det($K)
        push!(times_naive, t)
        compute_naive = t < tmax
    end
end

plot!(sizes.^2, times_kron, label="determinant", color=:green, lw=2)
plot!(sizes[1:length(times_naive)].^2, times_naive, color=:green, ls=:dash, label="", lw=2)

# squaring

times_kron = []
times_naive = []

v = rand(100)

@belapsed Kwarmup * Kwarmup

for s in sizes
    A = randn(s, s)
    B = rand(s, s)
    v = rand(s^2)
    K = A ⊗ B
    t = @belapsed $K * $K
    push!(times_kron, t)
    if s^2 < 5000
        K = collect(K)
        t = @belapsed $K * $K
        push!(times_naive, t)
        compute_naive = t < tmax
    end
end

plot!(sizes.^2, times_kron, label="squaring", color=:red, lw=2)
plot!(sizes[1:length(times_naive)].^2, times_naive, color=:red, ls=:dash, label="", lw=2)

# vector multiplication

times_kron = []
times_naive = []

for s in sizes
    A = randn(s, s)
    B = rand(s, s)
    v = rand(s^2)
    K = A ⊗ B
    t = @belapsed $K * $v
    push!(times_kron, t)
    if s^2 < 5000
        K = collect(K)
        t = @belapsed $K * $v
        push!(times_naive, t)
        compute_naive = t < tmax
    end
end

plot!(sizes.^2, times_kron, label="matrix-vector mult.", color=:orange, lw=2)
plot!(sizes[1:length(times_naive)].^2, times_naive, color=:orange, ls=:dash, label="", lw=2)

yaxis!(:log10)
ylabel!("CPU time (s)")
xaxis!(:log10)
xlabel!("Kronecker product size")
title!("Performance Kronecker.jl (-) vs. native code (--)")

savefig("benchmark.svg")


text:
using Test

@testset "BridgeOptimizer" begin
    include("bridgeoptimizer.jl")
end
@testset "LazyBridgeOptimizer" begin
    include("lazybridgeoptimizer.jl")
end
@testset "Separate bridges" begin
    include("flip_sign_bridge.jl")
    include("vectorizebridge.jl")
    include("scalarizebridge.jl")
    include("slackbridge.jl")
    include("functionize_bridge.jl")
    include("intervalbridge.jl")
    include("rsocbridge.jl")
    include("quadtosocbridge.jl")
    include("geomeanbridge.jl")
    include("square_bridge.jl")
    include("detbridge.jl")
    include("soctopsdbridge.jl")
    include("indicator_bridge.jl")
end


text:
module gridbox_

using IUP

export
    gridbox

#------------------------example html/examples/C/gridbox.c ----------------------------------
function gridbox()
    IupOpen()                #Initializes IUP
    GridBoxTest()
    IupMainLoop()           # Initializes IUP main loop
    IupClose()              # And close it when ready
end

# ------------------------------------------------------------------------------------------
function GridBoxTest()

    gbox = IupGridBox(
        IupSetAttributes(IupLabel(""), ""),
        IupSetAttributes(IupLabel("col1"), "FONTSTYLE=Bold"),
        IupSetAttributes(IupLabel("col2"), "FONTSTYLE=Bold"),

        IupSetAttributes(IupLabel("lin1"), "FONTSTYLE=Bold"),
        IupSetAttributes(IupLabel("lbl"), "XSIZE=50x12"),
        IupSetAttributes(IupButton("but"), "XSIZE=50"),

        IupSetAttributes(IupLabel("lin2"), "FONTSTYLE=Bold"),
        IupSetAttributes(IupLabel("label"), "XSIZE=x12"),
        IupSetAttributes(IupButton("button"), "XEXPAND=Horizontal"),

        IupSetAttributes(IupLabel("lin3"), "FONTSTYLE=Bold"),
        IupSetAttributes(IupLabel("label large"), "XSIZE=x12"),
        IupSetAttributes(IupButton("button large"), "")
    )
    fr1 = IupFrame(gbox)

    IupSetAttribute(gbox, "SIZECOL", "2")
    IupSetAttribute(gbox, "SIZELIN", "3")
    #IupSetAttribute(gbox, "HOMOGENEOUSLIN", "Yes");
    #IupSetAttribute(gbox, "HOMOGENEOUSCOL", "Yes");
    #IupSetAttribute(gbox, "EXPANDCHILDREN", "HORIZONTAL");
    #IupSetAttribute(gbox, "NORMALIZESIZE", "BOTH");
    #IupSetAttribute(gbox, "NORMALIZESIZE", "HORIZONTAL");

    IupSetAttribute(gbox, "NUMDIV", "3")
    #IupSetAttribute(gbox, "ORIENTATION", "VERTICAL")
    #IupSetAttribute(gbox, "NUMDIV", "2")
    #IupSetAttribute(gbox, "NUMDIV", "AUTO")

    IupSetAttribute(gbox, "ALIGNMENTLIN", "ACENTER")
    IupSetAttribute(gbox, "MARGIN", "10x10")
    IupSetAttribute(gbox, "GAPLIN", "5")
    IupSetAttribute(gbox, "GAPCOL", "5")

    dlg = IupDialog(IupHbox(fr1))

    IupSetAttribute(dlg, "TITLE", "IupGridBox Test")
    IupSetAttribute(dlg, "MARGIN", "10x10")
    IupSetAttribute(fr1, "MARGIN", "0x0")       # avoid attribute propagation */

    IupShowXY(dlg, IUP_CENTER, IUP_CENTER)      # Shows dlg in the center of the screen
end

end     # module


text:
module GaOperators

import MachineGa.Chromosome
import MachineGa.FloatOperators

function onepointcrossover(c1::Chromosome, c2::Chromosome)::Chromosome
    Chromosome(FloatOperators.onepointcrossover(c1.genes, c2.genes), Inf64)
end

function uniformcrossover(c1::Chromosome, c2::Chromosome)::Chromosome
    Chromosome(FloatOperators.uniformcrossover(c1.genes, c2.genes), Inf64)
end

function randombytemutation(c::Chromosome; mutateprob::Float64 = 0.10)::Chromosome
    Chromosome(FloatOperators.randombytemutation(c.genes, mutateprob = mutateprob), Inf64)
end

function nearbytemutation(c::Chromosome; mutateprob::Float64 = 0.10)::Chromosome
    Chromosome(FloatOperators.nearbytemutation(c.genes, mutateprob = mutateprob), Inf64)
end

function tournamentselection(cs::Array{Chromosome,1}; tournaments::Int = 2)::Chromosome
    csample = rand(cs, tournaments)
    bestindex = sortperm(cs, by = x -> x.cost) |> first
    best = cs[bestindex]
    return Chromosome(best.genes, best.cost)
end

end # module


text:
isdefined(:JuliaParser) && workspace()
using JuliaParser
using FactCheck

include("lexer.jl")
include("parser.jl")
include("diagnostics.jl")
include("LineNumbers.jl")
#include("allsrc.jl")

exitstatus()


text:
using FractalTools
using Makie

# Contruct a tesslation
tridln = Node(TriDelaunay([0., 0.], [1., 0], [0.5, 1.], addboundarypoints=true))

# Definea an observables
msh = lift(tridln) do val
    tomesh(val.delaunay)
end
vtx = lift(msh) do val
    val.position
end

# Test point
tpnt = Node(getpoint(tridln[]))
tripts = lift(tpnt) do val
    tidx = simplices(tridln[])[locate(tridln[], val), :]
    map(val -> Point(val...), eachrow(points(tridln[])[tidx, :]))
end

# Plot tesselation
fig = Figure(resolution=(500, 500), scale_plot = false, camera = Makie.campixel!)
ax = fig[1, 1] =  Axis(fig)
mesh!(ax, msh, color=:lightblue)
wireframe!(ax, msh, color=:black)
scatter!(ax, vtx, color=:black)
scatter!(ax, tripts, color=:red, markersize=10)
scatter!(ax, tpnt, color=:orange, markersize=10)

# Add mouse button interaction
on(ax.scene.events.mousebuttons) do buttons
   if ispressed(ax.scene, Mouse.left)
        area = pixelarea(ax.scene)[]
        mp = ax.scene.events.mouseposition[]
        mp = Point2f0(mp) .- minimum(area)
        pos = to_world(ax.scene, mp)
        @info "$pos added"
        tpnt[] = pos
   end
end

# Add button interaction.
on(ax.scene.events.keyboardbuttons) do button
    if ispressed(button, Keyboard.a)
        pnt = addpoint!(tridln[])
        @info "$pnt is added"
    end
    if ispressed(button, Keyboard.f)
        finegrain!(tridln[], 10)
        @info "delaunay is fine grained"
    end
    tridln[] = tridln[]   # Trigger change
end

# Display scene
display(fig)


text:
# state transition:
#   initialized => inprogress => finished => finalized

primitive type State 8 end

const initialized = reinterpret(State, 0x00)
const inprogress  = reinterpret(State, 0x01)
const finished    = reinterpret(State, 0x02)
const finalized   = reinterpret(State, 0x03)
# const errored?

function Base.show(io::IO, s::State)
    if s == initialized
        print(io, "initialized")
    elseif s == inprogress
        print(io, "inprogress")
    elseif s == finished
        print(io, "finished")
    elseif s == finalized
        print(io, "finalized")
    else
        @assert false
    end
end

# do state transition
macro trans(obj, ts)
    if ts.head == :call && ts.args[1] == :(=>)
        ts = :($(ts),)
    end
    @assert ts.head == :tuple
    if VERSION < v"0.7.0-beta.66" # changed calling convention for foldr
        foldr(:(error("invalid state: ", $(esc(obj)).state)), ts.args) do t, elblk
            @assert t.head == :call && t.args[1] == :(=>)
            from, to = t.args[2], t.args[3]
            quote
                if $(esc(obj)).state == $(esc(from))
                    $(esc(obj)).state = $(esc(to))
                else
                    $(elblk)
                end
            end
        end
    else
        foldr(ts.args, init = :(error("invalid state: ", $(esc(obj)).state))) do t, elblk
            @assert t.head == :call && t.args[1] == :(=>)
            from, to = t.args[2], t.args[3]
            quote
                if $(esc(obj)).state == $(esc(from))
                    $(esc(obj)).state = $(esc(to))
                else
                    $(elblk)
                end
            end
        end
    end
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Entropy and KLD" begin

    vtree_file = zoo_vtree_file("simple2.vtree")

    pc1 = read((zoo_psdd_file("simple2.1.psdd"), vtree_file), StructProbCircuit)
    pc2 = read((zoo_psdd_file("simple2.2.psdd"), vtree_file), StructProbCircuit)
    pc3 = read((zoo_psdd_file("simple2.3.psdd"), vtree_file), StructProbCircuit)

    @test entropy(pc1) ≈ 1.2899219826090118
    @test entropy(pc2) ≈ 0.9359472745536583

    @test kl_divergence(pc1, pc2) ≈ 0.5672800167911778
    @test kl_divergence(pc2, pc3) ≈ 0.38966506

end


text:
@doc doc"""
    Grassmann{n,k,F} <: Manifold

The Grassmann manifold $\operatorname{Gr}(n,k)$ consists of all subspaces spanned
by $k$ linear independent vectors $\mathbb F^n$, where
$\mathbb F \in \{\mathbb R, \mathbb C\}$ is either the real- (or complex-) valued vectors.
This yields all $k$-dimensional subspaces of $\mathbb R^n$ for the real-valued case and all
$2k$-dimensional subspaces of $\mathbb C^n$ for the second.

The manifold can be represented as

````math
\operatorname{Gr}(n,k) \coloneqq \bigl\{ \operatorname{span}(x)
: x \in \mathbb F^{n\times k}, \bar{x}^\mathrm{T}x = I_k\},
````

where ${\bar\cdot}^{\mathrm{T}}$ denotes the complex conjugate transpose and
$I_k$ is the $k\times k$ identity matrix. This means, that the columns of $x$
form an orthonormal basis of the subspace, that is a point on
$\operatorname{Gr}(n,k)$, and hence the subspace can actually be represented by
a whole equivalence class of representers.
Another interpretation is, that

````math
\operatorname{Gr}(n,k) = \operatorname{St}(n,k) / \operatorname{O}(k),
````

i.e the Grassmann manifold is the quotient of the [`Stiefel`](@ref) manifold and
the orthogonal group $\operatorname{O}(k)$ of orthogonal $k\times k$ matrices.

The tangent space at a point (subspace) $x$ is given by

````math
T_x\mathrm{Gr}(n,k) = \bigl\{
v \in \mathbb{F}^{n\times k} :
{\bar v}^{\mathrm{T}}x + {\bar x}^{\mathrm{T}}v = 0_{k} \bigr\},
````

where $0_{k}$ denotes the $k\times k$ zero matrix.

Note that a point $x\in \operatorname{Gr}(n,k)$ might be represented by
different matrices (i.e. matrices with orthonormal column vectors that span
the same subspace). Different representations of $x$ also lead to different
representation matrices for the tangent space $T_x\mathrm{Gr}(n,k)$

The manifold is named after
[Hermann G. Graßmann](https://en.wikipedia.org/wiki/Hermann_Grassmann) (1809-1877).

# Constructor

    Grassmann(n,k,F=ℝ)

Generate the Grassmann manifold $\operatorname{Gr}(n,k)$, where the real-valued
case $\mathbb F = \mathbb R$ is the default.
"""
struct Grassmann{n,k,F} <: Manifold end

Grassmann(n::Int, k::Int, F::AbstractNumbers = ℝ) = Grassmann{n,k,F}()

@doc doc"""
    check_manifold_point(M::Grassmann{n,k,F}, x)

Check whether `x` is representing a point on the [`Grassmann`](@ref) `M`, i.e. its
a `n`-by-`k` matrix of unitary column vectors and of correct `eltype` with respect to `F`.
"""
function check_manifold_point(M::Grassmann{n,k,F}, x; kwargs...) where {n,k,F}
    if (F === ℝ) && !(eltype(x) <: Real)
        return DomainError(
            eltype(x),
            "The matrix $(x) is not a real-valued matrix, so it does noe lie on the Grassmann manifold of dimension ($(n),$(k)).",
        )
    end
    if (F === ℂ) && !(eltype(x) <: Real) && !(eltype(x) <: Complex)
        return DomainError(
            eltype(x),
            "The matrix $(x) is neiter real- nor complex-valued matrix, so it does noe lie on the complex Grassmann manifold of dimension ($(n),$(k)).",
        )
    end
    if size(x) != representation_size(M)
        return DomainError(
            size(x),
            "The matrix $(x) is does not lie on the Grassmann manifold of dimension ($(n),$(k)), since its dimensions are wrong.",
        )
    end
    c = x' * x
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(x) does not lie on the Grassmann manifold of dimension ($(n),$(k)), because x'x is not the unit matrix.",
        )
    end
end

@doc doc"""
    check_tangent_vector(M::Grassmann{n,k,F}, x, v)

Check whether `v` is a tangent vector in the tangent space of `x` on the [`Grassmann`](@ref)
`M`, i.e. that `v` is of size and type as well as that

````math
    x^{\mathrm{H}}v + v^{\mathrm{H}}x = 0_k,
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian and $0_k$
denotes the $k\times k$ zero natrix.
"""
function check_tangent_vector(G::Grassmann{n,k,F}, x, v; kwargs...) where {n,k,F}
    t = check_manifold_point(G, x)
    t === nothing || return t
    if (F === ℝ) && !(eltype(v) <: Real)
        return DomainError(
            eltype(v),
            "The matrix $(v) is not a real-valued matrix, so it can not be a tangent vector to the Grassmann manifold of dimension ($(n),$(k)).",
        )
    end
    if (F === ℂ) && !(eltype(v) <: Real) && !(eltype(v) <: Complex)
        return DomainError(
            eltype(v),
            "The matrix $(v) is neiter real- nor complex-valued matrix, so it can not bea tangent vector to the complex Grassmann manifold of dimension ($(n),$(k)).",
        )
    end
    if size(v) != representation_size(G)
        return DomainError(
            size(v),
            "The matrix $(v) is does not lie in the tangent space of $(x) on the Grassmann manifold of dimension ($(n),$(k)), since its dimensions are wrong.",
        )
    end
    if !isapprox(x' * v + v' * x, zeros(k, k); kwargs...)
        return DomainError(
            norm(x' * v + v' * x),
            "The matrix $(v) is does not lie in the tangent space of $(x) on the Grassmann manifold of dimension ($(n),$(k)), since x'v + v'x is not the zero matrix.",
        )
    end
end

@doc doc"""
    distance(M::Grassmann, x, y)

Compute the Riemannian distance on [`Grassmann`](@ref) manifold `M`$= \mathrm{Gr}(n,k)$.

Let $USV = x^\mathrm{H}y$ denote the SVD decomposition of
$x^\mathrm{H}y$, where $\cdot^{\mathrm{H}}$ denotes the complex
conjugate transposed or Hermitian. Then the distance is given by
````math
d_{\mathrm{GR}(n,k)}(x,y) = \operatorname{norm}(\operatorname{Re}(b)).
````
where

$b_{i}=\begin{cases} 0 & \text{if} \; S_i \geq 1\\ \operatorname{acos}(S_i) & \, \text{if} \; S_i<1 \end{cases}.$
"""
function distance(M::Grassmann, x, y)
    x ≈ y && return zero(real(eltype(x)))
    a = svd(x' * y).S
    a[a.>1] .= 1
    return sqrt(sum((acos.(a)) .^ 2))
end

@doc doc"""
    exp(M::Grassmann, x, v)

Compute the exponential map on the [`Grassmann`](@ref) `M`$= \mathrm{Gr}(n,k)$ starting in
`x` with tangent vector (direction) `v`. Let $v = USV$ denote the SVD decomposition of $v$.
Then the exponential map is written using

````math
z = x V\cos(S)V^\mathrm{H} + U\sin(S)V^\mathrm{H},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
The cosine and sine are applied element wise to the diagonal entries of $S$.
A final QR decomposition $z=QR$ is performed for numerical stability reasons,
yielding the result as
````math
\exp_x v = Q.
````
"""
exp(::Grassmann, ::Any...)

function exp!(M::Grassmann, y, x, v)
    norm(M, x, v) ≈ 0 && return copyto!(y, x)
    d = svd(v)
    z = x * d.V * Diagonal(cos.(d.S)) * d.Vt + d.U * Diagonal(sin.(d.S)) * d.Vt
    # reorthonormalize
    return copyto!(y, Array(qr(z).Q))
end

injectivity_radius(::Grassmann) = π / 2

@doc doc"""
    inner(M::Grassmann, x, v, w)

Compute the inner product for two tangent vectors `v`, `w` from the
tangent space of `x` on the [`Grassmann`](@ref) manifold `M`.
The formula reads

````math
g_x(v,w) = \operatorname{trace}(v^{\mathrm{H}}w),
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
"""
inner(::Grassmann, x, v, w) = dot(v, w)

@doc doc"""
    inverse_retract(M::Grassmann, x, y, ::PolarInverseRetraction)

Compute the inverse retraction for the [`PolarRetraction`](@ref), on the
[`Grassmann`](@ref), i.e.,

````math
\operatorname{retr}_x^{-1}y = y*(x^\mathrm{H}y)^{-1} - x,
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
"""
inverse_retract(M::Grassmann, ::Any, ::Any, ::PolarInverseRetraction)

function inverse_retract!(::Grassmann, v, x, y, ::PolarInverseRetraction)
    return copyto!(v, y / (x' * y) - x)
end

@doc doc"""
    inverse_retract(M, x, y, ::QRInverseRetraction)

Compute the inverse retraction valid of the [`QRRetraction`](@ref)

````math
\operatorname{retr}_x^{-1}y = y*(x^\mathrm{H}y)^{-1} - x,
````
where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
"""
inverse_retract(::Grassmann, ::Any, ::Any, ::QRInverseRetraction)

inverse_retract!(::Grassmann, v, x, y, ::QRInverseRetraction) = copyto!(v, y / (x' * y) - x)

function isapprox(M::Grassmann, x, v, w; kwargs...)
    return isapprox(sqrt(inner(M, x, zero_tangent_vector(M, x), v - w)), 0; kwargs...)
end
isapprox(M::Grassmann, x, y; kwargs...) = isapprox(distance(M, x, y), 0.0; kwargs...)

@doc doc"""
    log(M::Grassmann, x, y)

Compute the logarithmic map on the [`Grassmann`](@ref) `M`$ = \mathcal M=\mathrm{Gr}(n,k)$,
i.e. the tangent vector `v` whose corresponding [`geodesic`](@ref) starting from `x`
reaches `y` after time 1 on `M`. The formula reads

````math
\log_xy = V\cdot \operatorname{atan}(S) \cdot U^\mathrm{H},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
$U$ and $V$ are the unitary matrices, and $S$ is a diagonal matrix containing
the singular values of the SVD-decomposition of
````math
USV = (y^\mathrm{H}x)^{-1} ( y^\mathrm{H} - y^\mathrm{H}xx^\mathrm{H} ).
````
In this formula the $\operatorname{atan}$ is meant elementwise.
"""
log(::Grassmann, ::Any...)

function log!(M::Grassmann, v, x, y)
    z = y' * x
    At = y' - z * x'
    Bt = z \ At
    d = svd(Bt')
    return copyto!(v, d.U * Diagonal(atan.(d.S)) * d.Vt)
end

@doc doc"""
    manifold_dimension(M::Grassmann)

Return the dimension of the [`Grassmann(n,k,𝔽)`](@ref) manifold `M`, i.e.

````math
\dim \operatorname{Gr}(n,k) = k(n-k) \dim_ℝ 𝔽,
````

where $\dim_ℝ 𝔽$ is the [`real_dimension`](@ref) of `𝔽`.
"""
manifold_dimension(M::Grassmann{n,k,𝔽}) where {n,k,𝔽} = k * (n - k) * real_dimension(𝔽)

"""
    mean(
        M::Grassmann,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/4);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::Grassmann{n,k,ℝ} where {n,k}, ::Any...)

function mean!(
    M::Grassmann{n,k,ℝ},
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
) where {n,k}
    return mean!(M, y, x, w, GeodesicInterpolationWithinRadius(π / 4); kwargs...)
end

@doc doc"""
    project_tangent(M::Grassmann, x, w)

Project the `n`-by-`k` `w` onto the tangent space of `x` on the [`Grassmann`](@ref) `M`,
which is computed by

````math
\operatorname{proj_x}(w) = w - xx^{\mathrm{H}}w,
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
"""
project_tangent(::Grassmann, ::Any...)

project_tangent!(M::Grassmann, v, x, w) = copyto!(v, w - x * x' * w)

@doc doc"""
    representation_size(M::Grassmann{n,k,F})

Return the represenation size or matrix dimension of a point on the [`Grassmann`](@ref)
`M`, i.e. $(n,k)$ for both the real-valued and the complex value case.
"""
@generated representation_size(::Grassmann{n,k}) where {n,k} = (n, k)

@doc doc"""
    retract(M::Grassmann, x, v, ::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref) on the
[`Grassmann`](@ref) `M`. With $USV = x + v$ the retraction reads
````math
\operatorname{retr}_x v = UV^\mathrm{H},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

    retract(M::Grassmann, x, v, ::QRRetraction )

Compute the QR-based retraction [`QRRetraction`](@ref) on the
[`Grassmann`](@ref) `M`. With $QR = x + v$ the retraction reads
````math
\operatorname{retr}_xv = QD,
````
where D is a $m\times n$ matrix with
````math
D = \operatorname{diag}( \operatorname{sgn}(R_{ii}+0,5)_{i=1}^n ).
````
"""
retract(::Grassmann, ::Any...)

function retract!(::Grassmann, y, x, v, ::PolarRetraction)
    s = svd(x + v)
    return mul!(y, s.U, s.Vt)
end
function retract!(::Grassmann{N,K}, y, x, v, ::QRRetraction) where {N,K}
    qrfac = qr(x + v)
    d = diag(qrfac.R)
    D = Diagonal(sign.(sign.(d .+ 0.5)))
    y .= zeros(N, K)
    y[1:K, 1:K] .= D
    return copyto!(y, Array(qrfac.Q) * D)
end

@doc doc"""
    zero_tangent_vector(M::Grassmann, x)

Return the zero tangent vector from the tangent space at `x` on the [`Grassmann`](@ref) `M`,
which is given by a zero matrix the same size as `x`.
"""
zero_tangent_vector(::Grassmann, ::Any...)

zero_tangent_vector!(::Grassmann, v, x) = fill!(v, 0)


text:
import Pkg
Pkg.add("Documenter")
using ComplexMixtures
using Documenter
using ComplexMixtures
push!(LOAD_PATH,"../src/")
makedocs(
    modules=[ComplexMixtures],
    sitename="ComplexMixtures.jl",
    pages = [
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Quick Guide" => "quickguide.md",
        "Full Example" => "examples.md",
        "Set solute and solvent" => "selection.md",
        "Loading the trajectory" => "trajectory.md",
        "Computing the MDDF" => "mddf.md",
        "Results" => "results.md",
        "Atomic and group contributions" => "contrib.md",
        "Save and load" => "save.md",
        "Multiple trajectories" => "multiple.md",
        "Parallel execution" => "parallel.md",
        "Options" => "options.md",
        "Tools" => "tools.md",
        "Help entries" => "help.md",
        "References" => "references.md"
    ]
)
deploydocs(
    repo = "github.com/m3g/ComplexMixtures.jl.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#" ],
)


text:
include("JuliaOpenGL.jl")
JuliaOpenGL.main()

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
## -------------------------------------------------------------------
let
    p = plot()
    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
        # (ChH.iJR904, ChH.HeerdenData),
    ]
        src = nameof(Data)
        marker = (6, source_markers[Data])
        scatter!(p, Data.val(:D), abs.(Data.val(:D) ./ Data.val(:uGLC));
            label = "", marker
        )
    end
    mysavefig(p, "exp_yield_vs_D_coor")
end


text:
#
# Provides default transforms, and the means for users to create their own.
#

"""
$(TYPEDEF)

A `supertype` for all astrodynamics coordinate frame transformations.
"""
abstract type OrbitalFrameTransform end # should this be a subtype of CoordinateTransformations.Transformation?

"""
$(TYPEDEF)

A generic frame transformation between two `OrbitalFrame` types.
Converts `F1` to`F2`.

$(TYPEDFIELDS)
"""
struct Transform{F1<:OrbitalFrame, F2<:OrbitalFrame, T1<:CoordinateTransformations.Transformation, T2<:CoordinateTransformations.Transformation} <: OrbitalFrameTransform
    transform_position::T1
    transform_velocity::T2
end

"""
$(SIGNATURES)

An outer constructor for `Transform` instances.
Provide position and velocity transformations,
and the initial and final reference frames via
function arguments. This is an alternate syntax
for `Transform{InitialFrame, FinalFrame}(position_transform, velocity_transform)`.
"""
function Transform(position_transform::CoordinateTransformations.Transformation, velocity_transform::CoordinateTransformations.Transformation, from::F1, to::F2) where {F1 <: OrbitalFrameTransform, F2 <: OrbitalFrameTransform}
    return Transform{F1, F2}(position_transform, velocity_transform)
end


text:
# Extension of some common Array operations to Tuple objects

import Base: isapprox, similar, zero

isapprox(x::Tuple, y::Tuple) = all(isapprox.(x, y))

similar(x::Tuple) = similar.(x)

zero(x::Tuple) = zero.(x)


text:
#' # [MNIST: Elastic Distortions](@id elastic)

#' In this example we are going to use Augmentor on the famous
#' **MNIST database of handwritten digits** [^MNIST1998] to
#' reproduce the elastic distortions discussed in [^SIMARD2003].
#' It may be interesting to point out, that the way Augmentor
#' implements distortions is a little different to how it is
#' described by the authors of the paper.
#' This is for a couple of reasons, most notably that we want the
#' parameters for our deformations to be independent of the size
#' of image it is applied on. As a consequence the
#' parameter-numbers specified in the paper are not 1-to-1
#' transferable to Augmentor.

#' If the effects are sensible for the dataset, then applying
#' elastic distortions can be a really effective way to improve
#' the generalization ability of the network.
#' That said, our implementation of [`ElasticDistortion`](@ref)
#' has a lot of possible parameters to choose from. To that end,
#' we will introduce a simple strategy for interactively
#' exploring the parameter space on our dataset of interest.

#md #' !!! note
#md #'
#md #'     This tutorial was designed to be performed in a
#md #'     [Juypter](https://jupyter.org/) notebook. You can
#md #'     find a link to the Juypter version of this tutorial
#md #'     in the top right corner of this page.

#' ## Loading the MNIST Trainingset

#' In order to access and visualize the MNIST images we employ
#' the help of two additional Julia packages. In the interest of
#' time and space we will not go into great detail about their
#' functionality. Feel free to click on their respective names to
#' find out more information about the utility they can provide.
#'
#' - [Images.jl](https://github.com/JuliaImages/Images.jl) will
#'   provide us with the necessary tools for working with image
#'   data in Julia.
#'
#' - [MLDatasets.jl](https://github.com/JuliaML/MLDatasets.jl)
#'   has an MNIST submodule that offers a convenience interface
#'   to read the MNIST database.

#' The function `MNIST.traintensor` returns the MNIST training
#' images corresponding to the given indices as a
#' multi-dimensional array. These images are stored in the native
#' horizontal-major memory layout as a single floating point
#' array, where all values are scaled to be between 0.0 and 1.0.

using Images, MLDatasets
train_tensor = MNIST.traintensor()
@show summary(train_tensor);
#md nothing # hide

#' This horizontal-major format is the standard way of utilizing
#' this dataset for training machine learning models.
#' In this tutorial, however, we are more interested in working
#' with the MNIST images as actual Julia images in vertical-major
#' layout, and as black digits on white background.

#' We can convert the "tensor" to a `Colorant` array using the
#' provided function `MNIST.convert2image`.
#' This way, Julia knows we are dealing with image data and can
#' tell programming environments such as Juypter how to visualize
#' it. If you are working in the terminal you may want to use the
#' package [ImageInTerminal.jl](https://github.com/JuliaImages/ImageInTerminal.jl)

train_images = MNIST.convert2image(train_tensor)
img_1 = train_images[:,:,1] # show first image
#md save("mnist_1.png",repeat(img_1,inner=(4,4))) # hide
#md nothing # hide

#md #' ![first image](mnist_1.png)

#' ## Visualizing the Effects

#' Before applying an operation (or pipeline of operations) on
#' some dataset to train a network, we strongly recommend
#' investing some time in selecting a decent set of hyper
#' parameters for the operation(s). A useful tool for tasks like
#' this is the package [Interact.jl](https://github.com/JuliaGizmos/Interact.jl).
#' We will use this package to define a number of widgets for
#' controlling the parameters to our operation.

#' Note that while the code below only focuses on configuring
#' the parameters of a single operation, specifically
#' [`ElasticDistortion`](@ref), it could also be adapted to tweak
#' a whole pipeline. Take a look at the corresponding section in
#' [High-level Interface](@ref pipeline) for more information
#' on how to define and use a pipeline.

# These two package will provide us with the capabilities
# to perform interactive visualisations in a jupyter notebook
using Augmentor, Interact, Reactive

# The manipulate macro will turn the parameters of the
# loop into interactive widgets.
@manipulate for
        unpaused = true,
        #ticks = fpswhen(Signal(unpaused), 5.), # non-boolean ... error
        image_index = 1:100,
        grid_size = 3:20,
        scale = .1:.1:.5,
        sigma = 1:5,
        iterations = 1:6,
        free_border = true
    op = ElasticDistortion(grid_size, grid_size, # equal width & height
                           sigma = sigma,
                           scale = scale,
                           iter = iterations,
                           border = free_border)
    augment(train_images[:, :, image_index], op)
end
#md nothing # hide

#md #' Executing the code above in a Juypter notebook will result
#md #' in the following interactive visualisation. You can now
#md #' use the sliders to investigate the effects that different
#md #' parameters have on the MNIST training images.
#md #'
#md #' !!! tip
#md #'
#md #'     You should always use your **training** set to do this
#md #'     kind of visualisation (not the test test!). Otherwise
#md #'     you are likely to achieve overly optimistic (i.e. biased)
#md #'     results during training.
#md #'
#md #' ![interact](https://user-images.githubusercontent.com/10854026/30867456-4afe0800-a2dc-11e7-90eb-800b6ea025d0.gif)

#' Congratulations! With just a few simple lines of code, you
#' created a simple interactive tool to visualize your image
#' augmentation pipeline. Once you found a set of parameters that
#' you think are appropriate for your dataset you can go ahead
#' and train your model.

#' ## References
#'
#' [^MNIST1998]: LeCun, Yan, Corinna Cortes, Christopher J.C. Burges. ["The MNIST database of handwritten digits"](http://yann.lecun.com/exdb/mnist/) Website. 1998.
#'
#' [^SIMARD2003]: Simard, Patrice Y., David Steinkraus, and John C. Platt. ["Best practices for convolutional neural networks applied to visual document analysis."](https://www.microsoft.com/en-us/research/publication/best-practices-for-convolutional-neural-networks-applied-to-visual-document-analysis/) ICDAR. Vol. 3. 2003.


text:
using TcpInstruments
using Test

pwr = initialize(VersatilePower)
@info "Successfully connected $(pwr.model) at $(pwr.address)"

"""
Spec:
enable_output()
disable_output()
set_voltage(x)
get_voltage()
set_current_limit(x)
get_current_limit()
remote_mode()
local_mode()
"""

remote_mode(pwr)
@testset "Output" begin
    @info get_output(pwr)

    enable_output(pwr)
    @test get_output(pwr) == true

    disable_output(pwr)
    @test get_output(pwr) == false
end

@testset "Current" begin
    function test_current(cur)
        set_current_limit(pwr, cur)
        @info get_current_limit(pwr)
        @test get_current_limit(pwr) == cur
    end

    test_current(0)
    test_current(1)
    test_current(1.5)
    test_current("1.6")
    test_current(1.72e0)
    test_current(2)
    test_current(3)
    test_current(4)

    set_current_limit(pwr, "MAX")
    @info get_current_limit(pwr)

    set_current_limit(pwr, "MIN")
    @info get_current_limit(pwr)

    test_current(5)
    test_current(0)
end

@testset "Voltage" begin
    @info get_voltage(pwr)
    set_voltage(pwr, 2)
    @info get_voltage(pwr)
    set_voltage(pwr, 4)
    @info get_voltage(pwr)
    set_voltage(pwr, 6)
    @info get_voltage(pwr)
    set_voltage(pwr, 8)
    @info get_voltage(pwr)
    set_voltage(pwr, 20)
    @info get_voltage(pwr)
    set_voltage(pwr, 25)
    @info get_voltage(pwr)
    set_voltage(pwr, 0)
    @info get_voltage(pwr)
    set_voltage(pwr, "MAX")
    @info get_voltage(pwr)
    set_voltage(pwr, "MIN")
    @info get_voltage(pwr)
end

local_mode(pwr)
terminate(pwr)

@info "Successfully disconnected"
@info "Goodbye"


text:
# ---
# title: 1312. Minimum Insertion Steps to Make a String Palindrome
# id: problem1312
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-insertion-steps-to-make-a-string-palindrome/description/>
# hidden: true
# ---
#
# Given a string `s`. In one step you can insert any character at any index of
# the string.
#
# Return _the minimum number of steps_ to make `s` palindrome.
#
# A  **Palindrome String**  is one that reads the same backward as well as
# forward.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "zzazz"
#     Output: 0
#     Explanation: The string "zzazz" is already palindrome we don't need any insertions.
#
#
# **Example 2:**
#
#
#
#     Input: s = "mbadm"
#     Output: 2
#     Explanation: String can be "mbdadbm" or "mdbabdm".
#
#
# **Example 3:**
#
#
#
#     Input: s = "leetcode"
#     Output: 5
#     Explanation: Inserting 5 characters the string becomes "leetcodocteel".
#
#
# **Example 4:**
#
#
#
#     Input: s = "g"
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: s = "no"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 500`
#   * All characters of `s` are lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
Return the width in miliseconds of the spike

  function spike_width(spike::Array{Float64,1},samplingrate::Real)
"""
function spike_width(spike::Array{Float64,1},samplingrate::Real)
	spikei = Grid.InterpGrid(spike,Grid.BCnil, Grid.InterpQuadratic)
	dt = length(spike)/samplingrate
	#upsample
	spikes = spikei[linspace(1,length(spike),1000)]
	idx = find(spikes .< 0.5*minimum(spikes))
	return dt*(idx[end]	- idx[1])
end

function spike_width(TF::TemplateFile,samplingrate::Real)
	#find the largest
	w = zeros(size(TF.templates,3))
	for i in 1:size(TF.templates,3)
		idx = indmin(TF.templates[:,:,i][:])
		pt,ch = ind2sub((size(TF.templates,1), size(TF.templates,2)),idx)
		w[i] = spike_width(TF.templates[:,ch,i],samplingrate)
	end
	w
end

spike_width(A) = spike_width(A, 40000)

function spike_dv_ratio(spike::Array{Float64,1},samplingrate::Real)
	spikei = Grid.InterpGrid(spike,Grid.BCnil, Grid.InterpQuadratic)
	dt = length(spike)/samplingrate
	#upsample
	spikes = spikei[linspace(1,length(spike),1000)]
	dv = diff(spikes)
	abs(maximum(dv)/minimum(dv))
end

spike_dv_ratio(spike::Array{Float64,1}) = spike_dv_ratio(spike,40000)

@compat function get_features(session::AbstractString)
	features = Array(Features,0)
	get_features!(features,session)
	features
end

@compat function get_features!(features::Array{Features,1},session::AbstractString)
	templatefiles = split(readchomp(`find . -name "$(session)*templates*.hdf5"`),"\n")
	return get_features!(features,templatefiles)
end

@compat function get_features{T<:AbstractString}(templatefiles::Array{T,1})
	features = Array(Features,0)
	get_features!(features,templatefiles)
	features
end

@compat function get_features!{T<:AbstractString}(features::Array{Features,1},templatefiles::Array{T,1})
	w = Float64[]
	isi = Float64[]
	for tf in templatefiles
		m = match(r"([[:alpha:][:digit:]_]*)_templatesg([[:digit:]]*)", tf)
		dd,pp = splitdir(tf)
		session = m.captures[1]
		group = m.captures[2]
		TF = get_valid_templates(tf)
		if TF == nothing
			continue
        end
		if !isempty(dd)
			f = "$(dd)/$(session)g$(group)Spiketrains.mat"
		else
			f = "$(session)g$(group)Spiketrains.mat"
		end
		if !isfile(f)
			continue
		end
		DD = MAT.matread(f)
		if length(DD) == size(TF.templates,3)
			_w = spike_width(TF)
			for c in 1:size(TF.templates,3)
				cc = @sprintf "%02d" c
				_isi = diff(DD["cluster$(cc)s"][:])
				spikes_in_bursts = (sum(_isi .< 5)+1)/length(_isi+1)
				m_isi = percentile(_isi,5)
				dvr = spike_dv_ratio(TF.templates[:,:,c][:])
				push!(features, Features(_w[c], spikes_in_bursts,m_isi,dvr))
			end
		end
	end
	return features
end