text:
@testset "  Urine PK                                                 " begin
    upkds = MetidaNCA.upkimport(upkdata, :st, :et, :conc, :vol, :subj; dosetime =  MetidaNCA.DoseTime(dose = 100))
    unca  = MetidaNCA.nca!(upkds)
    for k in keys(urefdict)
        @test  unca[1, k]  ≈ urefdict[k] atol=1E-4
    end
    #=
    @test unca[1, :Maxrate] == 4.0
    @test unca[1, :AUCall]  == 17.1250
    @test unca[1, :AUClast]  == 17.1250
    @test unca[1, :AUCpct]  ≈ 12.646069 atol=1E-6
    @test unca[1, :Tmax]    == 1.5
    @test unca[1, :AR]      == 16.0
    @test unca[1, :Vol]     == 11.0
    @test unca[1, :Prec]    == 16.0
    @test unca[1, :Rlast]    ≈ 0.333333333333 atol=1E-6
    @test unca[1, :HL]      ≈ 5.15525788157323 atol=1E-6
    @test unca[1, :ARsq]    ≈ 0.8109832410336681 atol=1E-6
    @test unca[1, :Rsq]   ≈  0.90549162 atol=1E-6
    @test unca[1, :Kel]     ≈ 0.13445441459631066 atol=1E-6
    @test unca[1, :AUCinf]  ≈ 19.60415499341648 atol=1E-6
    =#
    io = IOBuffer();
    @test_nowarn dsnca = MetidaNCA.nca!(upkds, verbose = 2, io = io)

    upkds = MetidaNCA.upkimport(upkdata, :st, :et, :conc, :vol; dosetime =  MetidaNCA.DoseTime(dose = 100))
    upkds = MetidaNCA.upkimport(upkdata[!, :st], upkdata[!, :et], upkdata[!, :conc], upkdata[!, :vol]; dosetime =  MetidaNCA.DoseTime(dose = 100))
    unca  = MetidaNCA.nca!(upkds)
    @test_nowarn show(io, upkds)
    @test_nowarn show(io, unca)
    @test_nowarn MetidaNCA.nca(upkdata, :st, :et, :conc, :vol; type = :ur, dosetime =  MetidaNCA.DoseTime(dose = 100))

    upkdatac = deepcopy(upkdata)
    upkdatac.st = float.(upkdatac.st)
    upkdatac[1, :st] = NaN
    @test_throws ErrorException MetidaNCA.upkimport(upkdatac, :st, :et, :conc, :vol, :subj; dosetime =  MetidaNCA.DoseTime(dose = 100))

    upkdatac = deepcopy(upkdata)
    upkdatac.et = float.(upkdatac.et)
    upkdatac[1, :et] = NaN
    @test_throws ErrorException MetidaNCA.upkimport(upkdatac, :st, :et, :conc, :vol, :subj; dosetime =  MetidaNCA.DoseTime(dose = 100))

    upkdatac = deepcopy(upkdata)
    upkdatac.et = float.(upkdatac.et)
    upkdatac[1, :et] = 1.5
    @test_throws ErrorException MetidaNCA.upkimport(upkdatac, :st, :et, :conc, :vol, :subj; dosetime =  MetidaNCA.DoseTime(dose = 100))
end


