text:
import Base: Matrix
import Base: size
using ArnoldiMethod

#using LightGraphs.LinAlg: eigs
# just so that we can assert equality of matrices
Matrix(nbt::Nonbacktracking) = Matrix(sparse(nbt))
@testset "Spectral" begin

    g3 = PathGraph(5)
    g4 = PathDiGraph(5)
    g5 = SimpleDiGraph(4)
    for g in testgraphs(g3)
        @test adjacency_matrix(g, Bool) == adjacency_matrix(g, Bool; dir=:out)
        @test adjacency_matrix(g)[3, 2] == 1
        @test adjacency_matrix(g)[2, 4] == 0
        @test laplacian_matrix(g)[3, 2] == -1
        @test laplacian_matrix(g)[1, 3] == 0
        @test laplacian_spectrum(g)[5] ≈ 3.6180339887498945
        @test adjacency_spectrum(g)[1] ≈ -1.732050807568878
    end

    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)
    for g in testdigraphs(g5)
        @test adjacency_matrix(g, Bool) == adjacency_matrix(g, Bool; dir=:out)
        @test laplacian_spectrum(g)[3] == laplacian_spectrum(g; dir=:both)[3] == 3.0
        @test laplacian_spectrum(g; dir=:in)[3] == 1.0
        @test laplacian_spectrum(g; dir=:out)[3] == 1.0
    end

    # check adjacency matrices with self loops
    gx = copy(g3)
    add_edge!(gx, 1, 1)
    for g in testgraphs(gx)
        @test adjacency_matrix(g)[1, 1] == 2
    end
    g = copy(g5)
    add_edge!(g, 1, 1)
    @test adjacency_matrix(g)[1, 1] == 1
    @test indegree(g) == sum(adjacency_matrix(g), dims=1)[1, :]
    @test outdegree(g) == sum(adjacency_matrix(g), dims=2)[:, 1]

    g10 = CompleteGraph(10)
    for g in testgraphs(g10)
        B, em = non_backtracking_matrix(g)
        @test length(em) == 2 * ne(g)
        @test size(B) == (2 * ne(g), 2 * ne(g))
        for i = 1:10
            @test sum(B[:, i]) == 8
            @test sum(B[i, :]) == 8
        end
        @test !issymmetric(B)

        v = ones(Float64, ne(g))
        z = zeros(Float64, nv(g))
        n10 = Nonbacktracking(g)
        @test size(n10) == (2 * ne(g), 2 * ne(g))
        @test size(n10,1) == n10.m
        @test size(n10,2) == n10.m
        @test eltype(n10) == Float64
        @test !issymmetric(n10)

        contract!(z, n10, v)

        zprime = contract(n10, v)
        @test z == zprime
        @test z == 9 * ones(Float64, nv(g))
    end

    for g in testdigraphs(g5)
        @test adjacency_spectrum(g)[3] ≈ 0.311 atol = 0.001
    end

    for g in testgraphs(g3)
        @test adjacency_matrix(g) ==
            adjacency_matrix(g, dir=:out) ==
            adjacency_matrix(g, dir=:in) ==
            adjacency_matrix(g, dir=:both)

        @test_throws ErrorException adjacency_matrix(g; dir=:purple)
    end

    #that call signature works
    for g in testdigraphs(g5)
        inmat   = adjacency_matrix(g, Int; dir=:in)
        outmat  = adjacency_matrix(g, Int; dir=:out)
        bothmat = adjacency_matrix(g, Int; dir=:both)

        #relations that should be true
        @test inmat' == outmat
        @test all((bothmat - outmat) .>= 0)
        @test all((bothmat - inmat)  .>= 0)

      #check properties of the undirected laplacian carry over.
        for dir in [:in, :out, :both]
            T = eltype(g)
            amat = adjacency_matrix(g, Float64; dir=dir)
            lmat = laplacian_matrix(g, Float64; dir=dir)
            @test isa(amat, SparseMatrixCSC{Float64,T})
            @test isa(lmat, SparseMatrixCSC{Float64,T})
            evals = eigvals(Matrix(lmat))
            @test all(evals .>= -1e-15) # positive semidefinite
            @test (minimum(evals)) ≈ 0 atol = 1e-13
        end
    end

    for g in testdigraphs(g4)
        # testing incidence_matrix, first directed graph
        @test size(incidence_matrix(g)) == (5, 4)
        @test incidence_matrix(g)[1, 1] == -1
        @test incidence_matrix(g)[2, 1] == 1
        @test incidence_matrix(g)[3, 1] == 0
    end

    for g in testgraphs(g3)
        # now undirected graph
        @test size(incidence_matrix(g)) == (5, 4)
        @test incidence_matrix(g)[1, 1] == 1
        @test incidence_matrix(g)[2, 1] == 1
        @test incidence_matrix(g)[3, 1] == 0

        # undirected graph with orientation
        @test size(incidence_matrix(g; oriented=true)) == (5, 4)
        @test incidence_matrix(g; oriented=true)[1, 1] == -1
        @test incidence_matrix(g; oriented=true)[2, 1] == 1
        @test incidence_matrix(g; oriented=true)[3, 1] == 0
    end
    # TESTS FOR Nonbacktracking operator.

    n = 10; k = 5
    pg = CompleteGraph(n)
    # ϕ1 = nonbacktrack_embedding(pg, k)'
    for g in testgraphs(pg)
        nbt = Nonbacktracking(g)
        B, emap = non_backtracking_matrix(g)
        Bs = sparse(nbt)
        @test sparse(B) == Bs
        @test eigs(nbt, which=LR(), nev=1)[1] ≈ eigs(B, which=LR(), nev=1)[1] atol = 1e-5

        # check that matvec works
        x = ones(Float64, nbt.m)
        y = nbt * x
        z = B * x
        @test norm(y - z) < 1e-8

        #check that matmat works and Matrix(nbt) == B

        @test norm(nbt * Matrix{Float64}(I, nbt.m, nbt.m) - B) < 1e-8

        #check that matmat works and Matrix(nbt) == B
        @test norm(nbt * Matrix{Float64}(I, nbt.m, nbt.m) - B) < 1e-8

        @test size(y) == size(x)

        B₁ = Nonbacktracking(g10)

        @test Matrix(B₁) == Matrix(B)
        @test  B₁ * ones(size(B₁)[2]) == B * ones(size(B)[2])
        @test size(B₁) == size(B)
        @test !issymmetric(B₁)
        @test eltype(B₁) == Float64
    end
    # END tests for Nonbacktracking

    # spectral distance checks
    for n = 3:10
        polygon = random_regular_graph(n, 2)
        for g in testgraphs(polygon)
            @test spectral_distance(g, g) ≈ 0 atol = 1e-8
            @test spectral_distance(g, g, 1) ≈ 0 atol = 1e-8
        end
    end
end


text:
# noticed something similar in Boost datetime

struct WrappingPaper{T} <: Signed
    minval::T  # minval - 1 == maxval
    maxval::T  # maxval + 1 == minval
    nvals::T   # maxval - minval + 1
end

function WrappingPaper(minval, maxval)
    minval, maxval = minmax(minval, maxval)
    nvals = maxval - minval + 1
    WrappingPaper(minval, maxval, nvals)
end

WrappingPaper(x::WrappingPaper) = x

function Base.show(io::IO, x::WrappingPaper)
    nt = (minval=x.minval, maxval=x.maxval, nvals=x.nvals)
    print(io, string(nt))
end

months28days = WrappingPaper(1, 28)
months29days = WrappingPaper(1, 29)
months30days = WrappingPaper(1, 30)
months31days = WrappingPaper(1, 31)

mutable struct WrappingInt{T} <: Signed
    x::T
    minval::T
    maxval::T
end

function WrappingInt(x::T, wp::WrappingPaper{T}) where {T}
    WrappingInt(x, wp.minval, wp.maxval)
end

ext: jl


text:
using CMAEvolutionStrategy
using Test, Statistics, LibGit2, PyCall
# using BlackBoxOptimizationBenchmarking
# const BBOB = BlackBoxOptimizationBenchmarking

const PYCMA_PATH = joinpath(@__DIR__, "pycma")
if !isdir(PYCMA_PATH)
    const PYCMAREPO = LibGit2.clone("https://github.com/CMA-ES/pycma", PYCMA_PATH)
    LibGit2.checkout!(PYCMAREPO, "025ef1fed91c86690a21e9ed81713062d29398ff")
end
pushfirst!(PyVector(pyimport("sys")."path"), PYCMA_PATH)
cma = pyimport("cma")

# TODO: reactivate these tests when BBOB registers new version
# https://github.com/jonathanBieler/BlackBoxOptimizationBenchmarking.jl/issues/9
# @testset "BBOB comparison with PyCMA" begin
#     pinit(D) = 10*rand(D) .- 5
#
#     struct CMAES end
#     function BBOB.optimize(::CMAES, f, D, run_length)
#         CMAEvolutionStrategy.minimize(f, pinit(D), 1., verbosity = 0, maxfevals = run_length)
#     end
#     BBOB.minimum(o::CMAEvolutionStrategy.Optimizer) = fbest(o)
#     BBOB.minimizer(o::CMAEvolutionStrategy.Optimizer) = xbest(o)
#     struct PyCMA end
#     function BBOB.optimize(m::PyCMA, f, D, run_length)
#         es = cma.CMAEvolutionStrategy(pinit(D), 1, cma.CMAOptions(verb_log = 0,
#                                                                   verb_disp = 0,
#                                                                   maxfevals = run_length))
#         mfit = es.optimize(f).result
#         (m, mfit[1], mfit[2])
#     end
#     BBOB.minimum(mfit::Tuple{PyCMA,Vector{Float64},Float64}) = mfit[3]
#     BBOB.minimizer(mfit::Tuple{PyCMA,Vector{Float64},Float64}) = mfit[2]
#
#     D = [3, 12]
#     lengths = [1_000, 20_000]
#     res1 = BBOB.benchmark(CMAES(), BBOB.BBOBFunctions, lengths, 10, D, 1e-6)
#     res2 = BBOB.benchmark(PyCMA(), BBOB.BBOBFunctions, lengths, 10, D, 1e-6)
#
#     @test sum(res1.success_counts) > .99 * sum(res2.success_counts)
#     @test sum(res1.minimum) < 1.1 * sum(res2.minimum)
# end

@testset "BoxConstaints" begin
    respycma = cma.fmin(cma.ff.rosen, zeros(3), .25,
                        cma.CMAOptions(bounds = [fill(-Inf, 3), fill(.5, 3)]))
    rescma = minimize(cma.ff.rosen, zeros(3), .25, upper = fill(.5, 3))
    @test respycma[1] ≈ xbest(rescma) atol = 1e-5
    @test respycma[2] ≈ fbest(rescma) atol = 1e-5

    import CMAEvolutionStrategy: BoxConstraints, transform, backtransform
    c = BoxConstraints(-1.2*ones(3), ones(3))
    x = [-1.3, 0, 1.1]
    @test (&)( (-1.2*ones(3) .<= transform(c, x) .<= ones(3))... )
    @test backtransform(c, transform(c, x)) ≈ x
    x = [-1.3 .4; 0 -1.2; 1.2 0]
    @test transform(c, x) == hcat([transform(c, x[:, i]) for i in 1:2]...)
end

@testset "Parallel Evaluation" begin
    import CMAEvolutionStrategy: Optimizer, evaluate
    o1 = Optimizer(ones(3), 1)
    o2 = Optimizer(ones(3), 1, parallel_evaluation = true)
    o3 = Optimizer(ones(3), 1, multi_threading = true)
    f1(x) = sum(abs2, x)
    f2(x) = [sum(abs2, x[:, i]) for i in 1:size(x, 2)]
    input = [1 3; 2 4; 3 5]
    @test o2.p.parallel_evaluation == true
    @test o3.p.multi_threading == true
    @test evaluate(o1.p, f1, input) == evaluate(o2.p, f2, input) == evaluate(o3.p, f1, input)
end

@testset "Noisy" begin
    f(x) = sum(abs2, x)
    res1 = minimize(f, zeros(3), .5, maxiter = 50)
    res2 = minimize(f, zeros(3), .5,
                    noise_handling = CMAEvolutionStrategy.NoiseHandling(3),
                    maxiter = 50)
    @test fbest(res2) < 1e-3
    @test CMAEvolutionStrategy.noisefevals(res2.p.noise_handling) > 0
end

@testset "Logging" begin
    f_log = let xhist = Vector{Float64}[], fhist = Float64[]
        x -> begin
            f = sum(abs2, x) + randn() * .1
            push!(xhist, copy(x))
            push!(fhist, f)
            f
        end
    end
    res = minimize(f_log, ones(4), .25, verbosity = 0,
                   noise_handling = CMAEvolutionStrategy.NoiseHandling(4),
                   maxiter = 50, lower = ones(4), upper = 2 * ones(4))
    idx = argmin(f_log.fhist)
    @test fbest(res) == f_log.fhist[idx]
    @test xbest(res) == f_log.xhist[idx]
end

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy


text:
module GoldenSequences
#using Random
using FixedPointNumbers
using Roots
using Base.Iterators
using Base.MathConstants
export GoldenIntSequence, GoldenSequence

"""
	nextodd(x)

Next odd integer or next fixpoint number with odd representation.
"""
nextodd(i::Integer) = i >> 1 << 1 + one(i)
nextodd(z::T) where {T<:Normed} = reinterpret(T, nextodd(z.i))
prevodd(i::Integer) = i >> 1 << 1 - one(i)
prevodd(z::T) where {T<:Normed} = reinterpret(T, prevodd(z.i))

struct GoldenSequence{T}
	x0::T
	coeff::T
end
Iterators.eltype(::Type{GoldenSequence{T}}) where {T} = T
Iterators.IteratorEltype(::GoldenSequence) = Base.HasEltype()

Iterators.IteratorSize(::GoldenSequence) = Iterators.IsInfinite()

struct GoldenIntSequence{T}
	x0::T
	coeff::T
end
Iterators.eltype(::Type{GoldenIntSequence{T}}) where {T} = T
Iterators.IteratorSize(::GoldenIntSequence) = Iterators.IsInfinite()

f(x, d) = x^(d+1) - (x + one(x))

"""
Compute or look up solution to x^(d+1) = x + 1
"""
function Phi(_, d)
	if d <= length(phis)
		phis[d]
	else
		find_zero(x->f(x, d), (1.0, 2.0))
	end
end

function GoldenSequence(x0)
	phi = Phi(eltype(T), n)
	t = ntuple(n) do i
	       T(phi^(-i))
	end
	GoldenSequence{typeof(t)}(x0, false*x0 .+ t)
end
function GoldenSequence(x0::NTuple{N,T}) where {N, T}
	phi = Phi(T, N)
	t = ntuple(N) do i
	       T(phi^(-i))
	end
	GoldenSequence{typeof(x0)}(x0, t)
end
GoldenSequence(n::Int) = GoldenSequence(ntuple(n) do i; 0.0 end)
GoldenSequence(T, n) = GoldenSequence(ntuple(n) do i; zero(T) end)

function GoldenSequence(x0::T) where {T <: Number}
	GoldenSequence{T}(x0, convert(T, golden))
end

function GoldenSequence(::Type{T}, n) where {T<:Normed}
	phi = Phi(T, n)
	t = ntuple(n) do i
	       prevodd(T(phi^(-i)))
	end
	GoldenSequence{typeof(t)}(t)
end

function GoldenIntSequence(::Type{T}, n) where {T<:Unsigned}
	phi = Phi(T, n)
	t = ntuple(n) do i
	       (round(T, (typemax(T)+1)*phi^(-i)))
	end
	GoldenIntSequence{typeof(t)}(t)
end

function Base.getindex(g::GoldenSequence{<:Tuple}, ind)
    ntuple(length(g.coeff)) do i
         (g.x0[i] + g.coeff[i]*ind) % true
    end
end
function Base.getindex(g::GoldenSequence, ind)
     (g.x0 .+ g.coeff.*ind) % true
end

function Base.iterate(g::GoldenSequence, val=g.x0)
	 state = (g.coeff .+ val) .% true
	 val, state
end
function Base.iterate(g::GoldenIntSequence, val=g.x0)
	 state = g.coeff .+ val # intentional overflow
	 val, state
end

const phis = [
1.618033988749895 # MathConstants.golden
1.324717957244746 # plastic number
1.2207440846057596
1.1673039782614187
1.1347241384015196
1.1127756842787055
1.09698155779856
1.085070245491451
1.0757660660868373
1.0682971889208415
1.0621691678642553
1.0570505752212285
1.0527109201475584
1.0489849347570346
1.045751024156342
1.0429177323017866
1.0404149477818474
1.0381880194364501
1.0361937171306836
1.034397396133807
1.032770966441043
1.0312914124792476
1.0299396966705234
1.0286999355282511
1.0275587723930215
1.0265048941462673
1.0255286547645814
1.0246217791365844
1.0237771278614651
1.0229885088662907
1.022250525317837
1.0215584519243472
1.020908133630809
1.0202959021164726
1.0197185065485594
1.0191730558310543
1.0186569701821966
1.0181679403286723
1.0177038929544147
1.0172629613133726
1.0168434601276772
1.0164438640594171
1.0160627891762175
1.0156989769358962
1.015351280299597
1.0150186516505193
1.0147001322501537
1.0143948430084788
1.0141019763809709
1.0138207892351239
1.013550596553797
1.0132907658630492
1.0130407122890395
1.0127998941626677
1.0125678091024284
1.0123439905158556
1.0121280044682868
1.0119194468747255
1.0117179409765626
1.0115231350700071
1.0113347004574071
1.01115232959636
1.010975734424683
1.0108046448420556
1.0106388073314982
1.0104779837058908
1.0103219499664957
1.0101704952619777
1.0100234209377514
1.00988053966664
1.0097416746528447
1.009606658902116
1.0094753345517855
1.009347552255012
1.009223170614191
1.0091020556590062
1.0089840803650738
1.0088691242095444
1.008757072760392
1.00864781729645
1.0085412544555385
1.0084372859082966
1.008335818055544
1.0082367617472274
1.0081400320211664
1.0080455478599986
1.0079532319648552
1.0078630105444437
1.0077748131183248
1.0076885723332831
1.0076042237917846
1.007521705891604
1.007440959675782
1.007361928692145
1.0072845588616808
1.0072087983551326
1.0071345974772095
1.0070619085578796
1.006990685850238
1.0069208854344955
]

end


text:
# audio.jl

@testset "audio.jl" begin
    # Example audio
    audio_one_channel = SampleBuf(rand(1000), 10)
    audio_two_channel = SampleBuf(rand(1000,2), 10)
    audio_multi_channel = SampleBuf(rand(1000,2), 10)

    @testset "mono" begin
        @test size(mono(audio_two_channel)) == (1000,)
        @test size(mono(audio_multi_channel)) == (1000,)

    end

    @testset "resample.jl" begin
        @test size(resample(audio_one_channel, 5))[1] == 500
        @test size(resample(audio_two_channel, 5))[1] == 500
        @test size(resample(audio_multi_channel, 5))[1] == 500
    end

    @testset "duration" begin
        @test duration(audio_one_channel) == 100
        @test duration(audio_two_channel) == 100
        @test duration(audio_multi_channel) == 100
    end

    # @testset "pitchshift" begin
    #     @test duration(pitchshift(audio_one_channel)) == duration(audio_one_channel)
    #     @test duration(pitchshift(audio_two_channel)) == duration(audio_two_channel)
    #     @test duration(pitchshift(audio_multi_channel)) == duration(audio_multi_channel)
    # end

    # @testset "speedup" begin
    #     @test duration(speedup(audio_one_channel, 2)) == duration(audio_one_channel) // 2
    #     @test duration(speedup(audio_two_channel, 2)) == duration(audio_two_channel) // 2
    #     @test duration(speedup(audio_multi_channel, 2)) == duration(audio_multi_channel) // 2
    # end

    # @testset "slowdown" begin
    #     @test duration(slowdown(audio_one_channel, 2)) == duration(audio_one_channel) * 2
    #     @test duration(slowdown(audio_two_channel, 2)) == duration(audio_two_channel) * 2
    #     @test duration(slowdown(audio_multi_channel, 2)) == duration(audio_multi_channel) * 2
    # end

end


text:
module MsgPackRpcClientSession

include("socks.jl")
include("const.jl")

type Session
  socks   :: MsgPackRpcClientSocks.Socks
  ptr     :: Int
  next_id :: Int

  create      :: Function
  destroy     :: Function
  set_socks   :: Function
  set_ptr     :: Function
  set_next_id :: Function
  get_sock    :: Function
  get_socks   :: Function
  get_ptr     :: Function
  get_next_id :: Function
  rotate      :: Function

  function Session(socks = MsgPackRpcClientSocks.Socks({}), ptr = 0, next_id = 1)
    this         = new()
    this.socks   = socks
    this.ptr     = ptr
    this.next_id = next_id
    this.create      = function() create(socks, next_id) end
    this.destroy     = function() destroy(this) end
    this.set_socks   = function() set_socks(this, socks) end
    this.set_ptr     = function() set_ptr(this, ptr) end
    this.set_next_id = function() set_next_id(this, next_id) end
    this.get_sock    = function() get_sock(this) end
    this.get_socks   = function() get_socks(this) end
    this.get_ptr     = function() get_ptr(this) end
    this.get_next_id = function() get_next_id(this) end
    this.rotate      = function() rotate(this) end
    this
  end
end

function create(socks::MsgPackRpcClientSocks.Socks = MsgPackRpcClientSocks.Socks({}),
                ptr = 0,
                next_id = 1)
  Session(socks, ptr, next_id)
end

function destroy(self::Session)
  try
    if length(self.socks.pool) > 0
      MsgPackRpcClientSocks.destroy(self.socks)
    end
    self.ptr = 1
    self.next_id = 1
  catch
    return false
  end
  true
end

function set_socks(self::Session, socks::Any)
  self.socks = socks
  self
end

function set_ptr(self::Session, ptr::Int)
  self.ptr = ptr
  self
end

function set_next_id(self::Session, next_id::Int)
  self.next_id = next_id
  self
end

function get_sock(self::Session)
  self.socks.pool[self.ptr]
end

function get_socks(self::Session)
  self.socks
end

function get_ptr(self::Session)
  self.ptr
end

function get_next_id(self::Session)
  self.next_id
end

function create_sock(self::Session, host::String = DEFAULT_HOST, port::Int = DEFAULT_PORT_NUMBER)
  sock = connect(host, port)
  MsgPackRpcClientSocks.enqueue!(self.socks, sock)
  if self.ptr == 0
    self.ptr = 1
  end
  self
end

function rotate(self::Session)
  if self.ptr == 0
    # TODO: raise
  end
  self.ptr = self.ptr + 1 > length(self.socks.pool) ? 1 : self.ptr + 1
  self.socks.pool[self.ptr]
end

end # module MagPackRpcClientSession

# module session
#
# export Session
#
# type Loop end
# type Tcp end
# type Method end
# type Reason end
# type Error end
# type Result end
#
# type Session
#   address::String
#   timeout::Int
# #  loop::Loop
# #  builder::Tcp
#   reconnect_limit::Int
#   pack_encoding::String
#   unpack_encoding::String
#
#   call::Function
#   call_async::Function
#   send_request::Function
#   notify::Function
#   close::Function
#   on_connect_failed::Function
#   on_response::Function
#   on_timeout::Function
#   step_timeout::Function
#   NoSyncIDGenerator::Function
#
#   request_table::Array
#
#   function Session(address, timeout, loop, builder, reconnect_limit, pack_encoding, unpack_encoding)
#     this                   = new()
#     this.address           = address
#     this.timeout           = timeout
# #    this.loop              = loop
# #    this.builder           = builder
#     this.reconnect_limit   = reconnect_limit
#     this.pack_encoding     = pack_encoding
#     this.unpack_encoding   = unpack_encoding
#
#     this.call              = function() call(this, method, args) end
#     this.call_async        = function() call_async(this, method, args) end
#     this.send_request      = function() send_request(this, method, args) end
#     this.notify            = function() notify(this, method, args) end
#     this.close             = function() close(this) end
#     this.on_connect_failed = function() on_connect_failed(self, reason) end
#     this.on_response       = function() on_response(self, msgid, error, result) end
#     this.on_timeout        = function() on_timeout(self, msgid) end
#     this.step_timeout      = function() step_timeout(self) end
#     this.generator         = function() NoSyncIDGenerator() end
#
#     this.request_table     = []
#
#     this
#   end
#
#   function call(self::Session, method::Method, args::Array)
#     # return self.send_request(method, args).get()
#   end
#
#   function call_async(self::Session, method::Method, args::Array)
#     # return self.send_request(method, args)
#   end
#
#   function send_request(self::Session, method::Method, args::Array)
#     # msgid = next(self.generator)
#     # future = Future(self.loop, self.timeout)
#     # self.request_table[msgid] = future
#     # self.transport.send_message([message.REQUEST, msgid, method, args])
#     # return future
#   end
#
#   function notify(self::Session, method::Method, args::Array)
#     # function callback
#     #   self.loop.stop()
#     # end
#     # self.transport.send_message([message.NOTIFY, method, args], callback = callback)
#     # self.loop.start()
#   end
#
#   function close(self::Session)
#     # if self.transport
#     #   self.transport.close()
#     # end
#     # self.transport = nothing
#     # self.request_table = {}
#   end
#
#   function on_connect_failed(self::Session, reason::Reason)
#   end
#
#   function on_response(self::Session, msgid::Int, error::Error, result::Result)
#     # if not msgid in self._request_table
#     #     return
#     # end
#     # future = self.request_table.pop(msgid)
#
#     # if error is not nothing
#     #   future.set_error(error)
#     # else
#     #   future.set_result(result)
#     # end
#     # self.loop.stop()
#   end
#
#   function on_timeout(self::Session, msgid::Int)
#     # future = self.request_table.pop(msgid)
#     # future.set_error("Request timed out")
#   end
#
#   function step_timeout(self::Session)
#   end
#
#   function NoSyncIDGenerator()
#     counter = 0
#     while true
#       counter += 1
#       if counter > (1 << 30)
#         counter = 0
#         break
#       end
#     end
#     return counter
#   end
#
# end
#
# end # module session

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
module CombinedUncertainDiffEq