text:
# ---
# title: 1568. Minimum Number of Days to Disconnect Island
# id: problem1568
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-number-of-days-to-disconnect-island/description/>
# hidden: true
# ---
#
# Given a 2D `grid` consisting of `1`s (land) and `0`s (water).  An _island_ is
# a maximal 4-directionally (horizontal or vertical) connected group of `1`s.
#
# The grid is said to be **connected** if we have **exactly one  island**,
# otherwise is said **disconnected**.
#
# In one day, we are allowed to change **any** single land cell `(1)` into a
# water cell `(0)`.
#
# Return _the minimum number of days_ to disconnect the grid.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/13/1926_island.png)**
#
#
#
#     Input: grid = [[0,1,1,0],[0,1,1,0],[0,0,0,0]]
#     Output: 2
#     Explanation: We need at least 2 days to get a disconnected grid.
#     Change land grid[1][1] and grid[0][2] to water and get 2 disconnected island.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[1,1]]
#     Output: 2
#     Explanation: Grid of full water is also disconnected ([[1,1]] -> [[0,0]]), 0 islands.
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[1,0,1,0]]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[1,1,0,1,1],
#                    [1,1,1,1,1],
#                    [1,1,0,1,1],
#                    [1,1,0,1,1]]
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: grid = [[1,1,0,1,1],
#                    [1,1,1,1,1],
#                    [1,1,0,1,1],
#                    [1,1,1,1,1]]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= grid.length, grid[i].length <= 30`
#   * `grid[i][j]` is `0` or `1`.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
import SVMLightLoader: load_svmlight_file

url = "http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a1a"
filename = basename(url)

if !isfile(filename)
    import HTTPClient: get
    contents = get(url)
    f = open(filename, "w")
    write(f, contents.body.data)
    close(f)
end

# The code is from below:
# http://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/
function benchmark(target, args...)
    # Any setup code goes here.

    # Run once, to force compilation.
    println("======================= First run:")
    @time target(args...)

    # Run a second time, with profiling.
    println("\n\n======================= Second run:")

    Profile.init()
    Profile.clear()
    Profile.clear_malloc_data()

    @profile @time target(args...)

    # Write profile results to profile.bin.
    f = open("trace.def", "w")
    Profile.print(f)
    close(f)
end

benchmark(load_svmlight_file, filename)


text:
struct ReplacementNode
    text::String
    leading_trivia::String
    trailing_trivia::String
end
ReplacementNode(text::String, orig::OverlayNode{T}) where {T} = ReplacementNode(text, leading_trivia(orig), trailing_trivia(orig))
fullspan_text(r::ReplacementNode) = string(r.leading_trivia, r.text, r.trailing_trivia)

struct ChildReplacementNode
    parent
    children::Vector{Any}
    onode
    function ChildReplacementNode(a,b,c)
        # For now, always force Trivia replacement nodes to wrap child replacement nodes
        @assert !isa(c, TriviaReplacementNode)
        new(a,b,c)
    end
end
AbstractTrees.printnode(io::IO, x::ChildReplacementNode) = AbstractTrees.printnode(io, x.onode)

struct TriviaInsertionNode
    trivia::String
end

struct TriviaReplacementNode
    parent
    onode
    leading_trivia::String
    trailing_trivia::String
end
TriviaReplacementNode(parent, node::TriviaReplacementNode, onode) = TriviaReplacementNode(parent, onode, node.leading_trivia, node.trailing_trivia)
function AbstractTrees.printnode(io::IO, x::TriviaReplacementNode)
    AbstractTrees.printnode(io, x.onode)
    print(io, " [", repr(x.leading_trivia), " ", repr(x.trailing_trivia), "]")
end
Base.show(io::IO, x::Union{ChildReplacementNode, TriviaReplacementNode}) = AbstractTrees.print_tree(io, x)

function findprev_str(testf::Function, A::String, start::Integer)
    i = start
    while i >= 1
        testf(A[i]) && return i
        i = prevind(A, i)
    end
    return 0
end

function line_pos(node::OverlayNode, buffer_pos)
    find_pos = findprev_str(i->i=='\n', node.buffer, buffer_pos)
    return prevind(node.buffer, buffer_pos - find_pos)
end

function countindent_ws(ws, indents)
    i = start(ws)
    while i <= sizeof(ws)
        c = ws[i]
        i = nextind(ws, i)
        if c == '\n'
            indent = 0
            while i <= sizeof(ws)
                c = ws[i]
                !isspace(c) && break
                c == ' ' && (indent += 1)
                c == '\t' && (indent += 4)
                i = nextind(ws, i)
            end
            push!(indents, indent)
        end
    end
end

function addindent_ws(ws, nchars)
    buf = IOBuffer()
    i = start(ws)
    while i <= sizeof(ws)
        c = ws[i]
        write(buf, ws[i])
        i = nextind(ws, i)
        if c == '\n'
            # Skip whitespace until nchars is reached
            remaining_chars = -nchars
            while i <= sizeof(ws) && remaining_chars > 0
                c = ws[i]
                !isspace(c) && break
                c == ' ' && (remaining_chars -= 1)
                c == '\t' && (remaining_chars -= 4)
                i = nextind(ws, i)
            end
            while remaining_chars < 0
                write(buf, ' ')
                remaining_chars += 1
            end
        end
    end
    String(take!(buf))
end

function setindent_ws(ws, nchars)
    buf = IOBuffer()
    i = start(ws)
    while i <= sizeof(ws)
        c = ws[i]
        write(buf, ws[i])
        i = nextind(ws, i)
        if c == '\n'
            # Skip whitespace until nchars is reached
            counted_ws = 0
            while i <= sizeof(ws) && counted_ws <= nchars
                c = ws[i]
                write(buf, c)
                !isspace(c) && break
                c == ' ' && (counted_ws += 1)
                c == '\t' && (counted_ws += 4)
                i = nextind(ws, i)
            end
            # Skip remaining whitespace
            if counted_ws == nchars
                while isspace(ws[i])
                    i = nextind(ws, i)
                end
            else
                # Add whitespace until we reach the desired indent
                while counted_ws < nchars
                    write(buf, ' ')
                    counted_ws += 1
                end
            end
        end
    end
    String(take!(buf))
end

function countindent_body(expr, indents = Int[], countlast = false)
    cs = children(expr)
    for (i, c) in enumerate(cs)
        # For leaves, count indent
        if isempty(children(c)) && (countlast || i != length(cs))
            countindent_ws(trailing_trivia(c), indents)
        else
            countindent_body(c, indents, countlast || i != length(cs))
        end
    end
    indents
end

function _format_addindent_body(expr, nexpr, nindent)
    for c in children(expr)
        # For leaves, try to unindent
        if isempty(children(c))
            push!(nexpr.children, TriviaReplacementNode(next, c, leading_trivia(c),
                addindent_ws(trailing_trivia(c), nindent)))
        else
            push!(nexpr.children, format_addindent_body(c, nindent, nexpr))
        end
    end
end

function _format_setindent_body(expr, nexpr, nindent)
    for c in children(expr)
        # For leaves, try to unindent
        if isempty(children(c))
            push!(nexpr.children, TriviaReplacementNode(next, c, leading_trivia(c),
                setindent_ws(trailing_trivia(c), nindent)))
        else
            push!(nexpr.children, format_setindent_body(c, nindent, nexpr))
        end
    end
end

function format_addindent_body(expr, nindent, parent = nothing, trailing=true)
    nexpr = expr
    if !isempty(children(expr))
        nexpr = ChildReplacementNode(parent, Any[], expr)
        _format_addindent_body(expr, nexpr, nindent)
    end
    TriviaReplacementNode(parent, nexpr, addindent_ws(leading_trivia(expr), nindent), trailing ? addindent_ws(trailing_trivia(expr), nindent) : trailing_trivia(expr))
end
function format_addindent_body(expr::TriviaReplacementNode, nindent, parent = nothing, trailing=true)
    nexpr = ChildReplacementNode(nothing, Any[], expr.onode)
    _format_addindent_body(expr, nexpr, nindent)
    TriviaReplacementNode(parent, nexpr, addindent_ws(leading_trivia(expr), nindent), trailing ? addindent_ws(trailing_trivia(expr), nindent) : trailing_trivia(expr))
end

function format_setindent_body(expr, nindent, parent = nothing)
    nexpr = ChildReplacementNode(parent, Any[], expr)
    _format_setindent_body(expr, nexpr, nindent)
    nexpr
end
function format_setindent_body(expr::TriviaReplacementNode, nindent, parent = nothing)
    nexpr = ChildReplacementNode(nothing, Any[], expr.onode)
    _format_setindent_body(expr, nexpr, nindent)
    TriviaReplacementNode(parent, nexpr, leading_trivia(expr), setindent_ws(trailing_trivia(expr), nindent))
end

CSTUtils.isexpr(x::ChildReplacementNode, ::Type{S}) where {S} = isexpr(x.onode, S)
CSTUtils.isexpr(x::TriviaReplacementNode, ::Type{S}) where {S} = isexpr(x.onode, S)

AbstractTrees.children(x::ChildReplacementNode) = x.children
AbstractTrees.children(x::TriviaReplacementNode) = children(x.onode)

trailing_trivia(x::ChildReplacementNode) = trailing_trivia(last(x.children))
trailing_trivia(x::TriviaReplacementNode) = x.trailing_trivia
trailing_trivia(x::ReplacementNode) = x.trailing_trivia
trailing_trivia(x::TriviaInsertionNode) = x.trivia

leading_trivia(x::ChildReplacementNode) = leading_trivia(last(x.children))
leading_trivia(x::TriviaReplacementNode) = x.leading_trivia
leading_trivia(x::ReplacementNode) = x.leading_trivia
leading_trivia(x::TriviaInsertionNode) = ""

print_replacement(io::IO, node::OverlayNode, leading_trivia, trailing_trivia) = print(io, text(node, leading_trivia, trailing_trivia))
function print_replacement(io::IO, node::ChildReplacementNode, leading_trivia, trailing_trivia)
    length(node.children) == 0 && return
    length(node.children) == 1 && return print_replacement(io, first(node.children), leading_trivia, trailing_trivia)
    print_replacement(io, first(node.children), leading_trivia, true)
    foreach(c->print_replacement(io, c, true, true), node.children[2:end-1])
    print_replacement(io, last(node.children), true, trailing_trivia)
end
function print_replacement(io::IO, node::TriviaReplacementNode, leading_trivia, trailing_trivia)
    leading_trivia && print(io, node.leading_trivia)
    print_replacement(io, node.onode, false, false)
    trailing_trivia && print(io, node.trailing_trivia)
end
print_replacement(io::IO, node::TriviaInsertionNode, leading_trivia, trailing_trivia) = print(io, node.trivia)
function print_replacement(io::IO, node::ReplacementNode, leading_trivia, trailing_trivia)
    leading_trivia && print(io, node.leading_trivia)
    print(io, node.text)
    trailing_trivia && print(io, node.trailing_trivia)
end

print_replacement(io::IO, node) = print_replacement(io, node, false, false)


text:
using Test, CompGrids, MPI, PETSc

@testset "PETSc grids" begin
    size_t = (32,32,32)
    length = (1, 2, 3)

    for idim=1:3
        @init_backend(PETSc, Base, true);
        grid = RegularRectilinearCollocatedGrid(size=size_t[1:idim], extent=length[1:idim])
        @test grid.L[1] == 1.0
        if idim>1
            @test grid.L[2] == 2.0
        end
        if idim>2
            @test grid.L[3] == 3.0
        end
    end

    # Create grid with fields
    grid = RegularRectilinearCollocatedGrid(size=(10,20), extent=(100.0,110), fields=(T=0,P=11.22))

    # function to transfer local cartesian indice to global indices
    for i in grid.corners.lower:grid.corners.upper
        # set T based on coordinates
        grid.fields.T[grid.idx(i)] = grid.Face[1][i[1]] + grid.Face[2][i[2]]
    end

    # Check
    T = PETSc.unsafe_localarray(grid.fields.T);
    T = PETSc.reshapelocalarray(T, grid.PETSc.da)[1,:,:,:];
    Text = grid.fields.T[grid.idx(CartesianIndex(2,3,1))]
    @test T[2,3] == grid.fields.T[grid.idx(CartesianIndex(2,3,1))]
    @test T[2,3] == grid.Face[1][2] + grid.Face[2][3]

end


text:
function parallelmap{S,T}(f, A::AbstractArray{S}, B::AbstractArray{T})
    F = similar(A, promote_type(S,T), promote_shape(size(A),size(B)))
    @parallel for i in eachindex(A,B)
        @inbounds F[i] = f(A[i], B[i])
    end
    return F
end

function simplemap{S,T}(f, A::AbstractArray{S}, B::AbstractArray{T})
    F = similar(A, promote_type(S,T), promote_shape(size(A),size(B)))
    for i in eachindex(A,B)
        @inbounds F[i] = f(A[i], B[i])
    end
    return F
end
function test()
    const a = rand(10^7)
    const b = rand(10^7)

    c = parallelmap(+, a, b)
    println("parallel:")
    @time parallelmap(+, a, b)
    d = parallelmap(Base.AddFun(), a, b)
    println("functor:")
    @time parallelmap(Base.AddFun(), a, b)

    e = simplemap(+, a, b)

    println("simplemap:")
    @time simplemap(+, a, b)
    f = simplemap(Base.AddFun(), a, b)
    println("functor:")
    @time simplemap(Base.AddFun(), a, b)

    g = map(+, a, b)
    println("map:")
    @time map(+, a, b)
    h = map(Base.AddFun(), a, b)
    println("functor:")
    @time map(Base.AddFun(), a, b)

    println("plus")
    i = a+b
    @time a+b

    @assert c==d && c==e && c==f && c==g && c==h && c==i
end
test()


text:
mutable struct MPState{Y<:Tensor{T,3} where T}
    As     :: Vector{Y}
    center :: Int
end

@inline center(mps::MPState) = mps.center
@inline tensortype(::MPState{Y}) where {Y} = Y
@inline Base.eltype(::MPState{Y}) where {Y} = eltype(Y)
@inline Base.length(mps::MPState) = length(mps.As)
#@inline checkbounds(mps, l) = 0 < l <= L || throw(BoundsError(mps.vectors), l)

@inline bondspace(mps::MPState, l::Int) = dual(space(mps.As[l], 3))
@inline bonddim(mps::MPState, l::Int) = size(mps.As[l], 3)
@inline bonddim(mps::MPState) =
    [dim(mps.As[1], 1); [bonddim(mps, l) for l in 1:length(mps)]]

@inline leftspace(mps::MPState) = space(mps.As[1], 1)
@inline rightspace(mps::MPState) = bondspace(mps, length(mps))

@inline sitespace(mps::MPState, l::Int) = space(mps.As[l], 2)

@inline function Base.push!(mps::MPState{Y}, A::Y) where{Y}
    lx = length(mps)
    Vl = space(A, 1)
    if lx > 0
        bondspace(mps, lx) == Vl || throw("SpaceMismatch()")
    else
        dim(Vl) == 1 || throw("SpaceMismatch()")
    end
    push!(mps.As, A)
    mps.center = 0
    mps
end

const U1MPS{T} = MPState{SymTensor{U1,T,3}}
const MPS{T} = MPState{Array{T, 3}}

####################
### CONSTRUCTORS ###
####################
MPState{Y}() where{Y} = MPState{Y}(Vector{Y}(), 0)

# equal probability constructor
function U1MPS(::Type{T},
               lx::Int,
               d::Int,
               m::U1Charge;
               noise::T=0.0,
               seed=1911) where {T}

    mps = U1MPS{T}()
    M = m.charge

    rng = MersenneTwister(1234);
    # Just find all the possible sectors for each site and put one(T)
    # This gives an MPS with a norm of binomial(lx,m) if no noise.
    Vd = U1Space(c=>1 for c in 0:d-1)
    Vl = U1Space(0=>1)
    for site in 1:lx
        cmin = max(0, M-(d-1)*(lx-site))
        cmax = min(M, (d-1)*site)
        Vr = U1Space(c=>1 for c in cmin:cmax)
        A = SymTensor((x,y)->ones(x,y) .+ noise .* randn(rng, x,y),
                      zero(U1), (Vl, Vd, dual(Vr)))
        push!(mps, A)
        Vl = Vr
    end
    center_at!(mps, 1)
    mps
end

# constructor with intial configuration vector
function U1MPS(::Type{T},
               lx::Int,
               d::Int,
               initconf::Vector{Int}) where {T}

    @assert all((0 .<= initconf) .& (initconf .< 2))
    mps = U1MPS{T}()
    Vd = U1Space(c=>1 for c in 0:d-1)
    lchr = 0
    for site=1:lx
        Vl = U1Space(lchr=>1)
        if initconf[site] == 0
            Vr = U1Space(lchr=>1)
        else
            Vr = U1Space(lchr+1=>1)
            lchr += 1
        end
        push!(mps, fill(one(T), (Vl, Vd, dual(Vr))))
    end
    center_at!(mps, 1)
    mps
end

#constructor from a ketstate given in Ising basis
function U1MPS(lx  ::Int,
               d   ::Int,
               ketstate::Vector{T};
               svtruncation::Bool=false) where {T}

    @assert length(ketstate) == binomial(lx, div(lx,2))
    charge = div(lx,2)
    A = SymVector(U1(charge), ketstate)

    mps = U1MPS{T}()
    ℓ = 1
    ex_leg_ℓ = U1Space(0=>1)
    Vd = U1Space(0=>1, 1=>1)
    charge_range = min(charge, lx-ℓ):-1:max(0, charge-ℓ)
    smleg_ℓ = U1Space(c=>binomial(lx-ℓ, c) for c in charge_range)

    A = splitleg(A, 1, (Vd, smleg_ℓ))
    u,s,v = _svd_(A)

    push!(mps, splitleg(u, 1, (ex_leg_ℓ, Vd)))

    for ℓ = 2:lx-1
        ex_leg_ℓ = space(s, 1)
        charge_range = min(charge, lx-ℓ):-1:max(0, charge-ℓ)
        smleg_ℓ = U1Space(c=>binomial(lx-ℓ, c) for c in charge_range)
        A = fuselegs(splitleg(s*v, 2, (Vd, smleg_ℓ)), 1, 2)
        u,s,v = _svd_(A)

        push!(mps, splitleg(u, 1, (ex_leg_ℓ, Vd)))
    end
    push!(mps, splitleg(s*v, 2, (Vd, U1Space(0=>1))))
    mps
end

### conversions
###############

MPState{Y}(mps::MPState{Y}) where {Y} = mps
function MPState{Y}(mps::MPState) where {Y}
    MPState{Y}(mps.As, mps.center)
end

function convert(::Type{MPState{Y1}}, mps::MPState{Y2}) where {Y1, Y2}
    MPState{Y1}([convert(Y1, mps.As[i]) for i in eachindex(mps.As)], mps.center)
end

function MPState{Y1}(mps::MPState{Y2}) where {Y1,Y2}
    MPState{Y1}([convert(Y1, mps.As[i]) for i in eachindex(mps.As)], mps.center)
end

### TOOLS
#########

function normalize!(mps::MPState)
    if mps.center < 1
        center_at!(mps, 1)
    end
    A = mps.As[mps.center]
    if mps.center > div(length(mps), 2)
        u,s,v = _svd_(fuselegs(A, 1, 2, false))
        normalize!(s)
        mps.As[mps.center] = splitleg(u*s*v, 1, space(A)[1:2])
    else
        u,s,v = _svd_(fuselegs(A, 2, 2, true))
        normalize!(s)
        mps.As[mps.center] = splitleg(u*s*v, 2, space(A)[2:3])
    end
    mps
end

@inline function _pushleft!(mps::MPState, l::Int)
    @boundscheck 1 < l <= length(mps) ||
        throw("can't pushleft at $l of $(length(mps))")
    u,s,v = _svd_(fuselegs(mps.As[l], 2, 2))
    mps.As[l] = splitleg(v, 2, space(mps.As[l])[2:3])
    mps.As[l-1] = contract(mps.As[l-1], (1, 2, -1), u*s, (-1, 3))
    nothing
end

@inline function _pushright!(mps::MPState, l::Int)
    @boundscheck 0 < l < length(mps) ||
        throw("can't pushright at $l of $(length(mps))")
    u,s,v = _svd_(fuselegs(mps.As[l], 1, 2))
    mps.As[l] = splitleg(u, 1, space(mps.As[l])[1:2])
    mps.As[l+1] = contract(s*v, (1,-1), mps.As[l+1], (-1,2,3))
    nothing
end

function center_at!(mps::MPState, center::Int)
    x1 = 1
    xL = length(mps)
    if mps.center > 0
        x1 = mps.center
        xL = mps.center
    end
    for p=x1:center-1
        _pushright!(mps, p)
    end
    for p=xL:-1:center+1
        _pushleft!(mps, p)
    end
    mps.center = center
    mps
end

function schmidtvalues(mps::MPState)
    lx = length(mps)
    values = Vector{Vector{Float64}}(undef, lx-1)
    center_at!(mps, 1)
    A = mps.As[1]
    for l = 1:lx-1
        u,s,v = _svd_(fuselegs(A, 1, 2))
        if typeof(s) <: Diagonal
            values[l] = diag(s)
        elseif typeof(s) <: SymDiagonal
            values[l] = sort(vcat([diag(blk) for (c,blk) in s.blocks]...), rev=true)
        else
            throw("unknown type for singular values!")
        end
        A = contract(s*v, (1, -1), mps.As[l+1], (-1,2,3))
    end
    values
end

entanglementspectrum(mps::MPState) = [-1 .* log.(svs) for svs in schmidtvalues(mps)]
entanglemententropy(mps::MPState) = [entropy(x.^2) for x in schmidtvalues(mps)]

"""
    measure(mps, ops...[, xs...])

Measures the expectation values of `ops` for the given `mps`. If the
sites `xs` are not indicated then the functions measures
everywhere. The `ops` can be single site operators (matrices) or MPOs as in
```
measure(mps, mpo)
```

It returns a dictionary of position of measurements as keys and
results as values.

"""
function measure end

function measure(mps::MPState{SymTensor{S,T,3}},
                 op::SymMatrix{S,T}, l::Int) where {S,T}
    lx = length(mps)
    @assert (l <= lx) && (l > 0)
    @assert sitespace(mps, l) == space(op, 2) == dual(space(op, 1))

    center_at!(mps, l)
    A = mps.As[l]
    v = contract(A, (-1, -2, -3),
                 contract(dual(A), (1, -1, 3), op, (-1, 2)), (-1, -2, -3))
    Dict{Int, T}(l=>v)
end

function measure(mps::MPState{SymTensor{S,T,3}},
                 op::SymMatrix{S,T}) where {S,T}

    Vd = space(op, 1)
    @assert space(op, 2) == Vd

    lx = length(mps)
    values = Dict{Int, T}()
    center_at!(mps, 1)
    L = SymTensor(ones, zero(S), (U1Space(0=>1), U1Space(0=>1)))

    for x = 1:lx
        A = mps.As[x]
        v = contract(contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                              (1, -1, 3), op, (-1, 2)),
                     (-1,-2,-3), A, (-1,-2,-3))
        values[x] = v
        L = contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                     (-1, -2, 1), A, (-1,-2, 2))
    end
    values
end

function measure(mps::MPState{SymTensor{S,T,3}},
                 op1::SymMatrix{S,T},
                 op2::SymMatrix{S,T},
                 x1::Int,
                 x2::Int) where {S,T}
    lx = length(mps)
    @assert (0 < x1 < x2) && (x2 <= lx)

    Vd1 = space(op1, 1)
    @assert space(op1, 2) == Vd1
    Vd2 = space(op2, 1)
    @assert space(op2, 2) == Vd2
    (op1.charge + op2.charge != zero(S)) &&
        error("operator charges don't add up to zero! $(op1.charge), $(op2.charge)")

    values = Dict{NTuple{2, Int}, T}()
    center_at!(mps, x1)
    A = mps.As[x1]

    L = contract(contract(dual(A), (1, -1, 3), op1, (-1, 2)),
                 (-1, -2, 1), A, (-1, -2, 2))
    for x = x1+1:x2-1
        A = mps.As[x]
        L = contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                     (-1, -2, 1), A, (-1,-2, 2))
    end
    A = mps.As[x2]
    v = contract(contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                          (1, -1, 3), op2, (-1, 2)),
                 (-1,-2,-3), A, (-1,-2,-3))
    values[(x1, x2)] = v
end

function measure(mps::MPState{SymTensor{S,T,3}},
                 op1::SymMatrix{S,T},
                 op2::SymMatrix{S,T}) where {S,T}

    lx = length(mps)
    Vd1 = space(op1, 1)
    @assert space(op1, 2) == dual(Vd1)
    Vd2 = space(op2, 1)
    @assert space(op2, 2) == dual(Vd2)
    (op1.charge + op2.charge != zero(S)) &&
        error("operator charges don't add up to zero! $(op1.charge), $(op2.charge)")

    values = Dict{NTuple{2,Int}, T}()
    for x1=1:lx-1
        center_at!(mps, x1)
        A = mps.As[x1]
        L = contract(contract(dual(A), (1, -1, 3), op1, (-1, 2)),
                     (-1, -2, 1), A, (-1, -2, 2))
        for x=x1+1:lx-1
            A = mps.As[x]
            v = contract(contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                                  (1, -1, 3), op2, (-1, 2)),
                         (-1,-2,-3), A, (-1,-2,-3))
            values[(x1, x)] = v
            L = contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                         (-1, -2, 1), A, (-1,-2, 2))
        end
        A = mps.As[lx]
        v = contract(contract(contract(L, (-1, 1), dual(A), (-1, 2, 3)),
                              (1, -1, 3), op2, (-1, 2)),
                     (-1,-2,-3), A, (-1,-2,-3))
        values[x1, lx] = v
    end
    values
end

function measure(mps::MPState{SymTensor{S,T1,3}},
                 op1::SymMatrix{S,T2},
                 op2::SymMatrix{S,T3}) where {S,T1,T2,T3}
    T = promote_type(T1,T2,T3)
    measure(convert(MPState{SymTensor{S,T,3}}, mps),
            convert(SymMatrix{S,T}, op1),
            convert(SymMatrix{S,T}, op2))
end

function measure(mps::MPState{SymTensor{S,T,3}},
                 mpo::MPOperator{SymTensor{S,T,4}}) where{S,T}
    @assert length(mpo) == length(mps)
    values = Dict{Int, T}()
    dummy = VectorSpace{S}(0=>1)
    L = fill(one(T), zero(S), (dummy, dummy, dummy))

    for x=1:length(mps)
        A = mps.As[x]
        W = mpo.Ws[x]
        L = contract(contract(contract(
            L, (-1,3,4), A, (-1,2,1)),
                              (1, -1, -2, 4), W, (-2, 3, 2,-1)),
                     (1,2,-1,-2), dual(A), (-2,-1,3))
    end
    values[1] = contract(L, (-1,-2,-3), fill(one(T), (dummy,dummy,dummy)), (-1,-2,-3))
    values
end

"""
    apply!(mps, op, x[, maxdim, pushto, svnormalize])

applies the `operator` which is `d x d x d x d` sym tensor to site `x`
and `x+1` of the `mps`. The order of indeces start from the bottom
left, bottom right, top left, top right. So the counterclockwise
convention is not assumed here!

The `maxdim` operator chooses the max possible size of dimension of
the new mps at bond between `x` and `x+1`, The singular values are
push to either left `:L` or right `:R` (default) matrices using the
argument `pushto`.

"""
function apply!(mps         :: MPState{SymTensor{S,T,3}},
                op          :: SymTensor{S,T,4},
                x           :: Int;
                maxdim      :: Int=bonddim(mps, x+1),
                pushto      :: Symbol=:R,
                svnormalize :: Bool=false) where {S,T}

    @boundscheck 0 < x < length(mps) || throw()

    if mps.center < x
        center_at!(mps, x)
    elseif mps.center > x+1
        center_at!(mps, x+1)
    end

    A1 = mps.As[x]
    A2 = mps.As[x+1]

    RR = contract(A1, (1,2,-1), A2, (-1,3,4))
    R = contract(RR,
                 (1,-1,-2,4), op, (2,3,-1,-2))

    u,s,v = svdtrunc(SymMatrix(R, [1,2], [3,4]),
                     maxdim=maxdim, tol=1.e-14)

    svnormalize && normalize!(s)

    _space = space(R)
    if (pushto == :R)
        mps.As[x] = splitleg(u, 1, _space[1:2])
        mps.As[x+1] = splitleg(s*v, 2, _space[3:4])
        mps.center = x+1
    elseif (pushto == :L)
        mps.As[x] = splitleg(u*s, 1, _space[1:2])
        mps.As[x+1] = splitleg(v, 2, _space[3:4])
        mps.center = x
    else
        error("invalid push_to :", pushto)
    end
    mps
end

function apply!(mps         :: MPState{SymTensor{S,T1,3}},
                op          :: SymTensor{S,T2,4},
                x           :: Int;
                maxdim      :: Int=bonddim(mps, x+1),
                pushto      :: Symbol=:R,
                svnormalize :: Bool=false) where {S,T1,T2}

    T = promote_type(T1, T2)
    apply!(convert(MPState{SymTensor{S,T,3}}, mps),
           convert(SymTensor{S,T,4}, op),
           x, maxdim=maxdim, pushto=pushto, svnormalize=svnormalize)
end

"""
    overlap(mps1, mps2)

calculates the overlap between two matrix product states `mps1` and
`mps2` that is to run the tensor contraction corresponding to
``⟨ψ_2|ψ_1⟩``. Note that it is not divided by the norm of the two
MPSs, so it returned value is the overlap of the two states multiplied
by the norm of each.

"""
function overlap(mps1::MPState{Y},
                 mps2::MPState{Y}) where {Y}

    lx = length(mps1)
    @assert length(mps2) == lx

    A = mps1.As[1]
    B = mps2.As[1]
    L = contract(A, (-1, -2, 1), dual(B), (-1, -2, 2))

    for x=2:lx-1
        A = mps1.As[x]
        B = mps2.As[x]
        L = contract(contract(L, (-1, 1), A, (-1, 2, 3)),
                     (-1,-2,1), dual(B), (-1,-2,2))
    end
    A = mps1.As[lx]
    B = mps2.As[lx]
    v = contract(contract(L, (-1, 1), A, (-1, 2, 3)),
                 (-1,-2,-3), dual(B), (-1,-2,-3))
    v
end

function overlap(mps1::MPState{Y1}, mps2::MPState{Y2}) where {Y1,Y2}
    Y = promote_type(Y1, Y2)
    overlap(convert(MPState{Y}, mps1), convert(MPState{Y}, mps2))
end

"""
    norm(mps)

calculates the norm of a matrix product state `mps` that is to
calculate the sqrt tensor contraction corresponding to ``⟨ψ|ψ⟩``.

"""
@inline norm(mps::MPState) = sqrt(real(overlap(mps, mps)))

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
export ungap, ungap!

@doc """
    ungap!(S[, chans=CC, m=true, tap=false])
    ungap(S[, chans=CC, m=true, tap=false])

Fill time gaps in each channel of S with the mean of the channel data.

    ungap!(C[, m=true, tap=false])
    ungap(C[, m=true, tap=false])

As above for GphysChannel object C.

### Keywords
* `chans=CC`: only ungap channels `CC`.
* `m=false`: this flag fills gaps with NaNs instead of the mean.
* `tap=true`: taper data before filling gaps.

!!! warning

    If channel segments aren't in chronological order, call `merge` before using `ungap`.
""" ungap!
function ungap!(C::GphysChannel; m::Bool=true, tap::Bool=false)
  if tap
    taper!(C)
  end
  N = size(C.t,1)-2
  (N < 0 || C.fs == 0) && return nothing
  (N == 0 && C.t[2,2] == 0) && return nothing
  gapfill!(C.x, C.t, C.fs, m=m)
  proc_note!(C, string("ungap!(C, m = ", m, ", tap = ", tap, ")"),
                string("filled ", N, " gaps (sum = ",
                        sum(C.t[2:end-1, 2]), " μs)"))
  C.t = [C.t[1:1,:]; [length(C.x) 0]]
  return nothing
end

@doc(@doc ungap!)
function ungap!(S::GphysData;
  chans::ChanSpec=Int64[],
  m::Bool=true,
  tap::Bool=false)

  chans = mkchans(chans, S, keepirr=false)

  if tap
    taper!(S, chans=chans)
  end

  for i in chans
    N = size(S.t[i],1)-2
    (N < 0 || S.fs[i] == 0) && continue
    (N == 0 && S.t[i][2,2] == 0) && continue
    gapfill!(S.x[i], S.t[i], S.fs[i], m=m)
    proc_note!(S, i, string("ungap!(S, chans=", i, ", m = ", m,
                            ", tap = ", tap, ")"),
                     string("filled ", N, " gaps (sum = ",
                             sum(S.t[i][2:end-1, 2]), " μs)"))
    S.t[i] = [S.t[i][1:1,:]; [length(S.x[i]) 0]]
  end
  return nothing
end

@doc (@doc ungap!)
ungap(S::GphysChannel; m::Bool=true, tap::Bool=false) = (T = deepcopy(S); ungap!(T, m=m, tap=tap); return T)

@doc(@doc ungap!)
ungap(S::GphysData;
  chans::ChanSpec=Int64[],
  m::Bool=true,
  tap::Bool=false) = (T = deepcopy(S); ungap!(T, chans=chans, m=m, tap=tap); return T)


text:
using Statistics,Serialization,StatsPlots, DataFrames, HypothesisTests

filename_prefix = "..//..//data//working_paper//main//xi06//"

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 582 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 84
# set periods for covid plots
inf_start_day = 1
inf_end_day = 400 # cut of after period 400
inf_start_week = 1
inf_end_week = 58 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("$(filename_prefix)beta5//alpha1//covid_par_ini.jl")