# Imported functions
import ReachabilityAnalysis
const RA = ReachabilityAnalysis
const ReachInterval = RA.Interval

using ComponentArrays: ComponentArrays, ComponentArray, getaxes, getdata
using ConcreteStructs: @concrete
using DifferentialEquations: DifferentialEquations, ODEProblem, EnsembleProblem, solve, remake, isinplace
using DiffEqUncertainty: DiffEqUncertainty, Koopman, MonteCarlo, expectation, _rand
using Distributions: Distribution, Uniform
using IntervalArithmetic: Interval, ±, (..), inf, sup
using ReachabilityAnalysis: HalfSpace, HPolyhedron, Hyperrectangle,
        ×, cartesian_product, concretize
using Setfield: @set!
using Statistics: mean
using Symbolics: @variables, Num
using UnPack: @unpack

# Exported functions
export ±, (..), Interval, Uniform

include("soft_functions.jl")

include("sim_utils.jl")
export Saturated, SoftSaturated, WithControls, reachable_constraint

include("conversions.jl")

include("solve.jl")
export mc_solve, mc_expectation, koopman_expectation, reachability_problem, reachability_solve

end


text:
###################################################################################################
# CONSTANTS
###################################################################################################
const LOG2 = log(2)                         # Ln(2)
const XCAL1 = [-1,1]                        # 𝒳 with N=1
const XCAL2 = [[-1,-1],[-1,1],[1,-1],[1,1]] # 𝒳 with N=2
###################################################################################################
# STRUCTS
###################################################################################################
mutable struct IntBitVec <: AbstractVector{Bool}
    data::UInt64
end
Base.size(::IntBitVec) = (64,)
@inline function Base.getindex(v::IntBitVec,i::Int)
   @boundscheck checkbounds(v, i)
   v.data & (1 << (i-1)) != 0
end
###################################################################################################
# USED FUNCTIONS
###################################################################################################
"""
`comp_ex([N1,...,NK],[α1,...,αK],β)`

Function that computes the mean methylation vector E[X] assuming the Ising model for a methylation
state vector of size `[N1,...,NK]` and parameters `[α1,...,αK]` and `β`.

# Examples
```julia-repl
julia> CpelAsm.comp_ex([4],[0.0],0.0)
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
```
"""
function comp_ex(n::Vector{Int64},a::Vector{Float64},b::Float64)::Vector{Float64}

    # Loop over all positions
    x = zeros(Int64,sum(n))
    y = zeros(Float64,sum(n))
    @inbounds for j in 1:sum(n)
        x[j] = 1
        y[j] = comp_lkhd(x,n,a,b)
        x[j] = 0
    end

    # Return
    return 2*y.-1

end # end comp_ex
"""
`comp_exx([N1,...,NK],[α1,...,αK],β;r=1)`

Function that computes mean methylation vector E[X_{i}X_{i+r}}] assuming the Ising model for a
methylation state vector of size `[N1,...,NK]` and parameters `[α1,...,αK]` and `β`.

# Examples
```julia-repl
julia> CpelAsm.comp_exx([4],[0.0],0.0)
3-element Array{Float64,1}:
 0.0
 2.220446049250313e-16
 0.0
```
"""
function comp_exx(n::Vector{Int64},a::Vector{Float64},b::Float64;r::Int64=1)::Vector{Float64}

    # Loop over all positions
    x = zeros(Int64,sum(n))
    y = zeros(Float64,sum(n)-r)
    @inbounds for j in 1:(sum(n)-r)
        x[j] = x[j+r] = 1
        y[j] = comp_lkhd(x,n,a,b) + comp_lkhd(-x,n,a,b)
        x[j] = x[j+r] = 0
    end

    # Return
    return 2*y.-1

end # end comp_exx
"""
    `comp_mml(Z,EX)`

Function that computes mean methylation level (MML) given at the CpG sites marked by binary vector
Z and given first order moment E[X] vector. This function is only used when not all CpG sites are
homozygous since it can compute MML over a subset of CpG sites.

# Examples
```julia-repl
julia> ex=CpelAsm.comp_ex([4],[0.0],0.0);
julia> CpelAsm.comp_mml(trues(4),ex)
0.5
```
"""
function comp_mml(z::BitArray{1},ex::Vector{Float64})::Float64

    # Return
    return abs(round(0.5/length(ex[z])*sum(ex[z])+0.5;digits=8))

end # end comp_mml
"""
    `comp_mml_∇([N1,...,NK],∇logZ)`

Function that computes mean methylation level (MML) over the entire haplotype. This function takes
advantage of the fact that E[S(X)]=Adot(θ), and it's much faster than `comp_mml`. This function is
only used when all CpG sites are homozygous.

# Examples
```julia-repl
julia> n=[1,1,1]; θ=[1.0,-1.0,1.0,0.0];
julia> ∇logZ = CpelAsm.get_grad_logZ(n,θ);
julia> CpelAsm.comp_mml_∇(n,∇logZ)
0.62693236
```
"""
function comp_mml_∇(n::Vector{Int64},∇logZ::Vector{Float64})::Float64

    # Return
    return abs(round(0.5*(1.0+1.0/sum(n)*sum(∇logZ[1:length(n)]));digits=8))

end # end comp_mml_∇
"""
    `comp_exlng(Z,[N1,...,NK],[α1,...,αK],β)`

Function that returns the third summand in the normalized methylation entropy (NME) defined
for the homozygous CpG sites, with positions determined by binary vector Z, assuming an
allele-specific vector with [N1,...,NK] CpG sites, with parameters [α1,...,αK] and β.

# Examples
```julia-repl
julia> n=[10]; z=trues(10); z[5]=false; α=[0.0]; β=0.0;
julia> CpelAsm.comp_exlng(z,n,α,β)
0.6931471805599452
```
"""
function comp_exlng(z::BitArray{1},n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Find changes to/from 0 in x vector
    hetst = findall(isequal(-1),abs.(z[2:end]) - abs.(z[1:(end-1)])) .+ 1
    hetend = findall(isequal(1),abs.(z[2:end]) - abs.(z[1:(end-1)])) .+ 1

    # Determine whether it starts/finishes with true/false
    z[1]==false && pushfirst!(hetst,1)
    z[end]==false && push!(hetend,sum(n)+1)

    # Find subregion label for each CpG site
    subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)

    # Loop over heterozygous stretches
    y = 0.0
    @inbounds for i=1:length(hetst)

        # Choose 𝒳 depending on boundary or not
        bound = hetst[i]==1 || hetend[i]==sum(n)+1 ? true : false
        xcal = bound ? XCAL1 : XCAL2

        # Figure out b (block IDs of heterozygous) and r (α indices)
        bid = subid[hetst[i]:(hetend[i]-1)]
        nhet = [count(x->x==id,bid) for id in unique(bid)]

        # Compute expectation
        for x in xcal

            # Find αp1 and αp2
            xaug = zeros(Int64,sum(n))
            if bound
                # If heterozygous CpG on boundary of X
                hetst[i]==1 ? xaug[hetend[i]]+=x : xaug[hetst[i]-1]+=x
                ap1 = hetst[i]==1 ? a[1] : 2.0*x*b+a[subid[hetst[i]]]
                ap2 = hetend[i]==sum(n)+1 ? a[end] : 2.0*x*b+a[subid[hetend[i]-1]]
            else
                # If heterozygous CpG not on boundary of X
                xaug[[hetst[i]-1,hetend[i]]]+=x
                ap1 = 2.0*x[1]*b+a[subid[hetst[i]]]
                ap2 = 2.0*x[2]*b+a[subid[hetend[i]-1]]
            end

            # Add contribution to total sum of expectations
            y += log(comp_g(nhet,a[unique(bid)],b,ap1,ap2))*comp_lkhd(xaug,n,a,b)

        end

    end

    # Return sum of all expectations
    return y

end # end comp_exlng
"""
    `comp_nme(Z,[N1,...,NK],[α1,...,αK],β,EX,EXX)`

Function that computes normalized methylation entropy (NME) over the CpG sites determined by binary
vector Z. This is done by assuming an Ising model for the allele-specific methylation state vector
of size `[N1,...,NK]`, parameters `[α1,...,αK]` and `β`. This function is used when not all CpG
sites are homozygous since it can compute NME over a subset of CpG sites.

# Examples
```julia-repl
julia> n=[10]; z=trues(sum(n)); z[5]=false; α=[0.0]; β=0.0;
julia> CpelAsm.comp_nme(z,n,α,β,CpelAsm.comp_ex(n,α,β),CpelAsm.comp_exx(n,α,β))
1.0
```
"""
function comp_nme(z::BitArray{1},n::Vector{Int64},a::Vector{Float64},b::Float64,ex::Vector{Float64},
                  exx::Vector{Float64})::Float64

    # Define output quantity
    h = 0.0

    # Vector of α for homozygous
    avec = vcat([a[i]*ones(Float64,n[i]) for i in 1:length(n)]...)[z]

    # 1st term (log partition function)
    h += log(comp_Z(n,a,b))

    # 2nd term (average potential energy over homozygous CpG sites)
    h -= avec'*ex[z]+b*sum(exx[z[1:(end-1)].*z[2:end]])

    # 3rd term if necessary
    h -= all(z) ? 0.0 : comp_exlng(z,n,a,b)

    # Return nme
    return abs(round(h/(sum(z)*LOG2);digits=8))