### store data in txt file
io = open("stat_tests_results.txt", "w")
println(io, "Statisitcal tests: p-values from Mann-Whitney U tests\n")

# POINT A, B, C in fig 2a
paras = ["beta5//alpha1","beta5//alpha25","beta30//alpha1"]

### get point A
p=paras[1]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_A = (results -> results[:togdploss]).(worker_results)
totcas_A = (results -> results[:totcas]).(worker_results)

### get point B
p=paras[2]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_B = (results -> results[:togdploss]).(worker_results)
totcas_B = (results -> results[:totcas]).(worker_results)

### get point C
p=paras[3]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_C = (results -> results[:togdploss]).(worker_results)
totcas_C = (results -> results[:totcas]).(worker_results)

println(io,"\n TESTS FOR FIGURE 2\n")

println(io,"\n TESTS FOR GDP FIGURE 2a\n")

HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_B)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_B))
println(io,"\n pvalue for GDP A1 vs B1: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_C)
p= pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_C))
println(io,"\n pvalue for GDP A1 vs C1: ", p)

println(io,"\n \n TESTS FOR TOT CAS FIGURE 2a\n")

HypothesisTests.MannWhitneyUTest(totcas_A, totcas_B)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_A, totcas_B))
println(io,"\n pvalue for CAS A1 vs B1: ", p)

HypothesisTests.MannWhitneyUTest(totcas_A, totcas_C)
p= pvalue(HypothesisTests.MannWhitneyUTest(totcas_A, totcas_C))
println(io,"\n pvalue for CAS A1 vs C1: ", p)

### include points from xi05
filename_prefix = "..//..//data//working_paper//main//xi05//"
# POINT A, B, C in fig 2a
paras2 = ["beta5//alpha1","beta5//alpha25","beta30//alpha1"]

### get point A
p=paras2[1]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_A2 = (results -> results[:togdploss]).(worker_results)
totcas_A2 = (results -> results[:totcas]).(worker_results)

### get point B
p=paras2[2]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_B2 = (results -> results[:togdploss]).(worker_results)
totcas_B2 = (results -> results[:totcas]).(worker_results)

### get point C
p=paras2[3]
data_prefix = string("$(filename_prefix)", p , "//")
worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
gdploss_C2 = (results -> results[:togdploss]).(worker_results)
totcas_C2 = (results -> results[:totcas]).(worker_results)

### btw xi05 points

println(io,"\n \n TESTS FOR GDP FIGURE 2b\n")

HypothesisTests.MannWhitneyUTest(gdploss_A2, gdploss_B2)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A2, gdploss_B2))
println(io,"\n pvalue for GDP A2 vs B2: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_A2, gdploss_C2)
p= pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A2, gdploss_C2))
println(io,"\n pvalue for GDP A2 vs C2: ", p)

println(io,"\n \n TESTS FOR TOT CAS FIGURE 2b\n")

HypothesisTests.MannWhitneyUTest(totcas_A2, totcas_B2)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_A2, totcas_B2))
println(io,"\n pvalue for CAS A2 vs B2: ", p)

HypothesisTests.MannWhitneyUTest(totcas_A2, totcas_C2)
p= pvalue(HypothesisTests.MannWhitneyUTest(totcas_A2, totcas_C2))
println(io,"\n pvalue for CAS A2 vs C2: ", p)

### compare btw xi05 and default

println(io,"\n \n TESTS btw xi06 (Fig 2a) and xi05 (Fig 2b)\n")

HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_A2)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_A, gdploss_A2))
println(io,"\n pvalue for GDP A1 vs A2: ", p)

HypothesisTests.MannWhitneyUTest(totcas_A, totcas_A2)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_A, totcas_A2))
println(io,"\n pvalue for CAS A1 vs A2: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_B, gdploss_B2)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_B, gdploss_B2))
println(io,"\n pvalue for GDP B1 vs B2: ", p)

HypothesisTests.MannWhitneyUTest(totcas_B, totcas_B2)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_B, totcas_B2))
println(io,"\n pvalue for CAS B1 vs B2: ", p)

HypothesisTests.MannWhitneyUTest(gdploss_C, gdploss_C2)
p = pvalue(HypothesisTests.MannWhitneyUTest(gdploss_C, gdploss_C2))
println(io,"\n pvalue for GDP C1 vs C2: ", p)

HypothesisTests.MannWhitneyUTest(totcas_C, totcas_C2)
p = pvalue(HypothesisTests.MannWhitneyUTest(totcas_C, totcas_C2))
println(io,"\n pvalue for CAS C1 vs C2: ", p)

### between different starting weeks

filename_prefix = "..//..//data//working_paper//adap_start//"

paras3 = ["week1","week2","week3","week4"]

include("$(filename_prefix)//week1/covid_par_ini.jl") # necessary for some basic parameter values (nhh)

gdploss_all_means = DataFrame()
totcas_all_means = DataFrame()
pubacc_all_means = DataFrame()
badpoltime_all_means = DataFrame()
polsw_all_means = DataFrame()

for (i,p) in enumerate(paras3)
    data_prefix = string("$(filename_prefix)", p , "//")
    if p == "week2"
        data_prefix = string("..//..//data//working_paper//main//xi06//beta5//alpha1//")
    end
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss_all = (results -> results[:togdploss]).(worker_results)
    totcas_all = 100*(results -> results[:totcas]).(worker_results)
    gdploss_all_means[!, Symbol.(paras3[i])] = gdploss_all
    totcas_all_means[!, Symbol.(paras3[i])] = totcas_all
end

println(io,"\n \n TEST BETWEEN DIFFERENT STARTING WEEKS\n")

println(io,"\n TESTS FOR GDP\n ")
week1 = gdploss_all_means[!, 1]
week2 = gdploss_all_means[!, 2]
week3 = gdploss_all_means[!, 3]
week4 = gdploss_all_means[!, 4]

HypothesisTests.MannWhitneyUTest(week1, week2)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week2))
println(io,"\n pvalue for week 1 vs 2: ", p)

HypothesisTests.MannWhitneyUTest(week1, week3)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week3))
println(io,"\n pvalue for week 1 vs 3: ", p)

HypothesisTests.MannWhitneyUTest(week1, week4)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week4))
println(io,"\n pvalue for week 1 vs 4: ", p)

println(io,"\n \n TESTS FOR TOT CAS\n")

week1 = totcas_all_means[!, 1]
week2 = totcas_all_means[!, 2]
week3 = totcas_all_means[!, 3]
week4 = totcas_all_means[!, 4]

HypothesisTests.MannWhitneyUTest(week1, week2)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week2))
println(io,"\n pvalue for week 1 vs 2: ", p)

HypothesisTests.MannWhitneyUTest(week1, week3)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week3))
println(io,"\n pvalue for week 1 vs 3: ", p)

HypothesisTests.MannWhitneyUTest(week1, week4)
p = pvalue(HypothesisTests.MannWhitneyUTest(week1, week4))
println(io,"\n pvalue for week 1 vs 4: ", p)

### close txt file
close(io)


text:
## The creation of *simple* GUIs usually involves 3 steps:
##
## * create the controls
## * layout the controls
## * propagate changes in a control to some output
##
## `GtkInteract` uses the widgets of the `Interact` package, adding a few, including output widgets
## `GtkInteract` uses layout containers using a style defined by the `Escher` package.
## `GtkInteract` uses the `Reactive` package to propagate changes from controls.

using GtkInteract, Reactive, Plots
immerse()

## To create controls is done by calling the constructor. The
## constructors are fairly primitive and have little customization.
## The available controls
## basically map some data into a GUI element:

n = slider(1:10, label="n")
rb = radiobuttons(["one", "two", "three"], label="radio")
cb = checkbox(true, label="checkbox")

## the `Interact.widget` function will try to read your mind:

Interact.widget(1:4) # a slider
Interact.widget(true) # a check box
Interact.widget(["one", "two", "three"]) # togglebuttons

## To layout the controls we first make a parent window, and then add
## widgets to that.  We have two types of parent windows `mainwindow`
## and `window`. The `mainwindow` container uses a fixed "form" layout
## for simplicity. One simply appends the controls to the
## mainwindow. The `label` property of the widget, is used to label
## the control.

w = mainwindow(title="Simple test");
n = slider(1:10, label="n")
m = slider(11:20, label="m")
cg = cairographic()
append!(w, [n,m,cg])
w

## Other layouts are available. The style is borrowed from Escher,
## though not all the functionality from Escher is implemented. For
## this use, instead of `mainwindow`, a `window` object is used.
## For example, we have:

n = slider(1:10, label="n:")
m = slider(11:20, label="m:")
cg = cairographic()
w = window(hbox(vbox(hbox(label(n.label), n),
                     hbox(label(m.label), m)
                    ),
               padding(5, grow(cg))); title = "Simple test")

## This uses boxes (`hbox` and `vbox`) to orangize child widgets. The
## structure above is that a `hbox` is used to hold a box that has the
## two slider controls and the graphic window. A vertical box is used
## to layout the slider controls. Unlike the layout with `mainwindow`,
## the labels must be managed by the programmer. Hence the construct
## `hbox(label(n.label), n)`. This pattern packs a label next to the
## control.
##
## The layout widgets are just containers. Child widgets are added at
## construction time, as illustrated above, or through the `push!` and
## `append!` methods.
##
## The `window` container is similar. It differs in rendering
## though. Once it is `display`ed, the child widgets are realized and
## can no longer have children appended.
##
##  To adjust layouts, there are a few
## attributes. Illustrated above is `padding`, which adds 5 pixels of
## "padding" around the widget used to display the graphic; and `grow`
## which instructs the child to grow to fill any allocated space.

## However the layout is done, once the widgets are made and layed
## out, they can be connected using Reactive signals:

map(n, m) do n,m
    push!(cg, plot(sin, n*pi, m*pi))
end

## the map function propagates changes to the underlying widgets to
## the function call.
##
## The `map` functon (`map(fn, widgets...)`) takes
## the widgets and allows them to be referenced by their values. That
## is, the use of `n` and `m` within the `plot` command uses the
## values in the controls `n` and `m`.
##
## With `GtkInteract`, there are
## two basic types of widgets: input widgets and output widgets. The
## graphics device is an output widget. For these, values are
## `push!`ed onto them. So the call `push!(cg, ...)` should update the
## graphic using the generated plot. Other output widgets include
## `immersefigure`, and `textarea`.

## A pattern where the update only is related to a button push, and
## not each control is desired if the update is expensive to compute
## (With sliders it may be computed up to 4 times during the move).

w = mainwindow(title="Simple test");   # suppress call to `display` with semicolon
n = slider(1:10, label="n")
m = slider(1:10, label="m")
btn = button("update")
append!(w, [n, m, btn])
w                                       # call display() to create window

## We can connect to the button to pass along the values of the other widget:
map(btn) do _
    println(Reactive.value(n.signal) * Reactive.value(m.signal))
end

## The above is cumbersome, as to extract a value from a widget
## requires grabbing the value from its
## signal property. This package provides an shorthant: `value(n)`.
##
## For a more elegant solution, we can leverage this pattern from
## Shasi. It is modifed from
## https://groups.google.com/forum/?fromgroups#!topic/julia-users/Ur5b2_dsJyA
## as there are changes needed for the newer Reactive.jl.
##

w = mainwindow(title="Simple test")
n = slider(1:10, label="n")
m = slider(1:10, label="m")
btn = button("update")
append!(w, [n, m, btn])
w

vals = map(tuple, n, m)  # not just (n,m), as map "lifts" values.

map(vals->println(join(vals, ", ")), Reactive.sampleon(btn.signal, vals))

## The `vals` is not just `(n,m)`, but instead `map` lifts the signals
## from the widgets `n` and `m` and passes their values on to a tuple
## whenever the signals changes. So `vals` reflects the current state
## of the two widgets.
##
## The `sampleon` function samples the value of `vals` when the
## button's signal is emitted, which happens when the button is
## clicked. So `Reactive.sampleon(btn.signal, vals)` makes a new
## signal which is emitted when the button is clicked and passes on a
## tuple of values reflecting the state of the widgets `m` and `n`.
##
## In this example, the values are just printed. In the next example a
## graph is drawn.

using Reactive, GtkInteract, Plots
immerse()

α = slider(1:10, label="α")
β = slider(1:10, label="β")
replot = button("replot")
cg = cairographic()

## display
window(hbox(vbox(hbox(label(α.label),α),
                 hbox(label(β.label), β),
                 replot),
            padding(5, cg));
       title="layout")

       ## Our action:
function draw_plot(α, β)
    push!(cg, plot(x -> α + sin(x + β), 0, 2pi))
end

## We can then connect the button as follows:
coeffs = Reactive.sampleon(replot.signal, map(tuple, α, β))
map(vals -> draw_plot(vals...), coeffs)

## Here, unlike with `@manipulate`, the graphic is only updated when
## the `replot` button is pressed. Changing the sliders only updates
## the `coeffs` values, which is propogated when the button is
## pressed.

## For this specific pattern with a button, `GtkInteract` extends Reactive's `sampleon` function, so that
## the first line can be just:
coeffs = sampleon(replot, α, β)

## kitchen sink of control widgets
##
using DataStructures
a = OrderedDict{Symbol, Int}()
a[:one]=1; a[:two] = 2; a[:three] = 3

l = Dict()
l[:slider]        = slider(1:10, label="slider")
l[:button]        = button("button label")
l[:checkbox]      = checkbox(true, label="checkbox")
l[:togglebutton]  = togglebutton(true, label="togglebutton")
l[:dropdown]      = dropdown(a, label="dropdown")
l[:radiobuttons]  = radiobuttons(a, label="radiobuttons")
l[:select]        = selectlist(a, label="select")  # aka Interact.select
l[:togglebuttons] = togglebuttons(a, label="togglebuttons")
l[:buttongroup]   = buttongroup(a, label="buttongroup") # non exclusive
l[:textbox]       = textbox("text goes here", label="textbox")

w = mainwindow()
append!(w, values(l))

## Another output widget
## progress bar
##
using GtkInteract
w = mainwindow();                       # don't call display until all done
b = button("press")
pb = progress()

append!(w, [pb, b])
w
# connect button press to update
map(b) do _
    push!(pb, floor(Integer, 100*rand()))
end

## There is some some support for toolbars and menubars. This silly example
## demonstrates the basic usage

btn = button("button")
btn1 = button("one")
btn2 = button("two")
tog = togglebutton(value=true, label="button")

mb = menu(menu(btn1, separator(),
               menu(btn1, btn2, label="submenu"), # submenus need `label=` argument
               btn2, label="one"),
          menu(btn1, btn2, btn, label="two"))

map(_ -> println("ouch"), btn1)         # some actin
map(_ -> println("kjapow"), btn2)

tb = toolbar(icon("window-close", btn), # add icon...
             separator(),
             (tog |> icon("window-open"))
             )

l = button("one")

window(mb, tb, grow(icon("window-close", l)))

## We see that a toolbar just takes child widgets as its
## arguments. The pattern `icon(icon_name, item)` adds an icon to the
## toolbar item.

## The menubar is similar to the toolbar. We add submenus using `menu`
## again. These submenus must have a `label` value specified through a
## keyword argument. (These labels show up on the toolbar.)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    @define_scenario(def)

Define a scenario type compatible with StochasticPrograms using the syntax
```julia
@define_scenario name = begin
    ...structdef...

    [@zero begin
        ...
        return zero(scenario)
    end]

    [@expectation begin
        ...
        return expected(scenarios)
     end]
end
```
The generated type is referenced through `name` and a default constructor is always generated. This constructor accepts the keyword `probability` to set the probability of the scenario occuring. Otherwise, any internal variables and specialized constructors are defined in the @define_scenario block as they would be in any Julia struct.

If possible, a `zero` method and an `expected` method will be generated for the defined type. Otherwise, or if the default implementation is not desired, these can be user provided through [`@zero`](@ref) and [`@expectation`](@ref).

The defined scenario type will be available on all Julia processes.

## Examples

The following defines a simple scenario ``ξ`` with a single value.

```jldoctest
@define_scenario ExampleScenario = begin
    ξ::Float64
end

ExampleScenario(1.0, probability = 0.5)

# output

ExampleScenario with probability 0.5
  ξ: 1.0

```

See also: [`@zero`](@ref), [`@expectation`](@ref), [`@sampler`](@ref)
"""
macro define_scenario(arg)
    @capture(arg, scenarioname_Symbol = scenariodef_) || error("Invalid syntax. Expected: scenarioname = begin scenariodef end")
    vars = Vector{Symbol}()
    vartypes = Vector{Union{Expr,Symbol}}()
    vardefs = Vector{Expr}()
    zerodefs = Vector{Expr}()
    expectdefs = Vector{Expr}()
    def = postwalk(prettify(scenariodef)) do x
        if @capture(x, var_Symbol::vartype_)
            push!(vars, var)
            push!(vartypes, vartype)
            push!(vardefs, x)
            return x
        end
        if @capture(x, @zero zerodef_)
            push!(zerodefs, zerodef)
            return @q begin end
        end
        if @capture(x, @expectation expectdef_)
            push!(expectdefs, expectdef)
            return @q begin end
        end
        return x
    end
    # Handle zero definition
    if length(zerodefs) > 1
        error("Only provide one zero implementation")
    end
    provided_zerodef, zerodef = if length(zerodefs) == 1
        zerodef = zerodefs[1]
        true, zerodef
    else
        false, Expr(:block)
    end
    # Handle expectation definition
    if length(expectdefs) > 1
        error("Only provide one expectation implementation")
    end
    provided_expectdef, expectdef = if length(expectdefs) == 1
        expectdef = expectdefs[1]
        expectdef = postwalk(prettify(expectdef)) do x
            @capture(x, constructor_Symbol(args__)) && constructor == scenarioname && return :(StochasticPrograms.ExpectedScenario($x))
            return x
        end
        true, expectdef
    else
        false, Expr(:block)
    end
    # Prepare automatic expectation definition
    combine = Expr(:call)
    push!(combine.args, scenarioname)
    push!(combine.args, :($(Expr(:parameters, :(probability = 1.0)))))
    for var in vars
        push!(combine.args, :(probability(s1)*s1.$var + probability(s2)*s2.$var))
    end
    # Prepare automatic scenariotext definition
    textdef = Expr(:block)
    for var in vars
        key = Meta.quot(var)
        push!(textdef.args, :(print(io, "\n  $($key): $(scenario.$var)")))
    end
    # Define scenario type
    code = @q begin
        if StochasticPrograms.supports_expected([$(vartypes...)], $provided_expectdef) || $provided_expectdef
            if StochasticPrograms.supports_zero([$(vartypes...)], $provided_zerodef) || $provided_zerodef
                struct $scenarioname <: AbstractScenario
                    probability::Probability
                    $def
                    function (::Type{$scenarioname})($(vardefs...); probability::AbstractFloat = 1.0)
                        return new(Probability(probability), $(vars...))
                    end
                end
                if $provided_zerodef
                    function Base.zero(::Type{$scenarioname})
                        $zerodef
                    end
                else
                    function Base.zero(::Type{$scenarioname})
                        return $scenarioname(zero.([$(vartypes...)])...; probability = 1.0)
                    end
                end
                if $provided_expectdef
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        $expectdef
                    end
                else
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        return StochasticPrograms.ExpectedScenario(reduce(scenarios) do s1, s2
                                                                   $combine
                                                                   end)
                    end
                end
                function StochasticPrograms.scenariotext(io::IO, scenario::$scenarioname)
                    $textdef
                    return io
                end
            else
                @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function zero(::Type{{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
            end
        else
            @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function expected(scenarios::Vector{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @zero(def)

Define the additive zero scenario inside a @scenario block using the syntax:
```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines a zero scenario for the example scenario defined in [`@define_scenario`](@ref)

```julia
@zero begin
    return ExampleScenario(0.0)
end
```

See also [`@define_scenario`](@ref)
"""
macro zero(def) @warn "@zero should be used inside a @define_scenario block." end
"""
    @expectation(def)

Define how to form the expected scenario inside a [`@define_scenario`](@ref) block. The scenario collection is accessed through the reserved keyword `scenarios`.

```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines expectation for the example scenario defined in [`@scenario`](@ref)

```julia
@expectation begin
    return ExampleScenario(sum([probability(s)*s.ξ for s in scenarios]))
end
```

See also [`@define_scenario`](@ref)
"""
macro expectation(def) @warn "@expectation should be used inside a @scenario block." end


text:
module Compatibility_test

using SpikeSorting

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

num_channels=1;
(buf,nums)=output_buffer(num_channels);

function make_voltage(num_channels)

    v=rand(1:1000, 1000, num_channels);
    count = 1
    for i=525:550
        v[i,1] = v[i,1] + count
        count = count + 10
    end
    for i=551:775
        v[i,1]+=1
        count -= 10
    end
    v
end

v=make_voltage(num_channels)

detects=subtypes(SpikeSorting.Detect)
aligns=subtypes(SpikeSorting.Align)
features=subtypes(SpikeSorting.Feature)
reduces=subtypes(SpikeSorting.Reduction)
clusters=subtypes(SpikeSorting.Cluster)
thresholds=subtypes(SpikeSorting.Threshold)

for d in detects
    for a in aligns
        for f in features
            for r in reduces
                for c in clusters
                    for t in thresholds

                        s1=create_multi(d(),c(),a(),f(),r(),t(),num_channels)

                        cal!(s1,v,buf,nums,0)
                        cal!(s1,v,buf,nums,1)
                        cal!(s1,v,buf,nums,2)
                        onlinesort!(s1,v,buf,nums)

                            @test s1[1].sigend == v[end-74:end,1]

                    end
                end
            end
        end
    end
end

end


text:
using GithubMarkdown
using Test

@testset "GithubMarkdown.jl" begin
    @testset "html" begin
        io = IOBuffer()
        rendergfm(io, joinpath(@__DIR__, "test1.md"); format = "html")
        @test String(take!(io)) == read("test1.html", String)

        rendergfm("testout.html", joinpath(@__DIR__, "test1.md"); format = "html")
        @test read("testout.html", String) == read("test1.html", String)
    end

    @testset "latex" begin
        io = IOBuffer()
        rendergfm(io, joinpath(@__DIR__, "test1.md"); format = "latex")
        @test String(take!(io)) == read("test1.tex", String)
    end

    @testset "error handling" begin
        @test_throws ErrorException rendergfm("dontwriteme", "idontexist.md")
        @test !isfile("dontwriteme")
        @test_throws ArgumentError rendergfm("dontwriteme", joinpath(@__DIR__, "test1.md"); format = "nahhh")
        @test !isfile("dontwriteme")
    end
end


text:
# routine tests for wrappers
# This section is just to check existence and sanity of methods.
# Serious accuracy verification is elsewhere.

using LinearAlgebra
using LinearAlgebra: sorteig!
using Test
using GenericSchur