end # end comp_nme
"""
    `comp_nme_∇([N1,...,NK],θ,∇logZ)`

Function that computes normalized methylation entropy (NME) using the ∇logZ information. This
function can only be used in `comp_tnull` since it does not allow computation on a subset of CpG
sites. This function is only used when all CpG sites are homozygous.

# Examples
```julia-repl
julia> n=[4]; α=[0.0]; β=0.0; θ=vcat([α,β]...); ∇logZ=CpelAsm.get_grad_logZ(n,θ);
julia> CpelAsm.comp_nme_∇(n,θ,∇logZ)
1.0
```
"""
function comp_nme_∇(n::Vector{Int64},θ::Vector{Float64},∇logZ::Vector{Float64})::Float64

    # Return
    return abs(round(1.0/(sum(n)*LOG2)*(log(comp_Z(n,θ[1:(end-1)],θ[end]))-θ'*∇logZ);digits=8))

end # end comp_nme_∇
"""
    `gen_x_mc(N,α,β)`

Function that generates a methylation vector from an Ising model with `[N1,...,NK]`, and parameters
`[α1,...,αK]` and `β`.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> CpelAsm.gen_x_mc([5],[0.0],0.0)
5-element Array{Int64,1}:
  1
  1
  1
 -1
  1
```
"""
function gen_x_mc(n::Vector{Int64},a::Vector{Float64},b::Float64)::Vector{Int64}

    # Find subregion label for each CpG site
    subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)

    # Sample first CpG site
    x = Vector{Int64}()
    p = comp_lkhd(pushfirst!(zeros(Int64,sum(n)-1),1),n,a,b)
    push!(x,sample(XCAL1,Weights([1-p,p])))
    sum(n)>1 || return x

    # Sequentially sample from Markov chain
    @inbounds for i=2:sum(n)

        # Add boundary function g() if necessary
        expaux = exp(a[subid[i]]+b*x[i-1])
        if i < sum(n)
            ap1p = 2.0*b + a[subid[i]]
            ap1m = -2.0*b + a[subid[i]]
            n_miss = [count(x->x==sr,subid[(i+1):sum(n)]) for sr in unique(subid[(i+1):sum(n)])]
            gp = comp_g(n_miss,a[unique(subid[(i+1):sum(n)])],b,ap1p,a[end])
            gm = comp_g(n_miss,a[unique(subid[(i+1):sum(n)])],b,ap1m,a[end])
            p = gp*expaux/(gp*expaux+gm/expaux)
        else
            p = expaux/(expaux+1.0/expaux)
        end

        # Add i-th CpG site
        push!(x,sample(XCAL1,Weights([1-p,p])))

    end

    # Return methylation vector
    return x

end # end gen_x_mc
"""
    `comp_nme_mix_mc(Z1,Z2,N1,N2,theta1,theta2)`

Function that estimates the entropy for the Ising mixture model using Monte Carlo. This is done
assuming an Ising model for the allele-specific methylation state vector of size `[N1,...,NK]`,
parameters `[α1,...,αK]` and `β`, for each allele. The homozygous part of each model is determined
by binary vector Z (i.e., via Hadamard product Z*X, where * is the Hadamard product).

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> α1=[0.0]; β1=0.0; α2=[1.0]; β2=1.0;
julia> n1=[10]; n2=[10]; z1=trues(sum(n1)); z2=trues(sum(n2));
julia> CpelAsm.comp_nme_mix_mc(z1,z2,n1,n2,vcat(α1,β1),vcat(α2,β2))
0.6278665195419888
```
"""
function comp_nme_mix_mc(z1::BitArray{1},z2::BitArray{1},n1::Vector{Int64},n2::Vector{Int64},
                         t1::Vector{Float64},t2::Vector{Float64};L::Int64=1000)::Float64

    # For loop to sample from mixture
    h = 0.0
    @inbounds for i=1:L

        # Sample allele
        y = sample(1:2,Weights([0.5,0.5]))

        # Sample x given allele
        x = y==1 ? gen_x_mc(n1,t1[1:(end-1)],t1[end]) : gen_x_mc(n2,t2[1:(end-1)],t2[end])
        x0 = y==1 ? x[z1] : x[z2]

        # Add contribution
        x1 = zeros(Int64,sum(n1))
        x2 = zeros(Int64,sum(n2))
        x1[z1] = x0
        x2[z2] = x0
        h += log(comp_lkhd(x1,n1,t1[1:(end-1)],t1[end])+comp_lkhd(x2,n2,t2[1:(end-1)],t2[end]))

    end

    # Return
    return min(1.0,max(0.0,1.0/sum(z1)*(1.0-h/(L*LOG2))))

end # end comp_nme_mix_mc
"""
    `comp_nme_mix_exact(Z1,Z2,N1,N2,theta1,theta2)`

Function that exactly computes the NME for the Ising mixture model. This is done by assuming an
Ising model for the allele-specific methylation state vector of size `[N1,...,NK]`, parameters
`[α1,...,αK]` and `β`. The homozygous part of the allele-specific vector is determined by binary
vector Z (i.e., via Hadamard product Z*X, where * is the Hadamard product). This is done by
traversing 𝒳h.

# Examples
```julia-repl
julia> n1=[10]; n2=[10];
julia> α1=[0.0]; β1=0.0; α2=[1.0]; β2=1.0;
julia> z1=trues(sum(n1)); z2=trues(sum(n2));
julia> CpelAsm.comp_nme_mix_exact(z1,z2,n1,n2,vcat(α1,β1),vcat(α2,β2))
0.6239489595165986
```
"""
function comp_nme_mix_exact(z1::BitArray{1},z2::BitArray{1},n1::Vector{Int64},n2::Vector{Int64},
                            t1::Vector{Float64},t2::Vector{Float64})::Float64

    # Loop over 𝒳h
    h = 0.0
    n = sum(z1)
    w1 = zeros(Int64,sum(n1))
    w2 = zeros(Int64,sum(n2))
    @inbounds for i=1:(2^n)
        z = IntBitVec(i)[1:n]
        x = -ones(Int64,n)
        x[z] .= 1
        w1[z1] =+ x
        w2[z2] =+ x
        lkhd1 = comp_lkhd(w1,n1,t1[1:(end-1)],t1[end])
        lkhd2 = comp_lkhd(w2,n2,t2[1:(end-1)],t2[end])
        h += (lkhd1+lkhd2) * log(lkhd1+lkhd2)
        w1[z1] =- x
        w2[z2] =- x
    end

    # Return
    return min(1.0,max(0.0,1.0/n*(1.0-0.5*h/LOG2)))

end # end comp_nme_mix_exact
"""
    `comp_uc(Z1,Z2,N1,N2,theta1,theta2,h1,h2)`

Function that exactly computes UC over the homozygous CpG sites. This is done by assuming an Ising
model for the allele-specific methylation state vector of size `[N1,...,NK]`, parameters
`[α1,...,αK]` and `β`. The homozygous part of the allele-specific vector is determined by binary
vector Z (i.e., via Hadamard product Z*X, where * is the Hadamard product).

# Examples
```julia-repl
julia> n1=[20]; n2=[20];
julia> z1=trues(sum(n1)); z2=trues(sum(n2));
julia> α1=[2.5]; α2=[2.5]; β1=0.0; β2=0.0;
julia> h1=CpelAsm.comp_nme(z1,n1,α1,β1,CpelAsm.comp_ex(n1,α1,β1),CpelAsm.comp_exx(n1,α1,β1));
julia> h2=CpelAsm.comp_nme(z2,n2,α2,β2,CpelAsm.comp_ex(n2,α2,β2),CpelAsm.comp_exx(n2,α2,β2));
julia> CpelAsm.comp_uc(z1,z2,n1,n2,vcat(α1,β1),vcat(α2,β2),h1,h2)
0.63304398
```
"""
function comp_uc(z1::BitArray{1},z2::BitArray{1},n1::Vector{Int64},n2::Vector{Int64},
                  t1::Vector{Float64},t2::Vector{Float64},h1::Float64,h2::Float64)::Float64

    # Compute h(x)
    h = sum(z1)<17 ? comp_nme_mix_exact(z1,z2,n1,n2,t1,t2) : comp_nme_mix_mc(z1,z2,n1,n2,t1,t2)
        # print_log("h: $(h)")

    # If H(X)≈0 then UC=0
    uc = isapprox(h,0.0,atol=1e-5) ? 0.0 : min(1.0,max(0.0,1.0-0.5*(h1+h2)/h))

    # Return
    return round(uc;digits=8)

end # end comp_uc
###################################################################################################
# DEV FUNCTIONS
###################################################################################################
"""
    `comp_mmlv([N1,...,NK],∇logZ)`

Function that computes the mean methylation level vector (MMLV) for each subregion. Reporting
output for subregion might be more informative than just an overall mean. This function takes
advantage of the fact that E[S(X)]=Adot(θ), which is much faster than computing first moment
and correlations as currently done.

# Examples
```julia-repl
julia> n=[1,1,1]; θ=[1.0,-1.0,1.0,0.0];
julia> ∇logZ = CpelAsm.get_grad_logZ(n,θ);
julia> CpelAsm.comp_mmlv(n,∇logZ)
3-element Array{Float64,1}:
 0.88079708
 0.11920292
 0.88079708
```
"""
function comp_mmlv(n::Vector{Int64},∇logZ::Vector{Float64})::Vector{Float64}

    # Return
    return abs.(round.(0.5.*(1.0./n.*∇logZ[1:length(n)].+1.0);digits=8))

end # end comp_mmlv
"""
    `comp_corr(EX,EXX)`

Function that returns the correlation vector between consecutive CpG sites given the vector of
first order moments E[X] and second order moments E[XX].

# Examples
```julia-repl
julia> CpelAsm.comp_corr(CpelAsm.comp_ex([4],[0.0],0.0),CpelAsm.comp_exx([4],[0.0],0.0))
3-element Array{Float64,1}:
 0.0
 0.0
 0.0
```
"""
function comp_corr(ex::Vector{Float64},exx::Vector{Float64})::Vector{Float64}

    # Initialize vector of 1's
    ov = ones(Float64,length(ex)-1)

    # Return ρ
    return round.((exx.-ex[1:(end-1)].*ex[2:end]) ./ sqrt.((ov.-ex[1:(end-1)].^2) .*
    (ov.-ex[2:end].^2));digits=8)

end # end comp_corr
"""
    `comp_cov([N1,...,NK],[α1,...,αK],β,EX,EXX)`

Function that returns the covariance matrix of a methylation vector given the `[N1,...,NK]`,
`[α1,...,αK]`, `β`, E[X], and E[XX].

# Examples
```julia-repl
julia> n=[4]; α=[1.0]; β=1.0;
julia> ex = CpelAsm.comp_ex(n,α,β);
julia> exx = CpelAsm.comp_exx(n,α,β);
julia> CpelAsm.comp_cov(n,α,β,ex,exx)
4×4 Array{Float64,2}:
 0.419974     0.0          4.44089e-16  3.33067e-16
 0.0          0.419974     1.11022e-16  4.44089e-16
 4.44089e-16  1.11022e-16  0.419974     0.0
 3.33067e-16  4.44089e-16  0.0          0.419974
```
"""
function comp_cov(n::Vector{Int64},a::Vector{Float64},b::Float64,ex::Vector{Float64},
                  exx::Vector{Float64})::Array{Float64,2}

    # Initialize matrix
    ntot = sum(n)
    cov = zeros(Float64,ntot,ntot)

    # Add first subdiagonal E[X_{i}X_{i+1}]
    cov[2:(ntot+1):((ntot-1)*ntot)] = exx

    # Add covariance terms E[X_{i}X_{i+k-1}]
    @inbounds for k=3:ntot
        cov[k:(ntot+1):((ntot-k+1)*ntot)] = comp_exx(n,a,b;r=k-1)
    end

    # Symmetrize
    cov += transpose(cov)

    # Substract E[X_i]E[X_j]
    cov -= ex * ex'

    # Add 1's to diagonal of matrix
    cov[1:(ntot+1):ntot^2] += ones(Float64,ntot)

    # Return Σ
    return cov

end # end comp_cov
"""
    `comp_evec(Σ)`

Function that returns the eigenvector associated with the largest eigenvalue of the covariance
matrix Σ. The components of the vector with the largest magnitude highlight the CpG sites in the
region with the largest amount of variability.

# Examples
```julia-repl
julia> n=[4]; α=[1.0]; β=1.0;
julia> ex = CpelAsm.comp_ex(n,α,β);
julia> exx = CpelAsm.comp_exx(n,α,β);
julia> Σ = CpelAsm.comp_cov(n,α,β,ex,exx);
julia> CpelAsm.comp_evec(Σ)
4-element Array{Float64,1}:
 0.6867257391267225
 0.16854601514378278
 0.16854601514379314
 0.6867257391267322
```
"""
function comp_evec(cov::Array{Float64,2})::Vector{Float64}

    # Return evec
    return abs.(eigvecs(cov)[:,size(cov)[1]])

end # end comp_evec
###################################################################################################
# VERIFICATION FUNCTIONS (NOT USED)
###################################################################################################
"""
    `comp_nme_xcal(Z,[N1,...,NK],[α1,...,αK],β,EX,EXX)`

Function that computes normalized methylation entropy (NME) over the homozygous CpG sites. This is
done by assuming an Ising model for the allele-specific methylation state vector of size
`[N1,...,NK]`, parameters `[α1,...,αK]` and `β`. The homozygous part of the allele-specific vector
is determined by binary vector Z (i.e., via Hadamard product Z*X, where * is the Hadamard product).

# Examples
```julia-repl
julia> n=[10]; z=trues(sum(n)); α=[0.0]; β=0.0;
julia> CpelAsm.comp_nme_xcal(z,n,α,β)
1.0
```
"""
function comp_nme_xcal(z::BitArray{1},n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Loop over 𝒳h
    h = 0.0
    w = zeros(Int64,sum(n))
    hom_ind = findall(y->y==true,z)
    xcal = generate_xcal(length(hom_ind))
    @inbounds for x in xcal
        w[hom_ind] =+ x
        lkhd = comp_lkhd(w,n,a,b)
        h += lkhd * log(lkhd)
        w[hom_ind] =- x
    end

    # Return
    return round(-h/(sum(z)*LOG2);digits=8)

end # end comp_nme_xcal
"""
    `comp_uc_xcal(Z1,Z2,N1,N2,theta1,theta2)`

Function that exactly computes the UC over the homozygous CpG sites. This is done by assuming an
Ising model for the allele-specific methylation state vector of size `[N1,...,NK]`, parameters
`[α1,...,αK]` and `β`. The homozygous part of the allele-specific vector is determined by binary
vector Z (i.e., via Hadamard product Z*X, where * is the Hadamard product).

# Examples
```julia-repl
julia> n1=[10]; n2=[10];
julia> z1=trues(sum(n1)); z2=trues(sum(n2));
julia> θ1=[2.0,0.0]; θ2=[-2.0,0.0];
julia> CpelAsm.comp_uc_xcal(z1,z2,n1,n2,θ1,θ2)
0.43482166257188215
```
"""
function comp_uc_xcal(z1::BitArray{1},z2::BitArray{1},n1::Vector{Int64},n2::Vector{Int64},
                      theta1::Vector{Float64},theta2::Vector{Float64})::Float64

    # Loop over 𝒳h
    num = 0.0
    den = 0.0
    n = sum(z1)
    w1 = zeros(Int64,sum(n1))
    w2 = zeros(Int64,sum(n2))
    hom_ind_1 = findall(y->y==true,z1)
    hom_ind_2 = findall(y->y==true,z2)
    @inbounds for i=1:(2^n)
        z = IntBitVec(i)[1:n]
        x = -ones(Int64,n)
        x[z] .= 1
        w1[hom_ind_1] =+ x
        w2[hom_ind_2] =+ x
        lkhd1 = comp_lkhd(w1,n1,theta1[1:(end-1)],theta1[end])
        lkhd2 = comp_lkhd(w2,n2,theta2[1:(end-1)],theta2[end])
        lkhd0 = 0.5*(lkhd1+lkhd2)
        den += lkhd0 * log(lkhd0)
        num += lkhd1 * log(lkhd1) + lkhd2 * log(lkhd2)
        w1[hom_ind_1] =- x
        w2[hom_ind_2] =- x
    end

    # Return
    return 1-0.5*num/den

end # end comp_uc_xcal
###################################################################################################
# COMPETING MODELS
###################################################################################################
"""
    `comp_lkhd_bin(XOBS,θ)`

Function that returns likelihood vector of XOBS, potentially partial, given parameter vector θ for
the non-parametric independent model (binomial).

"""
function comp_lkhd_bin(xobs::Array{Vector{Int64},1},θ::Vector{Float64})::Vector{Float64}

    # Return evec
    return [prod(θ.^(x.==1) .* (ones(length(xobs[1]))-θ).^(x.==-1)) for x in xobs]

end # end comp_lkhd_bin
"""
    `comp_lkhd_mult(XOBS,θ)`

Function that returns likelihood vector of XOBS, assuming full data, given parameter vector θ for
the non-parametric dependent model (multinomial).

"""
function comp_lkhd_mult(xobs::Array{Vector{Int64},1},θ::Vector{Float64})::Vector{Float64}

    # Generate 𝒳
    xcal = generate_xcal(length(xobs[1]))

    # Return evec
    return θ[[findfirst(y->y==x,xcal) for x in xobs]]

end # end comp_lkhd_mult


text:
function prepString(arr::Array{Float64,2})
  s = string(arr)
  ss = split(s,'\n')
  ss[length(ss)] = ss[length(ss)][1:(end-1)]
  fs = ""
  for i in 1:length(ss)
    ss[i]=replace(ss[i]," ",",")
    if ss[i][1] == ',' || ss[i][1] == '['
      ss[i] = ss[i][2:end]
    end
    ss[i] = string(ss[i],";")
    fs = string(fs,ss[i])
  end
  fs[1:(end-1)]
end

function parseInit!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  if length(sp2) == 3
    println("parseInit! -- received initialization point $(sp2) NOT CURRENTLY USED")
    x0 = parse(Float64,sp2[1])
    y0 = parse(Float64,sp2[2])
    th0 = parse(Float64,sp2[2])
  end
  prevn = initFactorGraph!(slam.fg, ready=USEREADY)
  println("init done")
  nothing
end

posecount = 1
function parseOdo!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  global posecount
  # println("parseOdo -- ")
  posecount += 1
  n = string("x", posecount)
  DX = [parse(Float64,sp2[3]);parse(Float64,sp2[4]);parse(Float64,sp2[5])]
  cov = diagm([parse(Float64,sp2[6]);parse(Float64,sp2[9]);parse(Float64,sp2[11])])
  cov[1,2] = parse(Float64,sp2[7])
  cov[1,3] = parse(Float64,sp2[8])
  cov[2,3] = parse(Float64,sp2[10])
  cov[2,1] = parse(Float64,sp2[7])
  cov[3,1] = parse(Float64,sp2[8])
  cov[3,2] = parse(Float64,sp2[10])
  addOdoFG!(slam.fg, Symbol(n), DX, cov, ready=USEREADY)
  return n
end

function parseAddLandmBR!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  poseid = floor(Int,parse(Float64, sp2[1]))
  pose = string(getVert(slam.fg,poseid).label)
  lmid = floor(Int,parse(Float64, sp2[2]))
  zbr = [parse(Float64,sp2[3]);parse(Float64,sp2[4])]
  cov = diagm([parse(Float64,sp2[5]);parse(Float64,sp2[7])])
  cov[1,2] = parse(Float64,sp2[6])
  cov[2,1] = parse(Float64,sp2[6])
  lm = string("")
  if !haskey(slam.fg.v, lmid)
    slam.lndmidx += 1
    lm = string(string('l',slam.lndmidx))
    projNewLandm!(slam.fg, pose, lm, zbr, cov, ready=USEREADY)
  else
    lm = string(getVert(slam.fg,lmid).label)
    addBRFG!(slam.fg, pose, lm, zbr, cov, ready=USEREADY)
  end
  return lm
end

function parseLandmBRMM!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  poseid = floor(Int,parse(Float64, sp2[1]))
  pose = string(getVert(slam.fg, poseid).label)
  lm1id = floor(Int,parse(Float64, sp2[2]))
  w1 = parse(Float64, sp2[3])
  lm2id = floor(Int,parse(Float64, sp2[4]))
  w2 = parse(Float64, sp2[5])
  zbr = [parse(Float64,sp2[6]);parse(Float64,sp2[7])]
  cov = diagm([parse(Float64,sp2[8]);parse(Float64,sp2[10])])
  cov[1,2] = parse(Float64,sp2[9])
  cov[2,1] = parse(Float64,sp2[9])
  # TODO --do the tests here
  lm1 = string(getVert(slam.fg,lm1id).label)
  lm2 = string(getVert(slam.fg,lm2id).label)
  addMMBRFG!(slam.fg, pose, [lm1;lm2], zbr, cov, w=[w1;w2], ready=USEREADY)
  # function addMMBRFG!(fg::G, pose::AbstractString,
  #                   lm::Array{AbstractString,1}, br::Array{Float64,1},
  #                   cov::Array{Float64,2}; w=[0.5;0.5]) where G <: AbstractDFG
  return "$(lm1), $(lm2)"
end

function parseLandmBRAuto!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  poseid = floor(Int,parse(Float64, sp2[1]))
  pose = string(getVert(slam.fg,poseid).label)
  lmid = sp2[2] != "*" ? floor(Int,parse(Float64, sp2[2])) : -1
  zbr = [parse(Float64,sp2[3]);parse(Float64,sp2[4])]
  cov = diagm([parse(Float64,sp2[5]);parse(Float64,sp2[7])])
  cov[1,2] = parse(Float64,sp2[6])
  cov[2,1] = parse(Float64,sp2[6])
  lm = string("")

  vlm, flm, slam.lndmidx = addAutoLandmBR!(slam.fg, pose, lmid, zbr, cov, slam.lndmidx, ready=USEREADY)

  println("parseLandmBRAuto! -- added $(vlm.label)")

  return vlm.label
end

function parseLandmarkXY!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  pose = getVert(slam.fg, map(Int, sp2[1]) ).label
  lmid = map(Int, sp2[2])
  zxy = [map(Float64,sp[3]);map(Float64,sp[4])]
  cov = diagm([map(Float64,sp[5]);map(Float64,sp[7])])
  cov[1,2] = map(Float64,sp[6])
  cov[2,1] = map(Float64,sp[6])
  error("parseLandmarkXY! -- not finished implementing yet")
  nothing
end

function batchSolve!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  println("batchSolve -- wiping tree and solving")
  slam.tree = wipeBuildNewTree!(slam.fg)
  @time inferOverTree!(slam.fg, slam.tree, N=100)
  nothing
end

function parseGetparticles(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  val = getVal(slam.fg, string(sp2[1]))
  retstr = prepString(val)
  return retstr
end

function parseLS(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  xx,ll = ls(slam.fg)
  str = string("")
  [str = string(str, x, ",") for x in xx]
  if length(str) > 0
    str = string(str[1:(end-1)], ';')
  end
  [str = string(str, x, ",") for x in ll]
  @show "LS", str
  return length(str) > 0 ? str : nothing
end

function parseGetID(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  id = slam.fg.IDs[sp2[1]]
  return "$(id)"
end

function parseGetNextID(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  id = slam.fg.id + 1
  return "$(id)"
end

function parseReset(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  @warn "parseReset -- resetting factor graph"
  slam = SLAMWrapper(emptyFactorGraph(), Union{}, 0)
  nothing
end

function parseDotFG(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  writeGraphPdf(slam.fg)
  nothing
end

function parseSetReady!(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  setDBAllReady!(slam.fg)
  nothing
end

function parseMongoFileSave(slam::SLAMWrapper, sp2::Array{SubString{AbstractString},1})
  id = slam.fg.IDs[sp2[1]]
  cgid = slam.fg.cgIDs[id]
  cv = CloudGraphs.get_vertex(slam.fg.cg, cgid)
  # for f in sp2[2:end]
  	fid = open(sp2[2],"r")
  	imageData = readbytes(fid) # imageData::Vector{UInt8}
  	close(fid)
	  bdei = CloudGraphs.BigDataElement("keyframe-image", imageData)
    push!(cv.bigData.dataElements, bdei);

    fid = open(sp2[3],"r")
  	imageData = readbytes(fid) # imageData::Vector{UInt8}
  	close(fid)
	  bdei = CloudGraphs.BigDataElement("depthframe-image", imageData)
    push!(cv.bigData.dataElements, bdei);

  CloudGraphs.save_BigData!(slam.fg.cg, cv)
  println("Finished writing to Mongo.")

  nothing
end

function parseTCP!(slam::SLAMWrapper, line::AbstractString)
    sp = split(line,' ')
    f = +
    cmd = sp[1]
    goahead = true
    if cmd == "INIT"
      f = parseInit!
    elseif cmd == "ODOMETRY"
      f = parseOdo!
    elseif cmd == "LANDMBR"
      f = parseAddLandmBR!
    elseif cmd == "LANDMARK"
      f = parseLandmarkXY!
    elseif cmd == "LANDMBRMM"
      f = parseLandmBRMM!
    elseif cmd == "LANDMBRAUTO"
      f = parseLandmBRAuto!
    elseif cmd == "SETALLREADY"
      f = parseSetReady!
    elseif cmd == "BATCHSOLVE"  # do not call batch when using DB solver
      f = batchSolve!
    elseif cmd == "GETPARTICLES"
      f = parseGetparticles
    elseif cmd == "LS"
      f = parseLS
    elseif cmd == "GETID"
      f = parseGetID
    elseif cmd == "GETNEXTID"
      f = parseGetNextID
    elseif cmd == "DRAWFGPDF"
      f = parseDotFG
    elseif cmd == "MONGOFILE"
      f = parseMongoFileSave
    elseif cmd == "RESET"
      f = parseReset
    elseif cmd == "QUIT"
      println("parseTCP -- should quit now")
      return false, string("")
    else
      @warn "parseTCP -- I don't know what $(cmd) means"
      goahead = false
    end
    retstr = nothing
    goahead ? retstr = f(slam, sp[2:end]) : nothing
    return true, retstr
end


text:
# Tests for Monetary values
@testset "Monetary" begin

# Basic arithmetic
@testset "Arithmetic" begin
    @test 1USD + 2USD == 3USD
    @test 1USD + 2USD + 3USD == 6USD
    @test 1.5USD * 3 == 4.5USD
    @test 1USD * 2 + 2USD == 4USD
    @test 2 * 1USD * 3 == 6USD
    @test 1.11USD * 999 == 1108.89USD
    @test -1USD == -(1USD)
    @test +20USD == 20USD
    @test 1USD == majorunit(USD)
    @test 1USD == USD
    @test one(25USD) ≡ 1

    @testset "Division" begin
        @test 10USD / 1USD == 10.0
        @test 10USD ÷ 3USD == 3
        @test 10USD % 3USD == 1USD
        @test -10USD ÷ 3USD == -3
        @test -10USD % 3USD == -1USD
        @test divrem(10USD, 3USD) == (10USD ÷ 3USD, 10USD % 3USD)

        @test mod(10USD, 3USD) == 1USD
        @test fld(10USD, 3USD) == 3
        @test mod(-10USD, 3USD) == 2USD
        @test fld(-10USD, 3USD) == -4
        @test fldmod(-10USD, 3USD) == (-4, 2USD)

        if VERSION ≥ v"0.5-"
            @test mod1(10USD, 3USD) == 1USD
            @test fld1(10USD, 3USD) == 4
            @test mod1(-10USD, 3USD) == 2USD
            @test fld1(-10USD, 3USD) == -3
            @test mod1(9USD, 3USD) == 3USD
            @test fld1(9USD, 3USD) == 3
            @test fldmod1(9USD, 3USD) == (3, 3USD)
        end
    end

    @testset "Number-like" begin
        @test abs(-10USD) == 10USD
        @test abs(0USD) == 0USD
        @test abs(10USD) == 10USD
        @test sign(-10USD) == -1
        @test sign(0USD) == 0
        @test sign(10USD) == 1
    end
end

# Type safety
@testset "Type Safety" begin
    @test_throws MethodError 1USD + 1CAD
    @test_throws MethodError 100JPY - 0USD  # even when zero!
    @test_throws MethodError 5USD * 10CAD
    @test_throws MethodError 5USD * 10USD
    @test_throws MethodError 1USD + 1
    @test_throws MethodError 2 - 1JPY
    @test_throws MethodError 1USD / 1CAD
    @test_throws MethodError 10USD ÷ 5        # meaningless
    @test_throws MethodError 10USD % 5        # meaningless
end

@testset "Comparison" begin
    # Comparisons
    @testset "Monetary — Homogenous Comparisons" begin
        @test 1EUR < 2EUR
        @test 3JPY > 2JPY
        @test 3JPY ≥ 3JPY
        @test 9USD ≤ 9.01USD
        @test -1USD ≠ 0USD
        @test sort([0.5EUR, 0.7EUR, 0.3EUR]) == [0.3EUR, 0.5EUR, 0.7EUR]

        # immutable types should be **equivalent**
        @test 10USD ≡ 10USD
        @test zero(USD) ≡ 0USD
        @test -majorunit(USD) ≡ +(-USD)
    end

    # Type safety for comparisons
    @testset "Monetary — Heterogenous Comparisons" begin
        @test 1EUR ≠ 1USD
        @test 5USD ≠ 5
        @test 5USD ≠ 500

        @test_throws MethodError EUR > USD
        @test_throws MethodError GBP ≥ USD
        @test_throws MethodError JPY < USD
    end
end

# Big int monetary
BI_USD = Monetary(:USD, BigInt(100))
BI_USD2 = majorunit(Monetary{:USD, BigInt})
BI_USD3 = Monetary(:USD; storage=BigInt)
I128_USD = majorunit(Monetary{:USD, Int128})
@testset "Custom Representation" begin
    @test BigInt(2)^100 * BI_USD + 10BI_USD == (BigInt(2)^100 + 10) * BI_USD

    # test **equality** — note equivalence is untrue because BigInt
    @test BI_USD == BI_USD2 == BI_USD3

    # wrapping behaviour (strange but documented)
    # TODO: this test is wrong
    # @test typemin(Int128) * I128_USD ≡ typemax(Int128) * I128_USD + I128_USD

    # mixing allowed
    # @test BI_USD + USD == 2USD # Commented for CurrenciesBase split
    @test typeof(BI_USD + USD) == typeof(BI_USD)
    @test typeof(I128_USD + USD) == typeof(I128_USD)
    # @test BI_USD - I128_USD == 0USD # Commented for CurrenciesBase split
    # @test typeof(BI_USD - I128_USD) == typeof(BI_USD) # Commented for CurrenciesBase split
    @test BI_USD / I128_USD == 1.0
    @test typeof(BI_USD / I128_USD) == BigFloat
end

# Custom decimals
@testset "Custom Precision" begin
    flatusd = majorunit(Monetary{:USD, Int, 0})
    millusd = majorunit(Monetary{:USD, Int, 3})

    # Constructor equivalence
    @test flatusd ≡ Monetary(:USD, 1; precision=0)
    @test flatusd ≡ Monetary(:USD; precision=0)
    @test millusd ≡ Monetary(:USD, 1000; precision=3)
    @test millusd ≡ Monetary(:USD; precision=3, storage=Int)
    @test flatusd ≢ millusd

    # Custom precision arithmetic
    @test 1.111millusd + 1.222millusd == 2.333millusd
    @test 0.1flatusd == 0flatusd
    @test 0.001millusd ≠ 0.002millusd

    # Mixing comparisons allowed (surprising behaviour?)
    @test zero(flatusd) ≢ 0USD
    # @test flatusd ≥ USD # Commented for CurrenciesBase split

    # and arithmetic also, but promotes to BigInt (surprising behaviour?)
    # @test flatusd + millusd == 2USD # Commented for CurrenciesBase split
    @test typeof(flatusd + millusd) == Monetary{:USD, BigInt, 3}
    @test flatusd / millusd == 1.0
    @test typeof(flatusd / millusd) == BigFloat

    # Special metals — precision required
    @test_throws ArgumentError eval(:(@usingcurrencies XAU))
    @test_throws ArgumentError Monetary(:XAU)
    @test_throws ArgumentError Monetary(:XAU, 100)

    @test Monetary(:XAU; precision=2) ≡ majorunit(Monetary{:XAU,Int,2})
    @test Monetary(:XSU; precision=0).val == 1
end

@testset "Constructors" begin
    # the grand constructor test!
    # split up for easy debugging
    @test USD ≡ majorunit(Monetary{:USD})
    @test USD ≡ majorunit(Monetary{:USD, Int, 2})
    @test USD ≡ Monetary(:USD)
    @test USD ≡ Monetary(:USD; storage=Int)
    @test USD ≡ Monetary(:USD; precision=2)
    @test USD ≡ Monetary(:USD; storage=Int, precision=2)
    @test USD ≡ Monetary(:USD, 100)
    @test USD ≡ Monetary(:USD, 100; precision=2)

    @test USD ≢ Monetary(:USD, 1)
    @test USD ≢ Monetary(:USD; precision=0)
    @test USD ≢ Monetary(:USD; precision=4)
    @test USD ≢ Monetary(:USD; storage=BigInt)

    # and the zero test
    @test 0USD ≡ zero(Monetary{:USD})
    @test 0USD ≡ zero(Monetary{:USD,Int,2})
    @test 0USD ≡ Monetary(:USD, 0)
    @test 0USD ≡ Monetary(:USD, 0; precision=2)

    @test 0USD ≢ USD
    @test 0USD ≢ Monetary(:USD, 0; precision=0)
    @test 0USD ≢ Monetary(:USD, BigInt(0))
end

end  # testset Monetary


text:
PT_ServerIP="0.0.0.0"
PT_ServerPort="19998"
fileReg="Data^Ver00^xxxx-004^20191117041823^ID00.RAW3"
baseFolder="/media/wang/34e67bec-7d42-4d87-824f-2f98672a41af/a/bin/sjz004/"
############## or set by enviroment

baseFolder=ENV["baseFolder"]
println(baseFolder )
fileReg=ENV["fileReg"]
println(fileReg )

PT_ServerIP=ENV["PT_ServerIP"]
println(PT_ServerIP )
PT_ServerPort=ENV["PT_ServerPort"]
println(PT_ServerPort )

SLEEPSeconds1=ENV["SLEEPSeconds"]
println(SLEEPSeconds1 )
SLEEPSeconds=parse(Int ,SLEEPSeconds1)

using Dates
using Printf
using Libdl
@static if Base.Sys.islinux()
   # RawProcDLLName = joinpath( pwd(),"../Jl-Aux/dll/SVRTest.so"  )
else
   # RawProcDLLName = joinpath( pwd(),"../Jl-Aux/dll/SVRTest.dll" )
end
#print(RawProcDLLName)
ghDLLRawProc = nothing

function dllLoad()
    global ghDLLRawProc
	RawProcDLLName="SVRTest.so"
    println( RawProcDLLName)
    println()
    ghDLLRawProc = Libdl.dlopen(RawProcDLLName)
end

function Running(PT_RAWFileName ,  PT_ServerIP, PT_ServerPort,SLEEPSeconds)
    dllLoad()
    cfunObj = Libdl.dlsym( Main.ghDLLRawProc,:mainJl )
    t = ccall(cfunObj,
    Int32,
    (Cstring ,Cstring, Cstring,Cint ),
    PT_RAWFileName,  PT_ServerIP, PT_ServerPort,SLEEPSeconds)
end

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

PT_RAWFileName=joinpath(baseFolder,fileReg);
Running(PT_RAWFileName ,  PT_ServerIP, PT_ServerPort,SLEEPSeconds)

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
import Base: print

# Base means.
print(io::IO, m::EmpiricalMean) = print(io, "EmpiricalMean")

# Base kenrels.
print(io::IO, k::ZeroKernel) = print(io, "0")
print(io::IO, k::ConstantKernel) = print(io, k.c)
print(io::IO, k::EmpiricalKernel) = print(io, "EmpiricalKernel")

# Composite Kernels.
const Composite{T} = Union{CompositeMean{T}, CompositeKernel{T}, CompositeCrossKernel{T}}
function print(io::IO, shunted::Shunted{T}) where T<:Composite
    k, shunt = shunted.x, shunted.shunt
    println(io, string(shunt) * "$T")
    print_shunted_list(io, extend(shunt), k.x)
end
print(io::IO, x::Composite) = print(io, dummy_shunt(x))

const Destructurable = Union{LhsCross, RhsCross, OuterCross, OuterKernel,
                            ITMean, ITKernel,
                            ConditionalMean, ConditionalKernel, ConditionalCrossKernel,
                            FiniteMean, FiniteKernel, FiniteCrossKernel,
                            LhsFiniteCrossKernel, RhsFiniteCrossKernel,
                            DeltaSumMean, DeltaSumKernel,
                            LhsDeltaSumCrossKernel, RhsDeltaSumCrossKernel}
function print(io::IO, shunted::Shunted{T}) where T<:Destructurable
    k, shunt = shunted.x, shunted.shunt
    println(io, string(shunt) * "$T")
    print_shunted_list(io, extend(shunt), destructure_fields(k))
end
print(io::IO, k::Destructurable) = print(io, dummy_shunt(k))

destructure_fields(k::Union{LhsCross, OuterCross, OuterKernel}) = (k.f, k.k)
destructure_fields(k::RhsCross) = (k.k, k.f)
destructure_fields(k::ITMean) = (k.μ, k.f)
destructure_fields(k::ITKernel) = (k.k, k.f)
destructure_fields(μ::ConditionalMean) = (μ.μg, μ.kfg)
destructure_fields(k::ConditionalKernel) = (k.kfg, k.kgg)
destructure_fields(k::ConditionalCrossKernel) = (k.kfg, k.kfh, k.kgh)
destructure_fields(μ::FiniteMean) = (μ.μ,)
destructure_fields(
    k::Union{
        FiniteKernel,
        FiniteCrossKernel,
        LhsFiniteCrossKernel,
        RhsFiniteCrossKernel,
    },
) = (k.k,)
destructure_fields(μ::DeltaSumMean) = (μ.ϕ, μ.μ, μ.g)
destructure_fields(
    k::Union{
        DeltaSumKernel,
        LhsDeltaSumCrossKernel,
        RhsDeltaSumCrossKernel,
    },
) = (k.ϕ, k.k)

# Block stuff.
function print(io::IO, shunted::Shunted{T}) where T<:BlockMean
    μ, shunt = shunted.x, shunted.shunt
    println(io, string(shunt) * "$T")
    for (j, x) in enumerate(μ.μ)
        println(io, Shunted(extend(shunt, 2), "$j: "))
        print(io, Shunted(extend(shunt), x))
        j != length(μ.μ) && print('\n')
    end
end
print(io::IO, μ::BlockMean) = print(io, dummy_shunt(μ))

function print(io::IO, shunted::Shunted{T}) where T<:BlockCrossKernel
    k, shunt = shunted.x, shunted.shunt
    println(io, string(shunt) * "$T")
    for p in 1:size(k.ks, 1), q in 1:size(k.ks, 2)
        println(io, Shunted(extend(shunt, 2), "($q, $p): "))
        print(io, Shunted(extend(shunt), k.ks[p, q]))
        (q * p) != length(k.ks) && print('\n')
    end
end
function print(io::IO, shunted::Shunted{T}) where T<:BlockKernel
    k, shunt = shunted.x, shunted.shunt
    println(io, string(shunt) * "$T")
    for p in 1:size(k.ks_off, 1), q in 1:size(k.ks_off, 2)
        println(io, Shunted(extend(shunt, 2), "($q, $p): "))
        if p != q
            print(io, Shunted(extend(shunt), k.ks_off[p, q]))
        else
            print(io, Shunted(extend(shunt), k.ks_diag[p]))
        end
        (q * p) != length(k.ks_off) && print('\n')
    end
end
print(io::IO, μ::Union{BlockKernel, BlockCrossKernel}) = print(io, dummy_shunt(μ))


text:
#! /usr/bin/env julia

using QuadGK

function Ei(x, minfloat=1e-3, maxfloat=100)
    f = t -> exp(-t) / t  # inline function
    if x > 0
        return quadgk(f, -x, -minfloat)[1]
             + quadgk(f, minfloat, maxfloat)[1]
    else
        return quadgk(f, -x, maxfloat)[1]
    end
end

X = linspace(-1, 1, 1000)   # 1000 points
Y = [ Ei(x) for x in X ]

using Winston
plot(X, Y)
title("The function Ei(x)")
xlabel("x"); ylabel("y")
savefig("figures/Ei_integral.png")

max_stars_repo_name: pierre-haessig/slides

max_issues_repo_name: pierre-haessig/slides

max_forks_repo_name: pierre-haessig/slides


text:
#! /usr/bin/julia

# Rosetta Code, File modification time

fname = "fool.txt"
tfmt = "%FT%T%z"

println("The modification time of ", fname, " is ")
println("    ", strftime(tfmt, mtime(fname)))

println("\nTouch this file.")
touch(fname)

println("The modification time of ", fname, " is now ")
println("    ", strftime(tfmt, mtime(fname)))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Test
using REPL.TerminalMenus
using PlayingCards
using TexasHoldem
import TexasHoldem
TH = TexasHoldem

function simulate_keystrokes(keys...)
    keydict =  Dict(:up => "\e[A",
                    :down => "\e[B",
                    :enter => "\r")
    for key in keys
        if isa(key, Symbol)
            write(stdin.buffer, keydict[key])
        elseif isa(key, Char)
            write(stdin.buffer, "$key")
        else
            write(stdin.buffer, key)
        end
    end
end

@testset "Test fold" begin
    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    simulate_keystrokes(:down, :down, :enter, 'd')
    TH.player_option!(game, players[1], CheckRaiseFold())
    simulate_keystrokes(:down, :down, :enter, 'd')
    TH.player_option!(game, players[1], CallRaiseFold())
    simulate_keystrokes(:down, :down, :enter, 'd')
    TH.player_option!(game, players[1], CallAllInFold())
    simulate_keystrokes(:down, :enter, 'd')
    TH.player_option!(game, players[1], CallFold())
end

@testset "Test check" begin
    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    simulate_keystrokes(:enter, 'd')
    TH.player_option!(game, players[1], CheckRaiseFold())
end

@testset "Test call" begin
    # No initial round contribution
    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    game.table.current_raise_amt = 10.0
    simulate_keystrokes(:enter, 'd', 10.0)
    TH.player_option!(game, players[1], CallRaiseFold())
    @test last(TH.action_history(players[1])) == Call{Float64}(10.0)

    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    game.table.current_raise_amt = 10.0
    simulate_keystrokes(:enter, 'd', 10.0)
    TH.player_option!(game, players[1], CallAllInFold())
    @test last(TH.action_history(players[1])) == Call{Float64}(10.0)

    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    game.table.current_raise_amt = 10.0
    simulate_keystrokes(:enter, 'd', 10.0)
    TH.player_option!(game, players[1], CallFold())
    @test last(TH.action_history(players[1])) == Call{Float64}(10.0)
end

@testset "Test All-in raise" begin
    # No initial round contribution
    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    game.table.current_raise_amt = 10.0
    simulate_keystrokes(:down, :enter, 'd')
    TH.player_option!(game, players[1], CallAllInFold())
    @test last(TH.action_history(players[1])) == Raise{Float64}(200.0)
end

#=
Ideally, we'd use `simulate_keystrokes` for all of
these tests. We can't seem to, however.

 - Fold, Check, and Call only require prompts, then user
   inputs, so we can get away with `simulate_keystrokes`.

 - Raises require multiple interleaved user input-prompts,
   so we need to use TerminalRegressionTests.
=#

TH.use_input_io() = true
@static if !Sys.iswindows()
    using TerminalRegressionTests
    TRT = TerminalRegressionTests
    path = "terminal_test_output"

    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    table = game.table
    file = joinpath(path,"input_raise_amt.output")
    TRT.automated_test(file, ["oops_string", 10]) do emuterm
        @testset "Human player options: input_raise_amt" begin
            @test TH.input_raise_amt(game.table, players[1], emuterm) == 10
        end
    end

    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    table = game.table
    table.current_raise_amt = 0
    file = joinpath(path,"CheckRaiseFold.output")
    simulate_keystrokes(:down, :enter) # select Raise
    TRT.automated_test(file, [10]) do emuterm
        @testset "Human player options: CheckRaiseFold" begin
            TH.player_option!(game, players[1], CheckRaiseFold(), emuterm)
        end
    end

    players = (Player(Human(), 1), Player(Bot5050(), 2))
    game = Game(players)
    table = game.table
    table.current_raise_amt = 10
    file = joinpath(path,"CallRaiseFold.output")
    simulate_keystrokes(:down, :enter) # select Raise
    TRT.automated_test(file, [20]) do emuterm
        @testset "Human player options: CallRaiseFold" begin
            TH.player_option!(game, players[1], CallRaiseFold(), emuterm)
        end
    end
end
TH.use_input_io() = false


text:
test_doc(obj, doc) = @test Documented(obj).p[:docstring] == doc
test_doc_broken(obj, doc) = @test_broken Documented(obj).p[:docstring] == doc

@testset "Generated documentation" begin
    @testset "Low-level structs" begin
        test_doc_broken(
            Constructor(StructDefinition{false}(struct_by_name(:VkExtent2D))),
            """
                _Extent2D(width::Integer, height::Integer)

            Arguments:
            - `width::Integer`
            - `height::Integer`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/VkExtent2D.html)
            """,
        )

        test_doc_broken(
            Constructor(StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo))),
            """
                _InstanceCreateInfo(enabled_layer_names::AbstractArray{<:AbstractString}, enabled_extension_names::AbstractArray{<:AbstractString}; next = C_NULL, flags = 0, application_info = C_NULL)

            Arguments:
            - `enabled_layer_names::AbstractArray{<:AbstractString}`
            - `enabled_extension_names::AbstractArray{<:AbstractString}`
            - `next`: defaults to `C_NULL`
            - `flags`: defaults to `0`
            - `application_info`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/VkInstanceCreateInfo.html)
            """,
        )

        test_doc_broken(
            Constructor(StructDefinition{false}(struct_by_name(:VkSubmitInfo2KHR))),
            """
                _SubmitInfo2KHR(wait_semaphore_infos::AbstractArray{_SemaphoreSubmitInfoKHR}, command_buffer_infos::AbstractArray{_CommandBufferSubmitInfoKHR}, signal_semaphore_infos::AbstractArray{_SemaphoreSubmitInfoKHR}; next = C_NULL, flags = 0)

            Extension: VK\\\\_KHR\\\\_synchronization2

            Arguments:
            - `wait_semaphore_infos::AbstractArray{_SemaphoreSubmitInfoKHR}`
            - `command_buffer_infos::AbstractArray{_CommandBufferSubmitInfoKHR}`
            - `signal_semaphore_infos::AbstractArray{_SemaphoreSubmitInfoKHR}`
            - `next`: defaults to `C_NULL`
            - `flags`: defaults to `0`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/VkSubmitInfo2KHR.html)
            """,
        )
    end

    @testset "High-level structs" begin
        test_doc(StructDefinition{true}(struct_by_name(:VkInstanceCreateInfo)),
            """
            High-level wrapper for VkInstanceCreateInfo.

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/VkInstanceCreateInfo.html)
            """,
        )
    end

    @testset "API functions" begin
        test_doc(
            APIFunction(func_by_name(:vkEnumerateInstanceExtensionProperties), false),
            """
                enumerate_instance_extension_properties(; layer_name = C_NULL)::ResultTypes.Result{Vector{ExtensionProperties}, VulkanError}

            Return codes:
            - Error:
              - `ERROR_OUT_OF_HOST_MEMORY`
              - `ERROR_OUT_OF_DEVICE_MEMORY`
              - `ERROR_LAYER_NOT_PRESENT`

            Arguments:
            - `layer_name`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/vkEnumerateInstanceExtensionProperties.html)
            """,
        )

        test_doc(
            APIFunction(func_by_name(:vkDestroyDevice), false),
            """
                destroy_device(device; allocator = C_NULL)::Cvoid

            Arguments:
            - `device` (externsync)
            - `allocator`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/vkDestroyDevice.html)
            """,
        )

        test_doc(
            APIFunction(func_by_name(:vkGetPipelineCacheData), false),
            """
                get_pipeline_cache_data(device, pipeline_cache)::ResultTypes.Result{Tuple{UInt, Ptr{Cvoid}}, VulkanError}

            Return codes:
            - Error:
              - `ERROR_OUT_OF_HOST_MEMORY`
              - `ERROR_OUT_OF_DEVICE_MEMORY`

            Arguments:
            - `device`
            - `pipeline_cache`

            !!! warning
                The pointer returned by this function holds memory owned by Julia. It is therefore **your** responsibility to free it after use (e.g. with `Libc.free`).

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/vkGetPipelineCacheData.html)
            """,
        )

        test_doc(
            APIFunction(func_by_name(:vkWriteAccelerationStructuresPropertiesKHR), false),
            """
                write_acceleration_structures_properties_khr(device, acceleration_structures::AbstractArray{AccelerationStructureKHR}, query_type::QueryType, data_size::Integer, data::Ptr{Cvoid}, stride::Integer)::ResultTypes.Result{Result, VulkanError}

            Extension: VK\\\\_KHR\\\\_acceleration\\\\_structure

            Return codes:
            - Error:
              - `ERROR_OUT_OF_HOST_MEMORY`
              - `ERROR_OUT_OF_DEVICE_MEMORY`

            Arguments:
            - `device`
            - `acceleration_structures::AbstractArray{AccelerationStructureKHR}`
            - `query_type::QueryType`
            - `data_size::Integer`
            - `data::Ptr{Cvoid}` (must be a valid pointer with `data_size` bytes)
            - `stride::Integer`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/vkWriteAccelerationStructuresPropertiesKHR.html)
            """,
        )

        test_doc(
            APIFunction(func_by_name(:vkEnumerateInstanceLayerProperties), false),
            """
                enumerate_instance_layer_properties()::ResultTypes.Result{Vector{LayerProperties}, VulkanError}

            Return codes:
            - Error:
              - `ERROR_OUT_OF_HOST_MEMORY`
              - `ERROR_OUT_OF_DEVICE_MEMORY`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/vkEnumerateInstanceLayerProperties.html)
            """,
        )
    end
end


text:
@testset "Test Joins" begin
    y = rand(10)
    z = rand(10)

    t  = table((x=1:10,   y=y), pkey=:x)
    t2 = table((x=1:2:20, z=z), pkey=:x)

    @testset "how = :inner" begin
        t_inner = table((x = 1:2:9, y = y[1:2:9], z = z[1:5]), pkey = :x)
        @test isequal(join(t, t2; how=:inner), t_inner)
    end
    @testset "how = :left" begin
        # Missing
        z_left = Union{Float64,Missing}[missing for i in 1:10]
        z_left[1:2:9] = z[1:5]
        t_left = table((x = 1:10, y = y, z = z_left))
        @test isequal(join(t, t2; how=:left), t_left)

        # DataValue
        z_left2 = [DataValue{Float64}() for i in 1:10]
        z_left2[1:2:9] = z[1:5]
        t_left2 = table((x=1:10, y = y, z = z_left2))
        @test isequal(join(t, t2, how=:left, missingtype=DataValue), t_left2)
    end
    @testset "how = :outer" begin
        # Missing
        x_outer = union(1:10, 1:2:20)
        y_outer = vcat(y, fill(missing, 5))
        z_left = Union{Float64,Missing}[missing for i in 1:10]
        z_left[1:2:9] = z[1:5]
        z_outer = vcat(z_left, z[6:10])
        t_outer = table((x=x_outer, y=y_outer, z=z_outer); pkey=:x)
        @test isequal(join(t, t2; how=:outer), t_outer)

        # DataValue
        y_outer2 = vcat(y, fill(DataValue{Float64}(), 5))
        z_left2 = [DataValue{Float64}() for i in 1:10]
        z_left2[1:2:9] = z[1:5]
        z_outer2 = vcat(z_left2, z[6:10])
        t_outer2 = table((x=x_outer, y=y_outer2, z=z_outer2); pkey=:x)
        @test isequal(join(t, t2; how=:outer, missingtype=DataValue), t_outer2)
    end
    @testset "how = :anti" begin
        t_anti = table((x=2:2:10, y=y[2:2:10]), pkey=:x)
        @test isequal(join(t, t2; how=:anti), t_anti)
    end
    @testset "stringarray" begin
        a = table((sa=StringArray(["a","b","c"]), id=1:3))
        b = table((x=1:3, id=1:2:5))
        t_outer = table((id=[1,2,3,5], sa=["a","b","c",missing], x=[1,missing,2,3]), pkey=:id)
        t_outer2 = table((id=[1,2,3,5], sa=["a","b","c",NA], x=[1,NA,2,3]), pkey=:id)
        @test isequal(join(a, b, how=:outer, lkey=:id, rkey=:id), t_outer)
        @test isequal(join(a, b, how=:outer, lkey=:id, rkey=:id, missingtype=DataValue), t_outer2)
    end
    @testset "empty" begin
        t1 = table((key = ["a","b","c"], ))
        t2 = table((key = ["a","b"], value = [1,2]))
        res = join(t1, t2, how = :left, lkey = :key, rkey = :key)
        @test isequal(res, table((key = ["a","b","c"], value = [1,2,missing])))
    end
end


text:
using DiffEqBase, DiffEqJump
using Test, Statistics

# using Plots; plotlyjs()
doplot = false
# using BenchmarkTools
# dobenchmark = false

dotestmean   = true
doprintmeans = false

# SSAs to test
SSAalgs = (Direct(), SortingDirect(), NRM(), RSSA(), DirectCR()) #, DirectFW(), FRM(), FRMFW()),

# numerical parameters
Nsims        = 8000
tf           = 1000.0
u0           = [1,0,0,0]
expected_avg = 5.926553750000000e+02
reltol       = .01

# average number of proteins in a simulation
function runSSAs(jump_prob)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[3,end]
    end
    mean(Psamp)
end

# MODEL SETUP

# DNA repression model DiffEqBiological
# using DiffEqBiological
# rs = @reaction_network dtype begin
#     k1, DNA --> mRNA + DNA
#     k2, mRNA --> mRNA + P
#     k3, mRNA --> 0
#     k4, P --> 0
#     k5, DNA + P --> DNAR
#     k6, DNAR --> DNA + P
# end k1 k2 k3 k4 k5 k6

# model using mass action jumps
# ids: DNA=1, mRNA = 2, P = 3, DNAR = 4
reactstoch =
[
    [1 => 1],
    [2 => 1],
    [2 => 1],
    [3 => 1],
    [1 => 1, 3 => 1],
    [4 => 1]
]
netstoch =
[
    [2 => 1],
    [3 => 1],
    [2 => -1],
    [3 => -1],
    [1 => -1, 3 => -1, 4 => 1],
    [1 => 1, 3 => 1, 4 => -1]
]
spec_to_dep_jumps = [[1,5],[2,3],[4,5],[6]]
jump_to_dep_specs = [[2],[3],[2],[3],[1,3,4],[1,3,4]]
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
majumps = MassActionJump(rates, reactstoch, netstoch)

# TESTING:
prob = DiscreteProblem(u0, (0.0, tf), rates)

# plotting one full trajectory
if doplot
    plothand = plot(reuse=false)
    for alg in SSAalgs
        jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)
        sol = solve(jump_prob, SSAStepper())
        plot!(plothand, sol.t, sol[3,:], seriestype=:steppost)
    end
    display(plothand)
end

# test the means
if dotestmean
    means = zeros(Float64,length(SSAalgs))
    for (i,alg) in enumerate(SSAalgs)
        jump_prob = JumpProblem(prob, alg, majumps, save_positions=(false,false), vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)
        means[i]  = runSSAs(jump_prob)
        relerr = abs(means[i] - expected_avg) / expected_avg
        if doprintmeans
            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ", relerr)
        end

        # if dobenchmark
        #     @btime (runSSAs($jump_prob);)
        # end

        @test abs(means[i] - expected_avg) < reltol*expected_avg
    end
end

# # benchmark performance
# if dobenchmark
#     # exact methods
#     for alg in SSAalgs
#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)
#         @btime solve($jump_prob, SSAStepper())
#     end
#     println()
# end


text:
import Laplacian
import ConjugateGradient
import p2

# Run question 2 with grid size N
function q2(N::Integer, μs::Array)

    # Boundary conditions as matrices
    Ψb, B = p2.boundaryconditions2(N)

    # Reshape conditions into vectors
    ψb = reshape(Ψb, (N+1)^2)
    β  = reshape(B,  (N+1)^2)

    # Laplacian matrix
    Δ = Laplacian.lap(N)
    Δf, Δb = Laplacian.opsplit(Δ, β)
    Δe = Laplacian.deltaeff(Δf, β)

    # First guess should be zero for boundary points
    b0 = ones(Int64,(N+1)^2) .* (1 .- β)  / √sum(1.-β)

    # Compute Bμ (in parallel, of course)
    Bs = Array(Float64, length(μs))
    @sync @parallel for i in 1:length(μs)
        # Inverse iteration matrix A
        println("On μ = $(μs[i])")
        A = p2.invit(-5Δe, μs[i], β)

        # Grip it and rip it
        b = ConjugateGradient.conjgrad(A, b0)
        Bs[i] = dot(b,b)
    end

    # return the B values
    return Bs, Δe

end

# Find eigenvectors corresponding to our eigenvalues
function q2vectors(N::Integer, μ::Float64)

    # Boundary conditions as matrices
    Ψb, B = p2.boundaryconditions2(N)

    # Reshape conditions into vectors
    ψb = reshape(Ψb, (N+1)^2)
    β  = reshape(B,  (N+1)^2)

    # Laplacian matrix
    Δ = Laplacian.lap(N)
    Δf, Δb = Laplacian.opsplit(Δ, β)
    Δe = Laplacian.deltaeff(Δf, β)

    # First guess should be zero for boundary points
    ψ0 = ones(Int64,(N+1)^2) .* (1 .- β)  / √sum(1.-β)

    # Grip it and rip it; calculate ϕfree, use it to find total ψ
    A = p2.invit(-5Δe, μ, β)
    ψf = ConjugateGradient.conjgrad(A, ψ0)
    ψf /= norm(ψf)
    ψ = β.*ψb + (1.-β).*ψf

    # Reshape ψ to give the physical grid and return that value
    return reshape(ψ, N+1, N+1)

end

# Return X, Y, and Z coordinates of ϕ gridi matrix, assuming 1×1 grid size
function xyzcoords(ϕ)

    N = size(ϕ)[1] - 1
    a = 1/N

    ii = [0:N]  .* ones(N+1)'
    jj = [0:N]' .* ones(N+1)

    xs = a .* ii
    ys = a .* jj
    zs = ϕ

    X = reshape(xs, (N+1)^2)
    Y = reshape(ys, (N+1)^2)
    Z = reshape(zs, (N+1)^2)

    return X, Y, Z

end

function writexyxcoords(prefix,ϕ)

    X, Y, Z = xyzcoords(ϕ)

    writedlm("$prefix-phi-X.dat", X)
    writedlm("$prefix-phi-Y.dat", Y)
    writedlm("$prefix-phi-Z.dat", Z)

end


text: apple() = :orange


text:
# Run CRD using default parameters on the snap datasets

using MAT

include("FlowSeed.jl")
using PyCall
@pyimport localgraphclustering as lgc

function snapexp(percentofset::Int64,IndividualSeeds::Bool = false)

datasets = ["Amazon","DBLP","LiveJournal","Orkut","Youtube"]

for graph = datasets

    if IndividualSeeds
        outputstring = "Output"*string(percentofset)*"/"*graph*"_crd_individual.txt"
        outputmatrix = "Output"*string(percentofset)*"/"*graph*"_crd_individual.mat"
    else
        outputstring = "Output"*string(percentofset)*"/"*graph*"_crd.txt"
        outputmatrix = "Output"*string(percentofset)*"/"*graph*"_crd.mat"
    end

    open(outputstring,"w") do f
         write(f,"CRD. We assume we know $percentofset percent of the target cluster.\n")
    end

    # Load the edgelist of the graph
    tic()
    g = lgc.GraphLocal(homedir()*"/data/snap-top10/"*graph*".edgelist","edgelist"," ")
    toc()

    # Load the communities
    mat = matread(homedir()*"/data/snap-top10/"*graph*"-top10.mat")
    A = mat["A"]

    # Load the seed sets
    matS = matread(homedir()*"/data/snap-top10/"*graph*"-seed-starter.mat")
    if percentofset == 2
         S = matS["S2"]
    elseif percentofset == 3
         S = matS["S3"]
    elseif percentofset == 5
         S = matS["S5"]
    else
         println("Choose 2,3, or 5 percent of the target set. Defaulting to 5.")
         S = mat["S5"]
    end

    comm = mat["C"]
    Conds = mat["Conds"]
    d = sum(A,1)'
    volA = sum(nonzeros(A))

    n = size(A,1)
    d = sum(A,1)'
    println("nodes = $n")
    numcom = size(comm,2)

    crd_sets = spzeros(n,numcom)
    crd_b_sets = spzeros(n,numcom)
    crd_stats = zeros(6,numcom)
    crd_b_stats = zeros(6,numcom)

    for commID = collect(1:numcom)

        ## Choose one community to be the target cluster
        Target = find(comm[:,commID])
        TarSize = length(Target)
        TargetStats = set_stats(A,vec(Target),volA)
        Tcond = round(TargetStats[4],3)
        Tsize = length(Target)
        println("Community $commID has $Tsize nodes and a conductance of $Tcond")

        ## Get the starter for the seed set for this community
        Rstart = find(S[:,commID])
        R = neighborhood(A,Rstart,1)      # Grow by a one-hope neighborhood
        Rn = neighborhood(A,R,1)    # get the immediate neighbors of R...
        Rn = setdiff(Rn,R)          # ...but we exclude R itself
        inRc = ones(n)
        inRc[R] = 0
        Rc = find(inRc)             # complement of R
        numR = length(R)
        volR = sum(d[R])
        STATS = set_stats(A,R,volA)
        condR = round(STATS[4],4)
        pr,re,f1 = round.(PRF(Target,R),4)
        println("Seed Set: |R| = $numR, condR = $condR, PR = $pr, RE = $re, F1 = $f1")

        open(outputstring,"a") do f
             write(f,"\nCommunity $commID has $Tsize nodes and a conductance of $Tcond \n")
             write(f, "Seed Set: |R| = $numR, condR = $condR, PR = $pr, RE = $re, F1 = $f1 \n")
        end

        if IndividualSeeds
            # Run CRD trying each known target node as a seed
            # (Takes a long time and returns poor results)
            ## Next do CRD with default parameters and the full seed set
            condCRD = 1
            crd_S = []
            tic()
            for rr = Rstart
                 crd = lgc.flow_clustering(g,[rr-1],method="crd")
                 crd_temp = convert(Array{Int64},crd[1])+1
                 STATS = set_stats(A,crd_temp,volA)
                 cond_temp = round(STATS[4],3)

                 if cond_temp < condCRD
                      condCRD = cond_temp
                      crd_S = crd_temp
                 end
            end
            tCRD = toc()

            # Output some stats
            sizeCRD = length(crd_S)
            prf = PRF(Target,crd_S)
            prCRD = round(prf[1],3)
            reCRD = round(prf[2],3)
            f1CRD = round(prf[3],3)
            crd_sets[crd_S,commID] = 1
            if isnan(f1CRD)
                 f1CRD = 0
            end
            crd_stats[:,commID] = [tCRD; sizeCRD; condCRD; prCRD; reCRD; f1CRD]
            println("CRD Indiv: \tPR = $prCRD \t RE = $reCRD \t F1 = $f1CRD, Size = $sizeCRD, Cond = $condCRD \t Time = $tCRD")

            open(outputstring,"a") do f
                 write(f, "CRD Indiv: \tPR = $prCRD \t RE = $reCRD \t F1 = $f1CRD \t Size = $sizeCRD \t Cond = $condCRD \t Time = $tCRD\n")
            end

        else

            ## Run CRD with default parameters and the full seed set (best results)
            tic()
            crd = lgc.flow_clustering(g,R-1,method="crd",U=3,h=10,w=2,iterations=20)
            tCRD = toc()
            crd_S = convert(Array{Int64},crd[1])+1
            STATS = set_stats(A,crd_S,volA)
            condCRD = round(STATS[4],4)
            sizeCRD = length(crd_S)

            # Output some stats
            prf = PRF(Target,crd_S)
            prCRD = round(prf[1],3)
            reCRD = round(prf[2],3)
            f1CRD = round(prf[3],3)
            crd_sets[crd_S,commID] = 1
            if isnan(f1CRD)
                 f1CRD = 0
            end
            crd_stats[:,commID] = [tCRD; sizeCRD; condCRD; prCRD; reCRD; f1CRD]
            println("CRD Stand: \tPR = $prCRD \t RE = $reCRD \t F1 = $f1CRD, Size = $sizeCRD, Cond = $condCRD \t Time = $tCRD")

            open(outputstring,"a") do f
                 write(f, "CRD Stand: \tPR = $prCRD \t RE = $reCRD \t F1 = $f1CRD \t Size = $sizeCRD \t Cond = $condCRD \t Time = $tCRD\n")
            end
        end

    end
    matwrite(outputmatrix, Dict("crd_sets" => crd_sets,"crd_stats" => crd_stats))

end

end

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: Union{<:Rotation{2},<:RotMatrix{2}}
    # We use the awkward Union{<:Rotation{2},<:RotMatrix{2}} here rather than `Rotation{2}`
    # to make this efficient 2D method more specific than the method for general `RotMatrix`.
    T = eltype(R)
    if T == Any
        T = Float64
    end

    R(2π * rand(rng, T))
end

# A random rotation can be obtained easily with unit quaternions
# The unit sphere in R⁴ parameterizes quaternion rotations according to the
# Haar measure of SO(3) - see e.g. http://math.stackexchange.com/questions/184086/uniform-distributions-on-the-space-of-rotations-in-3d
function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: Union{<:Rotation{3},<:RotMatrix{3}}
    # We use the awkward Union{<:Rotation{3},<:RotMatrix{3}} here rather than `Rotation{3}`
    # to make this efficient 3D method more specific than the method for general `RotMatrix`.
    T = eltype(R)
    if T == Any
        T = Float64
    end

    q = QuatRotation(randn(rng, T), randn(rng, T), randn(rng, T), randn(rng, T))
    return R(q)
end

# A random rotation can be obtained via random matrix and nearest_rotation.
# This is slower than the implementation for 2 or 3 dimensions, but it is possible to create a random matrix on SO(n)
function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: RotMatrix{N} where N
    T = eltype(R)
    if T == Any
        T = Float64
    end

    m = randn(rng, SMatrix{N,N,T})
    return nearest_rotation(m)
end

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: Union{RotX,RotY,RotZ}
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return R(2π*rand(rng, T))
end

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: Union{RotXY,RotYZ,RotZX, RotXZ, RotYX, RotZY}
    T = eltype(R)
    if T == Any
        T = Float64
    end

    # Not really sure what this distribution is, but it's also not clear what
    # it should be! rand(RotXY) *is* invariant to pre-rotations by a RotX and
    # post-rotations by a RotY...
    return R(2π*rand(rng, T), 2π*rand(rng, T))
end

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{RP}) where RP <: MRP
    RP(rand(rng, QuatRotation))
end

@inline function Random.rand(rng::AbstractRNG, ::Random.SamplerType{RP}) where RP <: RodriguesParam
    RP(rand(rng, QuatRotation))
end

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{<:QuatRotation{T}}) where T
    _normalize(QuatRotation{T}(randn(rng,T), randn(rng,T), randn(rng,T), randn(rng,T)))
end