# allow for different phase factors and perhaps ordering
# v1 is the reference value against which we are checking v2
function _chkeigvecs(v1::AbstractMatrix{Tref}, v2::AbstractMatrix{T},
                     isnrml, sameorder=true
                     ) where {Tref,T}
    m,n = size(v2)
    if isnrml
        @test norm(v2' * v2 - I) < 50 * n * eps(real(T))
    end
    tmat = v1' * v2
    zT = zero(real(T))
    nvecs_ok = 0
    Tcheck = (precision(real(Tref)) < precision(real(T))) ? real(Tref) : real(T)
    myeps = eps(Tcheck)
    for j in 1:m
        # largest and next-largest projections
        if sameorder
            p1 = abs(tmat[j,j])
            tmat[j,j] = 0
            p2 = maximum(abs.(tmat[:,j]))
        else
            p1,p2 = zT,zT
            for i in 1:m
                ta = abs(tmat[i,j])
                if ta > p1
                    p1 = ta
                elseif ta > p2
                    p2 = ta
                end
            end
        end
        δ = abs(p1-1)
        ok = (δ < sqrt(myeps))
        if (verbosity[] > 0) && (δ > 100*myeps)
            @warn "vector projection error is $δ angle $(acos(clamp(p1,0,1)))"
        end
        if isnrml
            ok &= (abs(p2) < sqrt(myeps))
        end
        nvecs_ok += ok
    end
    @test nvecs_ok == m
end

let T = BigFloat
    n = 10
    # FIXME: should protect against accidental poor condition
    A = rand(T,n,n)
    for (w,f) in zip([:bare, :hermitian, :symmetric],[identity, Hermitian, Symmetric])
        @testset "wrappers $w $T" begin
            Awrk = f(A)
            E = eigen(Awrk)
            @test norm(Awrk*E.vectors - E.vectors * Diagonal(E.values)) < sqrt(eps(T))
            v = eigvecs(Awrk)
            _chkeigvecs(E.vectors, v, w != :bare)
            λ = eigvals(Awrk)
            if (verbosity[] > 0) && (! (λ ≈ E.values))
                @warn "eigval order comparison (eigvals, eigen, diff): "
                    display(hcat(λ, E.values, λ .- E.values))
                    println()
            end
            @test λ ≈ E.values
        end
    end
end

let T = Complex{BigFloat}
    n = 10
    # FIXME: should protect against accidental poor condition
    A = rand(T,n,n)
    for (w,f) in zip([:bare, :hermitian],[identity, Hermitian])
        @testset "wrappers $w $T" begin
            Awrk = f(A)
            E = eigen(Awrk)
            @test norm(Awrk*E.vectors - E.vectors * Diagonal(E.values)) < sqrt(eps(real(T)))
            v = eigvecs(Awrk)
            _chkeigvecs(E.vectors, v, w != :bare)
            λ = eigvals(Awrk)
            @test λ ≈ E.values
            if (verbosity[] > 0) && (! (λ ≈ E.values))
                @warn "eigval order comparison (eigvals, eigen, diff): "
                    display(hcat(λ, E.values, λ .- E.values))
                    println()
            end
        end
    end
end

function _chkrcond(x,y,rtol,atol=1e3*eps(eltype(x)))
    ok = true
    for (xi,yi) in zip(x,y)
        x1 = clamp(xi,0,1)
        y1 = clamp(yi,0,1)
        if abs(x1 - y1) > rtol * abs(y1) + atol
            ok = false
        end
    end
    if (verbosity[] > 0 && !ok) || (verbosity[] > 1)
        if !ok
            @warn "rcond mismatch (trial, ref, ratio): "
        else
            @info "rcond comparison (trial, ref, ratio): "
        end
        display(hcat(x, y, x ./ y))
        println()
    end
    return ok
end

# This section tests compatibility with Julia PR #38483 which was reverted.
# It should be rewritten but is kept to facilitate study of questionable cases.
if ((VERSION > v"1.7.0-DEV.976") && (VERSION < v"1.7.0-beta3.37"))
    for (T,Tref) in ((Complex{BigFloat},ComplexF64),(BigFloat, Float64))
        @testset "extended eigen() $T" begin
            n = 10
            Aref = rand(Tref,n,n)
            A = T.(Aref)
            old = precision(real(T))
            @assert old >= 53
            if T <: Real
                Eref1 = eigen(Aref, jvl=true)
                # we convert to complex for condition nrs, so this is a fairer check
                Eref2 = eigen(Aref .+ 0im, jvl=true, jce=true, jcv=true)
            else
                Eref1 = eigen(Aref, jvl=true, jce=true, jcv=true)
                Eref2 = Eref1
            end
            setprecision(real(T), 53) do
                E = eigen(A, jvl=true, jce=true, jcv=true)
                @test norm(E.vectorsl' * A - Diagonal(E.values) * E.vectorsl') < sqrt(eps(real(T)))
                if verbosity[] > 0
                    @info "eigval comparison (trial, ref, diff): "
                    display(hcat(E.values, Eref1.values, E.values .- Eref1.values))
                    println()
                end
                @test sorteig!(E.values) ≈ sorteig!(Eref1.values)
                # stdlib inverts condition nrs for real matrices (WTF?)
                # if we switch to the real forms, we would need this:
                # rconde = (Tref <: Real) ? (1.0 ./ Eref1.rconde) : Eref1.rconde
                # rcondv = (Tref <: Real) ? (1.0 ./ Eref1.rcondv) : Eref1.rcondv
                rconde = Eref2.rconde
                rcondv = Eref2.rcondv
                @test _chkrcond(E.rconde, rconde, 0.1)
                # CHECKME: we should not need to be so lenient.
                # The nasty cases are not common (empirically 1 vector in < 10% of 10x10),
                # but nightlys on Github are good at finding them.
                # This accuracy is sufficient for typical applications.
                @test _chkrcond(E.rcondv, rcondv, 5.0)
                _chkeigvecs(Eref1.vectors, E.vectors, false)
                _chkeigvecs(Eref1.vectorsl, E.vectorsl, false)
            end
        end
    end
end


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "Differentials" begin
    p = [1.,0.,0.]
    q = [0.,1.,0.]
    M = Sphere(2)
    X = log(M,p,q)
    @testset "Differentials on Sn(2)" begin
        @test differential_log_basepoint(M,p,p,X) == -X
        @test differential_log_basepoint(M,p,q,X) == -X
        @test differential_log_argument(M,p,p,X) == X
        @test differential_log_argument(M,p,q,X) == zero_tangent_vector(M,q)
        @test differential_exp_basepoint(M,p,zero_tangent_vector(M,p),X) == X
        @test norm(M,q, differential_exp_basepoint(M,p,X,X) - [-π/2, 0., 0.]) ≈ 0 atol=6*10^(-16)
        @test differential_exp_argument(M,p,zero_tangent_vector(M,p),X) == X
        @test norm(M,q,differential_exp_argument(M,p,X,zero_tangent_vector(M,p))) ≈ 0
        for t in [0,0.15,0.33,0.66,0.9]
            @test differential_geodesic_startpoint(M,p,p,t,X) == (1-t)*X
            @test norm(M,p,differential_geodesic_endpoint(M,p,p,t,X) - t*X)  ≈ 0 atol=10.0^(-16)
        end
    end
    @testset "Differentials on Power of Sn(2)" begin
        N = PowerManifold(M, NestedPowerRepresentation(), 3)
        x = [p,q,p]
        y = [p,p,q]
        V = [X, zero_tangent_vector(M,p), -X]
        @test norm(
                N,
                x,
                differential_forward_logs(N,x,V)
                - [-X, [π/2, 0., 0.],zero_tangent_vector(M,p)]
            ) ≈ 0 atol=8*10.0^(-16)
        @test differential_log_argument(N, x, y, V) == [ V[1], V[2], V[2] ]
    end
    @testset "Differentials on SPD(2)" begin
        #
        # Single differentials on Hn
        M2 = SymmetricPositiveDefinite(2)
        p2 = [1. 0.; 0. 1.]
        X2 = [0.5 1.;1. 0.5]
        q2 = exp(M2,p2,X2)
        # Text differentials (1) Dx of Log_xy
        @test norm(M2, p2, differential_log_basepoint(M2, p2, p2, X2) + X2) ≈ 0 atol=4*10^(-16)
        @test norm(M2, q2, differential_log_argument(M2, p2, q2, zero_tangent_vector(M2,p2))) ≈ 0 atol=4*10^(-16)
        @test norm(M2, p2, differential_exp_basepoint(M2, p2, zero_tangent_vector(M2,p2), X2) - X2) ≈ 0 atol=4*10^(-16)
        @test norm(M2, p2, differential_exp_argument(M2,p2,zero_tangent_vector(M2,p2),X2) - X2) ≈ 0 atol=4*10^(-16)
        for t in [0,0.15,0.33,0.66,0.9]
            @test norm(M2, p2, differential_geodesic_startpoint(M2, p2, p2, t, X2) - (1-t)*X2 ) ≈ 0 atol=4*10^(-16)
            @test norm(M2, p2, differential_geodesic_endpoint(M2, p2, p2, t, X2) - t*X2) ≈ 0 atol=4*10.0^(-16)
        end
        @test norm(M2, q2, differential_geodesic_startpoint(M2, p2, q2, 1., X2)) ≈ 0 atol=4*10.0^(-16)
        @test norm(M2, q2, differential_exp_basepoint(M2, p2, X2, zero_tangent_vector(M2,p2) )) ≈ 0 atol=4*10.0^(-16)
        @test norm(M2,q2,differential_exp_argument(M2,p2,X2,zero_tangent_vector(M2,p2))) ≈ 0 atol=4*10.0^(-16)
        # test coeff of log_basepoint, since it is not always expicitly used.
        @test βdifferential_log_basepoint(-1.0, 1.0, 2.0) ≈ -2*cosh(2.0)/sinh(2.0)
    end
    @testset "Differentials on Euclidean(2)" begin
        M3 = Euclidean(2)
        x3 = [1., 2.]
        ξ3 = [1.,0.]
        @test norm(M3, x3, differential_exp_basepoint(M3, x3, ξ3, ξ3) - ξ3) ≈ 0 atol=4*10.0^(-16)
    end
    @testset "Differentials on the Circle" begin
        M = Circle()
        p = 0
        q = π/4
        X = π/8
        @test differential_log_argument(M,p,q,X) == X
    end
end


text:
module MLJLibs

using MLJ
using MLJScikitLearnInterface,
    MLJDecisionTreeInterface,
    MLJClusteringInterface,
    MLJGLMInterface,
    MLJLIBSVMInterface,
    MLJMultivariateStatsInterface,
    MLJXGBoostInterface
export MLJ, MLJScikitLearnInterface,MLJDecisionTreeInterface,MLJClusteringInterface,MLJGLMInterface,MLJLIBSVMInterface,MLJMultivariateStatsInterface,MLJXGBoostInterface

end

max_stars_repo_name: andreakropp/datarobot-user-models

max_issues_repo_name: andreakropp/datarobot-user-models

max_forks_repo_name: andreakropp/datarobot-user-models


text:
const logfact = [0,[sum(log(i) for i=1:n) for n=1:50]...];

"""
    f_logfact(n)

The function returns logarithm of n!
"""
function f_logfact(n)
    (n < 0 || n > 50) && error("n < 0 || n > 50. Modify if needed.")
    @inbounds return logfact[n+1]
end

# special function used in the code:
# (two_n/2)! for even numbers
# (two_n-1/2)! for odd numbers
const logfact2 = [f_logfact(div(two_n,2)) for two_n=1:100];
function f_logfact2(two_n)
    (two_n < 0 || two_n > 100) && error("two_n < 0 || two_n > 100. Modify if needed.")
    @inbounds return logfact2[two_n+1]
end


text:
using MLStyle
using MLStyle.Record
export trans

@as_record Token
@as_record Struct_mainprogram
@as_record Struct_ifstmt
@as_record Struct_gate
@as_record Struct_gatedecl
@as_record Struct_decl
@as_record Struct_barrier_ids
@as_record Struct_reset
@as_record Struct_measure
@as_record Struct_iduop
@as_record Struct_u
@as_record Struct_cx
@as_record Struct_idlist
@as_record Struct_mixedlist
@as_record Struct_argument
@as_record Struct_explist
@as_record Struct_pi
@as_record Struct_fnexp
@as_record Struct_neg
@as_record Struct_bin

struct Session
    n     :: Ref{Int}
end

function declare_q!(sess::Session, reg::Symbol, amount::Int)
    rng(a, b) = a+1 : a + b
    qbits = collect(rng(sess.n.x, amount))
    sess.n.x += amount
    :($reg = Int[$(qbits...)])
end

macro const_eval(x)
    __module__.eval(x)
end

fake_line_node = LineNumberNode(1, "fake")

get_qbits(sess) = Expr(:macrocall, Symbol("@eval"), fake_line_node, :($sess.n.x))

const compare_bits(a :: BitArray{1}, b :: I) where I <: Integer =
    begin
        for (bit, ai) in enumerate(a)
            cur = 2 ^ (bit - 1)
            if ((cur & b) !== 0) !== ai
                return false
            end
        end
        true
    end

init_sess() = Session(Ref(0))
trans(qasm) = trans(qasm, init_sess())
function trans(qasm, sess::Session)
    function app(op, args...)
        args = map(rec, args)
        op = Symbol(op)
        :($op($(args...)))
    end

    rec(qasm) = trans(qasm, sess)
    @match qasm begin
        Struct_pi(_) => Base.pi
        Token{:id}(str=str) => Symbol(str)
        Token{:real}(str=str) => parse(Float64, str)
        Token{:nninteger}(str=str) => parse(Int64, str)
        Struct_neg(value=value) => :(-$(rec(value)))
        Struct_bin(l=l, op=Token(str=op), r=r) => app(op, l, r)
        Struct_fnexp(fn = Token(str=fn), arg=arg) =>
            let fn = @match fn begin
                        "sin" => sin
                        "cos" => cos
                        "tan" => tan
                        "exp" => exp
                        "ln"  => log
                        "sqrt"=> sqrt
                        _     => error("not impl yet")
                    end
                app(fn, arg)
            end

        Struct_idlist(hd=Token(str=hd), tl=nothing) => [Symbol(hd)]
        Struct_idlist(hd=Token(str=hd), tl=tl) => [Symbol(hd), rec(tl)...]

        Struct_explist(hd=hd, tl=nothing) => [rec(hd)]
        Struct_explist(hd=hd, tl=tl) => [rec(hd), rec(tl)...]

        Struct_mixedlist(hd=hd, tl=nothing) => [rec(hd)]
        Struct_mixedlist(hd=hd, tl=tl) => [rec(hd), rec(tl)...]

        Struct_argument(id=Token(str=id), arg=nothing) => Symbol(id)
        Struct_argument(id=Token(str=id), arg=Token(str=int)) =>
            let ind = parse(Int, int) + 1 # due to julia 1-based index
                :($(Symbol(id))[$ind])
            end

        Struct_cx(out1=out1, out2=out2) =>
            let ref1 = rec(out1),
                ref2 = rec(out2)
                :($CX($(get_qbits(sess)), $ref1, $ref2))
            end
        Struct_u(in1=in1, in2=in2, in3=in3, out=out) =>
            let (a, b, c) = map(rec, (in1, in2, in3)),
                ref = :($(rec(out))[1])
                :($U($(get_qbits(sess)), $ref, $a, $b, $c))
            end
        Struct_iduop(gate_name = Token(str=gate_name), args=nothing, outs=outs) =>
            let refs = rec(outs),
                gate_name = Symbol(gate_name)
                :($gate_name((), $(refs...)))
            end
        Struct_iduop(gate_name = Token(str=gate_name), args=exprlist, outs=outs) =>
            let refs = rec(outs),
                exprs = Expr(:tuple, rec(exprlist)...),
                gate_name = Symbol(gate_name)
                :($gate_name($exprs, $(refs...)))
            end
        Struct_gate(
            decl = Struct_gatedecl(
                id=Token(str=fid),
                args= nothing && Do(args=[]) ||
                      idlist  && Do(args = rec(idlist)),
                outs=outs
            ),
            goplist = nothing && Do(goplist=[]) ||
                      goplist && Do(goplist = map(rec, goplist))
         ) =>
            let out_ids :: Vector{Symbol} = rec(outs),
                fid = Symbol(fid)
                quote
                    function $fid(($(args...), ), $(out_ids...))
                        chain($(goplist...))
                    end
                end
            end

        Struct_decl(
            regtype = Token(str=regtype),
            id = Token(str=id),
            int = Token(str = n)
        ) =>
            let id = Symbol(id),
                n = parse(Int, n)

                if regtype == "qreg"
                    declare_q!(sess, id, n)
                else
                    :($id = BitArray(undef, $n))
                end
            end
        Struct_ifstmt(l=Token(str=l), r=r, body=body) =>
            let l = Symbol(l),
                r = rec(r),
                body = rec(body)
                :(if $compare_bits($l, $r); $body end)
            end
        Struct_mainprogram(
            prog = stmts
        ) =>
            let stmts = map(rec, stmts)
                Expr(:block, stmts...)
            end
    end
end


text:
@testset "Sudoku" begin

@testset "Sudoku from opensourc.es" begin
    com = CS.init()

    grid = zeros(Int, (9,9))
    grid[1,:] = [0,2,1,0,7,9,0,8,5]
    grid[2,:] = [0,4,5,3,1,0,0,0,9]
    grid[3,:] = [0,7,0,0,4,0,0,1,0]
    grid[4,:] = [0,0,0,1,0,8,0,3,6]
    grid[5,:] = [0,6,0,0,0,0,2,0,8]
    grid[6,:] = [0,0,0,0,0,3,0,0,4]
    grid[7,:] = [6,0,8,0,0,0,0,0,0]
    grid[8,:] = [0,9,4,0,0,7,8,0,0]
    grid[9,:] = [2,0,0,5,0,0,0,4,0]

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "Hard sudoku" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 5 4 6 0 0 9]
    grid[2,:] = [0 2 0 0 0 0 0 0 7]
    grid[3,:] = [0 0 3 9 0 0 0 0 4]
    grid[4,:] = [9 0 5 0 0 0 0 7 0]
    grid[5,:] = [7 0 0 0 0 0 0 2 0]
    grid[6,:] = [0 0 0 0 9 3 0 0 0]
    grid[7,:] = [0 5 6 0 0 8 0 0 0]
    grid[8,:] = [0 1 0 0 3 9 0 0 0]
    grid[9,:] = [0 0 0 0 0 0 8 0 6]

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "Hard sudoku infeasible" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 5 4 6 0 0 9]
    grid[2,:] = [0 2 0 0 0 0 0 0 7]
    grid[3,:] = [0 0 3 9 0 0 0 0 4]
    grid[4,:] = [9 0 5 0 0 0 0 7 3]
    grid[5,:] = [7 0 0 0 0 0 0 2 0]
    grid[6,:] = [0 0 0 0 9 3 0 0 0]
    grid[7,:] = [0 5 6 0 0 8 0 0 0]
    grid[8,:] = [0 1 0 0 3 9 0 0 0]
    grid[9,:] = [0 0 0 0 0 0 8 0 6]

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Infeasible
    @test !fulfills_sudoku_constr(com_grid)
end

@testset "Hard fsudoku repo" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 0 0 0 0 0 0]
    grid[2,:] = [0 1 0 6 2 0 0 9 0]
    grid[3,:] = [0 0 2 0 0 9 3 1 0]
    grid[4,:] = [0 0 4 0 0 6 0 8 0]
    grid[5,:] = [0 0 8 7 0 2 1 0 0]
    grid[6,:] = [0 3 0 8 0 0 5 0 0]
    grid[7,:] = [0 6 9 1 0 0 4 0 0]
    grid[8,:] = [0 8 0 0 7 3 0 5 0]
    grid[9,:] = [0 0 0 0 0 0 0 0 0]

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "Hard fsudoku repo 0-8" begin
    com = CS.init()

    grid = zeros(Int,(9,9))
    grid[1,:] = [0 0 0 0 0 0 0 0 0]
    grid[2,:] = [0 1 0 6 2 0 0 9 0]
    grid[3,:] = [0 0 2 0 0 9 3 1 0]
    grid[4,:] = [0 0 4 0 0 6 0 8 0]
    grid[5,:] = [0 0 8 7 0 2 1 0 0]
    grid[6,:] = [0 3 0 8 0 0 5 0 0]
    grid[7,:] = [0 6 9 1 0 0 4 0 0]
    grid[8,:] = [0 8 0 0 7 3 0 5 0]
    grid[9,:] = [0 0 0 0 0 0 0 0 0]
    grid .-= 1

    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        if val == -1
            if ind == 81 # bottom right
                # some other values are possible there
                com_grid[ind] = add_var!(com, 9, 11)
            elseif ind == 80 # one above (will be 9 in the end)
                com_grid[ind] = add_var!(com, 7, 11)
            else
                com_grid[ind] = add_var!(com, 0, 8)
            end
        else
            com_grid[ind] = add_var!(com, 0, 8; fix=val)
        end
    end

    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
end

@testset "top95 some use backtracking" begin
    grids = sudokus_from_file("data/top95")
    c = 0
    for grid in grids
        com = CS.init()

        com_grid = create_sudoku_grid!(com, grid)
        add_sudoku_constr!(com, com_grid)

        @test solve!(com) == :Solved
        @test fulfills_sudoku_constr(com_grid)
        c += 1
    end
    # check that actually all 95 problems were tested
    @test c == 95
end

@testset "Number 7 in top95.txt w/o backtracking" begin
    com = CS.init()

    grid = Int[6,0,2,0,5,0,0,0,0,0,0,0,0,0,3,0,4,0,0,0,0,0,0,0,0,0,0,4,3,0,0,0,8,0,
              0,0,0,1,0,0,0,0,2,0,0,0,0,0,0,0,0,7,0,0,5,0,0,2,7,0,0,0,0,0,0,0,0,0,
              0,0,8,1,0,0,0,6,0,0,0,0,0]
    grid = transpose(reshape(grid, (9,9)))

    com_grid = create_sudoku_grid!(com, grid)
    add_sudoku_constr!(com, com_grid)

    solve!(com; backtrack=false)
    @test !com.info.backtracked
    @test com.info.backtrack_fixes == 0
    @test com.info.in_backtrack_calls == 0
    @show com.info
end

end


text:
mkpath("output/out_fh")
outdir="output/out_fh/"

using Bridge, StaticArrays, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using DataFrames
using CSV

T = 2.0
dt = 1/5000
τ(T) = (x) ->  x * (2-x/T)
tt = τ(T).(0.:dt:T)

sk = 0 # skipped in evaluating loglikelihood

generate_data = false
νHparam = false

# settings in case of νH - parametrisation
ϵ = 10^(-3)
Σdiagel = 10^(-10)#10^(-10)

# settings sampler
# iterations = 5*10^4
# skip_it = 1000
# subsamples = 0:skip_it:iterations

# specify observation scheme
L = @SMatrix [1. 0.]
Σdiagel = 10^(-10)
Σ = @SMatrix [Σdiagel]

# specify target process
struct FitzhughDiffusion <: ContinuousTimeProcess{ℝ{2}}
    ϵ::Float64
    s::Float64
    γ::Float64
    β::Float64
    σ::Float64
end

Bridge.b(t, x, P::FitzhughDiffusion) = ℝ{2}((x[1]-x[2]-x[1]^3+P.s)/P.ϵ, P.γ*x[1]-x[2] +P.β)
Bridge.σ(t, x, P::FitzhughDiffusion) = ℝ{2}(0.0, P.σ)
Bridge.constdiff(::FitzhughDiffusion) = true

P = FitzhughDiffusion(0.1, 0.0, 1.5, 0.8, 0.3) # Ditlevsen-Samson
#P = FitzhughDiffusion(0.1, 0.0, 1.5, 0.8, 1.0)
x0 = ℝ{2}(-0.5, -0.6)

# Generate Data
if generate_data
     include("/Users/Frank/Sync/DOCUMENTS/onderzoek/code/diffbridges/truepaths_fh.jl")
end

# specify auxiliary process
struct FitzhughDiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    ϵ::Float64
    s::Float64
    γ::Float64
    β::Float64
    σ::Float64
    t::Float64
    u::Float64
    T::Float64
    v::Float64
end

function uv(t, P::FitzhughDiffusionAux)
    λ = (t - P.t)/(P.T - P.t)
    P.v*λ + P.u*(1-λ)
end

for k1 in (1:3)
    for k2 in (1:2)
        Random.seed!(4)# this is what i used all the time
        Random.seed!(44)
        aux_choice = ["linearised_end" "linearised_startend"  "matching"][k1]
        endpoint = ["first", "extreme"][k2]

        # settings sampler
        iterations =  !(k1==3) ? 5*10^4 : 10*10^4
        skip_it = 1000
        subsamples = 0:skip_it:iterations
        printiter = 100

        if endpoint == "first"
            #v = -0.959
            v = -1
        elseif endpoint == "extreme"
            #v = 0.633
            v = 1.1
        else
            error("not implemented")
        end

        if aux_choice=="linearised_end"
            Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ-3*P.v^2/P.ϵ  -1/P.ϵ; P.γ -1.0]
            Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ+2*P.v^3/P.ϵ, P.β)
            ρ = endpoint=="extreme" ? 0.9 : 0.0
        elseif aux_choice=="linearised_startend"
            Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ-3*uv(t, P)^2/P.ϵ  -1/P.ϵ; P.γ -1.0]
            Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ+2*uv(t, P)^3/P.ϵ, P.β)
            ρ = endpoint=="extreme" ? 0.98 : 0.0
        else
            Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ  -1/P.ϵ; P.γ -1.0]
            Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ-(P.v^3)/P.ϵ, P.β)
            ρ = 0.99
        end

        Bridge.σ(t, P::FitzhughDiffusionAux) = ℝ{2}(0.0, P.σ)
        Bridge.constdiff(::FitzhughDiffusionAux) = true

        Bridge.b(t, x, P::FitzhughDiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
        Bridge.a(t, P::FitzhughDiffusionAux) = Bridge.σ(t,P) * Bridge.σ(t, P)'

        Pt = FitzhughDiffusionAux(P.ϵ, P.s, P.γ, P.β, P.σ, tt[1], x0[1], tt[end], v)

        # Solve Backward Recursion
        Po = νHparam ? Bridge.PartialBridgeνH(tt, P, Pt, L, ℝ{1}(v),ϵ, Σ) : Bridge.PartialBridge(tt, P, Pt, L, ℝ{1}(v), Σ)

        ####################### MH algorithm ###################
        # initalisation
        W = sample(tt, Wiener())
        X = solve(Euler(), x0, W, P)
        Xo = copy(X)
        solve!(Euler(),Xo, x0, W, Po)

        solve!(Euler(),X, x0, W, Po)
        ll = llikelihood(Bridge.LeftRule(), X, Po,skip=sk)

        # further initialisation
        Wo = copy(W)
        W2 = copy(W)
        XX = Any[]
        if 0 in subsamples
            push!(XX, copy(X))
        end

        acc = 0

        for iter in 1:iterations
            # Proposal
            sample!(W2, Wiener())

            Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy

            #    i0 = length(W.yy)÷4
            #    Wo.yy[1:i0] .= W.yy[1:i0]
            #    Wo.yy[i0+1:end] .= W.yy[i0] + ρ*(W.yy[i0+1:end] - W.yy[i0+1]) + sqrt(1-ρ^2)*(W2.yy[i0+1:end]-W2.yy[i0+1])

            solve!(Euler(),Xo, x0, Wo, Po)

            llo = llikelihood(Bridge.LeftRule(), Xo, Po)

            if mod(iter, printiter) == 0
                print(iter," ll $ll $llo, diff_ll: ",round(llo-ll, digits=3))#, X[10], " ", Xo[10])
            end
            accept = false
            if log(rand()) <= llo - ll
                X.yy .= Xo.yy
                W.yy .= Wo.yy
                ll = llo
                accept = true
                acc +=1
            end

            if iter in subsamples
                push!(XX, copy(X))
            end
            if mod(iter, printiter) == 0
                accept && print("✓")
                println()
            end
        end

        @info "Done."*"\x7"^6

        # write mcmc iterates to csv file

        fn = outdir*"iterates-"*endpoint*"-"*aux_choice*".csv"
        f = open(fn,"w")
        head = "iteration, time, component, value \n"
        write(f, head)
        iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:50:length(X), (i,s) in enumerate(subsamples) ][:]
        #writecsv(f,iterates)
        writedlm(f,iterates,',')
        close(f)

        ave_acc_perc = 100*round(acc/iterations, digits=2)

        # write info to txt file
        fn = outdir*"info-"*endpoint*"-"*aux_choice*".txt"
        f = open(fn,"w")
        write(f, "Choice of auxiliary process: ",aux_choice,"\n")
        write(f, "Choice of endpoint: ",endpoint,"\n\n")
        write(f, "Number of iterations: ",string(iterations),"\n")
        write(f, "Skip every ",string(skip_it)," iterations, when saving to csv","\n\n")
        write(f, "Starting point: ",string(x0),"\n")
        write(f, "End time T: ", string(T),"\n")
        write(f, "Endpoint v: ",string(v),"\n")
        write(f, "Noise Sigma: ",string(Σ),"\n")
        write(f, "L: ",string(L),"\n\n")
        write(f,"Mesh width: ",string(dt),"\n")
        write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n")
        write(f, "Average acceptance percentage: ",string(ave_acc_perc),"\n")
        close(f)

       println("Average acceptance percentage: ",ave_acc_perc,"\n")

    end
end

println("Parametrisation of nu and H? ", νHparam)


text:
module HW5

using DMUStudent
using Compose
using StaticArrays
using CommonRLInterface
using Random: GLOBAL_RNG, MersenneTwister, SamplerTrivial, AbstractRNG
using Obfuscatee
using ProgressMeter: @showprogress
using JSON
using IntervalSets

export mc

Base.rand(rng::AbstractRNG, s::SamplerTrivial{<:AbstractInterval}) = minimum(s[]) + rand(rng)*(maximum(s[]) - minimum(s[]))

@binclude(".bin/hw5_eval")

mc = MountainCar()

end


text:
function main()
    order = [50,48,46,44,42,40,38,36,34,32,30,31,33,35,37,39,41,43,45,47,49]
    l = order[1] + order[21] + sum(sqrt(200 * (order[i] + order[i+1] - 50)) for i in 1:20)
    println(Int(round(l * 1000)))
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
function tcm_transactionsToCounts(N, transactions)
	T = size(transactions,1)
	assortments = unique(transactions[:,1:N], 1)
	M = size(assortments,1)

	transaction_counts = zeros(Int64,M,N)

	for t = 1 : T
	    for m in 1:M
	    	if (transactions[t,1:N] == assortments[m,:] )
	    		transaction_counts[m, transactions[t,N+1] ] += 1
	    	end
	    end
	end

	return assortments, transaction_counts
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using GerryChain
using LightGraphs
using PyPlot
using LinearAlgebra
using StatsBase

"""
    load_graph_from_edge_list(graph_file_path::String, weight_file_path::String)

Construct basegraph object using a graph file and a weight file.
Graph file describes basic structure of graph using edge list.
Weight file stores population assigned to each node

"""
function load_graph_from_edge_list(graph_file_path::String, weight_file_path::String)

    graph, edge_list = open(graph_file_path) do file
        edge_list = []
        nodes = Set()
        for ln in eachline(file)
            ln = split(ln)
            first = parse(Int64,ln[1])
            second = parse(Int64,ln[2])
            push!(edge_list,[first,second])
            push!(nodes, first)
            push!(nodes, second)
        end
        if nodes != Set(1:length(nodes))
            throw(ArgumentError("Error loading input graph. Graph must be connected"))
        end
        graph = SimpleGraph(length(nodes))
        return graph, edge_list
    end

    for i in edge_list
        add_edge!(graph, i[1], i[2])
    end

    edge_src, edge_dst = GerryChain.edges_from_graph(graph)
    adj_matrix = GerryChain.adjacency_matrix_from_graph(graph)
    neighbors = GerryChain.neighbors_from_graph(graph)

    total_population, populations = open(weight_file_path) do file
        node_populations = []
        weights = read(file, String)
        weights = split(weights)
        for i in weights
            push!(node_populations, parse(Int64,i))
        end
        return sum(node_populations), node_populations
    end

    return BaseGraph(
    nv(graph),
    ne(graph),
    total_population,
    populations,
    adj_matrix,
    edge_src,
    edge_dst,
    neighbors,
    graph,
    [],
    )
end

"""
    load_partition_from_line(graph, assignments::String)::Partition

Load a partition from one assignment
Numbers in the assignment string correspond to the part that node is assigned to
Use BaseGraph generated by load_graph_from_edge_list()

"""
function load_partition_from_line(graph, assignments::String)::Partition
    populations = graph.populations
    if occursin(",", assignments)
        input = split(assignments, ",")
        weight = input[1]
        assignments_string = input[2:end]
        for i in eachindex(assignments_string)
            assignments_string[i] = strip(assignments_string[i], [' ', '[', ']'])
        end
    else
        assignments_string = split(assignments)
    end
    assignments = zeros(Int64, length(assignments_string))
    districts = []

    for i in eachindex(assignments_string)
        district_id = parse(Int64, assignments_string[i])
        assignments[i] = district_id
        if !(district_id in districts)
            push!(districts, district_id)
        end
    end

    sort!(districts)

    # districts are reassigned to 1 -> n-1 following the original order
    reassign_district = Dict(zip(districts, 1:length(districts)))

    for i in eachindex(assignments)
        assignments[i] = reassign_district[assignments[i]]
    end
    num_districts = length(districts)

    # get cut_edges, district_adjacencies
    dist_adj, cut_edges = GerryChain.get_district_adj_and_cut_edges(graph, assignments, num_districts)
    num_cut_edges = sum(cut_edges)

    # get district populations
    dist_populations =
        get_district_populations(assignments, populations, graph.num_nodes, num_districts)
    # get district_nodes
    dist_nodes = get_district_nodes(assignments, graph.num_nodes, num_districts)
    # return Partition with no parent by default
    return Partition(
        num_districts,
        num_cut_edges,
        assignments,
        dist_populations,
        cut_edges,
        dist_adj,
        dist_nodes,
        nothing,
    )

end

"""
    spanning_tree_distribution(graph, enum_file::String)::Dict

Calculate spanning tree distribution of a full enumeration

"""
function spanning_tree_distribution(graph, enum_file::String)
    cut_edges = Dict{Int64,Float64}()
    lines = 0
    cut_edges, total_score = open(enum_file) do file
        for line in eachline(file)
            lines += 1
            partition = load_partition_from_line(graph,line)
            score = 1
            for district in partition.dist_nodes
                vlist = Vector{Int64}()
                for k in district
                    push!(vlist, k)
                end
                subgraph, vector = LightGraphs.induced_subgraph(graph.simple_graph, vlist)
                L = LightGraphs.laplacian_matrix(subgraph)
                L = L[1:end .!= 1, 1:end .!= 1]
                spanning_tree = det(L)
                score *= spanning_tree
            end
            if partition.num_cut_edges in keys(cut_edges)
                cut_edges[partition.num_cut_edges] += score
            else
                cut_edges[partition.num_cut_edges] = score
            end
        end
        return cut_edges, sum(values(cut_edges))
    end
    return Dict(k => v / total_score for (k, v) in cut_edges)
end

"""
    calculate_benchmark(graph, benchmark_type::String)::Dict

Calculate the desired benchmark (cut edges, etc.) for input enumeration
Output a dictionary of summary statistics

"""
function calculate_benchmark(graph, benchmark_type)::Dict
    @assert benchmark_type == "cut_edges" "Only the `cut_edges` benchmark is currently implemented"
    result = Dict{Int64,Int64}()
    counter = 0
    for ln in eachline(stdin)
        counter += 1
        if counter > 5
            break
        end
        partition = load_partition_from_line(graph, ln)
        cut_edge = partition.num_cut_edges
        if cut_edge in keys(result)
            result[cut_edge] += 1
        else
            result[cut_edge] = 1
	    end
    end
    sum_values = sum(values(result))
    return Dict(k => v / sum_values for (k, v) in result)
end

"""
    calculate_Kullback_Leibler_distance(target, test)

Calculate the Kullback Leibler distance between the target distribution and the testing distribution
Input: dictionary output from calculate_benchmark and spanning_tree_distribution

"""
function calculate_Kullback_Leibler(target, test)
    p = sort(collect(keys(target)))
    q = collect(keys(test))
    final_q = zeros(Int64, length(p))
    dummy_value = 0.0001
    for (index, value) in enumerate(p)
        if value in q
            final_q[index] = value
        else
            final_q[index] = dummy_value
        end
    end
    for value in q
        if !(value in p)
            push!(final_q, value)
            push!(p, dummy_value)
        end
    end
    kl = StatsBase.kldivergence(p, final_q)
    return kl
end

"""
    plot_distribution(results, color)

Plot the summary benchmarks that we need.
Results is a dictionary with count of each benchmark level

"""
function plot_distribution(results, color = "maroon")
    pygui(true)
    x = [key for key in keys(results)]
    y = [value for value in values(results)]
    bar(x, y, color = color, width = 0.4, alpha = 0.3)
    title("Cut Edges")
end

max_stars_repo_name: mggg/benchmarks

max_issues_repo_name: mggg/benchmarks

max_forks_repo_name: mggg/benchmarks


text:
const WILDCARD_REPR='␣'

struct MillString{T <: Optional{BitVector}, S <: AbstractString} <: AbstractString
    s::S
    m::T

    MillString(s::S, m::Nothing) where S <: AbstractString = new{Nothing, S}(s, m)
    function MillString(s::S, m::BitVector) where S <: AbstractString
        @assert length(s) == length(m)
        ns = _replace(s, m)
        new{BitVector, S}(ns, m)
    end
end

MillString(s::AbstractString) = MillString(s, nothing)

Flux.@forward MillString.s Base.length, Base.firstindex, Base.lastindex, Base.nextind, Base.reverse,
    Base.ncodeunits, Base.codeunit

macro mill_str(s::String)
    quote
        sc = collect($s)
        m = isequal.(WILDCARD_REPR, sc)
        if any(m)
            MillString(_replace!(sc, m), m)
        else
            MillString($s, nothing)
        end
    end
end

function _replace!(s::AbstractVector{<:Char}, m)
    s[m] .= Char(wildcard_code())
    String(s)
end
_replace(s::AbstractString, m) = _replace!(collect(s), m)

Base.show(io::IO, s::MillString{Nothing}) = print(io, s.s)

function Base.show(io::IO, s::MillString)
    print(io, "\"")
    for (i, si) in eachindex(s.s) |> enumerate
        if s.m[i]
            print(io, WILDCARD_REPR)
        else
            print(io, s.s[si])
        end
    end
    print(io, "\"")
end

Base.codeunit(s::MillString, i::Integer) = codeunit(s.s, i)