text:
module CQL

####################################################################
# CQLConnection
####################################################################

type CQLConnection
  socket  :: Base.TCPSocket
  buffer  :: IOBuffer
  msg_id  :: UInt8
  replies :: Dict
  pending :: Int
end

####################################################################
# Helpers
####################################################################

function ticks()
	return round(Int, time()*1000)
end

function die(;pid=getpid())
	ccall( (:kill, "libc"), Int32, (Int32,Int32), pid, 9) #Need this to force Julia to not hang
end

####################################################################
# Connect & Disconnect
####################################################################

function connect(server::AbstractString = "localhost", port::Int = 9042)
  con = CQLConnection(Base.connect(server, port),
                      IOBuffer(), 1, Dict(), 0);
  sendMessage(con, 0x01, Dict("CQL_VERSION" => "3.0.0"));
  version, id, opcode, len = readServerMessage(con.socket);
  con.pending = 0;
  @async handleServerMessages(con);
  @assert version == 0x82 || version == 0x84 || version == 0x83
  @assert opcode  == 0x02

  con
end

function disconnect(con::CQLConnection)
  while 0 < con.pending
    yield();
  end
  close(con.socket);
  con.socket  = Base.TCPSocket();
  con.buffer  = IOBuffer();
  con.msg_id  = 1;
  con.replies = Dict();
  con.pending = 0;
  con
end

####################################################################
# Handle Server Messages
####################################################################

function readServerMessage(socket::Base.TCPSocket)
  version = read(socket, UInt8);
  flags   = read(socket, UInt8);
  id      = read(socket, UInt16);
  opcode  = read(socket, UInt8);
  len     = ntoh(read(socket, UInt32));
  (version, id, opcode, len)
end

function handleServerMessage(con::CQLConnection)
  version, id, opcode, len = readServerMessage(con.socket);
  if id > 0 then
    put!(pop!(con.replies, id),
         (opcode, readbytes(con.socket, len)));
  elseif opcode == 0x00 then
    kind   = ntoh(read(con.socket, Int32));
    strlen = ntoh(read(con.socket, Int16));
    print("ERROR (HANDLE) [$kind] : "*
            bytestring(readbytes(con.socket, strlen)));
    throw(InterruptException())
  else
    for _ in 1:len
      read(con.socket, UInt8);
    end
  end
  con.pending -= 1;
  nothing
end

function handleServerMessages(con::CQLConnection)

  while !eof(con.socket)
    try
      handleServerMessage(con);
    catch err

      if !isa(err, EOFError) then
	close(con.socket)
	yield()
	#die() #bit harsh forcing thread to die - moved outside, might work with spawnat if julia hangs again
	Base.throwto(current_task(), InterruptException())
      else
	nothing
      end

    end
  end
  nothing
end

####################################################################
# Timestamp
####################################################################

type Timestamp
  milliseconds::Int64
end

const timeOffset = time(Libc.strptime("%F %T %z %Z",
                                 "1970-01-01 00:00:00 +0000 UTC"));

function Timestamp(str::AbstractString)
  Timestamp(1000 * int64(time(Libc.strptime("%F %T", str)) + timeOffset))
end

function Base.print(io::IO, t::Timestamp)
  print(io, "ts\"");
  print(io, strftime("%F %T", div(t.milliseconds, 1000)
                              - timeOffset));
  print(io, "\"");
  nothing
end

show(io::IO, t::Timestamp) = print(io, t);

macro ts_str(p)
  Timestamp(p)
end

####################################################################
# Decoding
####################################################################

function decodeString(s::IO)
  strlen = ntoh(read(s, Int16));
  bytestring(readbytes(s, strlen));
end

function decodeResultSubColumn(s::IO, typ)
  nrOfBytes = ntoh(read(s, Int32));
  decodeValue(s, nrOfBytes, (typ, nothing, nothing))
end

function decodeList(s::IO, typ)
  nrOfElements = ntoh(read(s, UInt32));
  ar = Array(Any, nrOfElements);
  for ix in 1:nrOfElements
    ar[ix] = decodeResultSubColumn(s, typ);
  end
  ar
end

function decodeMap(s::IO, val_type)
  Set(decodeList(s, val_type))
end

function decodeDict(s::IO, key_type, val_type)
  d = Dict();
  nrOfElements = Int(ntoh(read(s, UInt32)));
  for i in 1:nrOfElements
    key = decodeResultSubColumn(s, key_type);
    val = decodeResultSubColumn(s, val_type);
    d[key] = val;
  end
  d
end

function decodeValue(s::IO, nrOfBytes::Integer, types)
  type_kind, val_type, key_type = types;
  nrOfBytes < 0     ? nothing :
  type_kind == 0x02 ? ntoh(read(s, UInt64)) :
  type_kind == 0x09 ? Int(ntoh(read(s, UInt32))) :
  type_kind == 0x0B ? (Timestamp(ntoh(read(s, UInt64)))) :
  type_kind == 0x0C ? (Base.Random.UUID(ntoh(read(s, UInt128)))) :
  type_kind == 0x0D ? bytestring(readbytes(s, nrOfBytes)) :
  type_kind == 0x20 ? decodeList(s, val_type) :
  type_kind == 0x21 ? decodeDict(s, key_type, val_type) :
  type_kind == 0x22 ? decodeMap(s, val_type) :
            ("*NYI*", type_kind, readbytes(s, nrOfBytes))
end

function decodeResultRowTypes(s::IOBuffer)
  flags   = ntoh(read(s, UInt32));
  colcnt  = ntoh(read(s, UInt32));
  gl_spec = (flags & 0x0001) != 0;

  if gl_spec then
    gl_ksname = decodeString(s);
    gl_tablename = decodeString(s);
  end
  if (flags & 0x0002) != 0 then
    println("d >> ");
  end
  if (flags & 0x0004) != 0 then
    println("e >> ");
  end

  types = Array(Tuple{Int16,Any,Any}, colcnt);
  for col in 1:colcnt
    ksname     = gl_spec ? gl_ksname : decodeString(s);
    tablename  = gl_spec ? gl_tablename : decodeString(s);
    name       = decodeString(s);
    kind       = ntoh(read(s, UInt16));
    key_type   = kind == 0x21 ?  ntoh(read(s, UInt16)) : nothing;
    sub_type   = kind in UInt8[0x20, 0x21, 0x22] ? ntoh(read(s, UInt16)) : nothing;
    types[col] = (kind, sub_type, key_type);
  end
  types
end

function decodeResultColumn(s::IO, typ)
  nrOfBytes = ntoh(read(s, Int32));
  decodeValue(s, nrOfBytes, typ)
end

function decodeResultRows(s::IOBuffer)
  types = decodeResultRowTypes(s);
  nrOfColumns = length(types);
  nrOfRows = ntoh(read(s, UInt32));
  rows = Array(Array{Any}, nrOfRows)
  for rix = 1:nrOfRows
    rows[rix] = Array(Any, nrOfColumns);
    for cix in 1:nrOfColumns
      rows[rix][cix] = decodeResultColumn(s, types[cix])
    end
  end
  rows
end

function decodeResultMessage(buffer::Array{UInt8})
  s = IOBuffer(buffer)
  kind = ntoh(read(s, UInt32))
  return kind == 0x01 ? ("void") :
  kind == 0x02 ? decodeResultRows(s) :
  kind == 0x03 ? ("set keyspace", decodeString(s)) :
  kind == 0x04 ? ("prepared") :
  kind == 0x05 ? ("schema change", decodeString(s), decodeString(s), decodeString(s)) :
    ("???");
end

function decodeErrorMessage(buffer::Array{UInt8})
  s = IOBuffer(buffer);
  kind = ntoh(read(s, UInt32));
  errmsg = decodeString(s);
  println("ERROR (DECODE) [$kind]: ", errmsg);
  ("ERROR", kind, errmsg)
end

function decodeMessage(opcode::UInt8, buffer::Array{UInt8})
  opcode == 0x08 ? decodeResultMessage(buffer) :
  opcode == 0x00 ? decodeErrorMessage(buffer) :
         (opcode, buffer);
end

####################################################################
# Encoding
####################################################################

function cql_encode_string(buf::IOBuffer, str::AbstractString)
  encStr = bytestring(isvalid(UTF8String, str) ? str : utf8(str));
  write(buf, hton(UInt16(sizeof(encStr))));
  write(buf, encStr);
  nothing
end

function cql_encode_long_string(buf::IOBuffer, str::AbstractString)
  encStr = bytestring(isvalid(UTF8String, str) ? str : utf8(str));
  write(buf, hton(UInt32(sizeof(encStr))));
  write(buf, encStr);
  nothing
end

function cql_encode(buf::IOBuffer, dict::Dict)
  write(buf, hton(UInt16(length(dict))));
  for (k,v) in dict
    cql_encode_string(buf, k);
    cql_encode_string(buf, v);
  end
  nothing
end

function cql_encode(buf::IOBuffer, query::AbstractString)
  cql_encode_long_string(buf, query);
  write(buf, 0x00);
  write(buf, 0x04);
  write(buf, 0x00);
  nothing
end

####################################################################
# Sending Message to the server
####################################################################

function sendMessageBody(con::CQLConnection, msg)
  buf = con.buffer;
  truncate(buf, 0);
  cql_encode(buf, msg);
  write(con.socket, hton(UInt32(buf.size)));
  write(con.socket, takebuf_array(buf));
end

function sendMessage(con::CQLConnection, kind::UInt8, msg,
                     id::UInt16 = 0x0000)
  con.pending += 1;
  write(con.socket, 0x03);
  write(con.socket, 0x00);
  write(con.socket, id);
  write(con.socket, kind);
  sendMessageBody(con, msg);

  flush(con.socket);
  yield();
  nothing
end

function nextReplyID(con::CQLConnection)
  id::UInt8 = con.msg_id;
  con.msg_id = 1 + (id + 1) % 99;
  while haskey(con.replies, id)
    yield();
  end
  reply = RemoteRef();
  con.replies[id] = reply;
  (id, reply)
end

####################################################################
# Queries
####################################################################

function query(con::CQLConnection, msg::AbstractString)
  sync(con);
  getResult(asyncQuery(con, msg))
end

function command(con::CQLConnection, msg::AbstractString)
  sync(con);
  asyncCommand(con, msg)
  sync(con);
  nothing
end

function asyncQuery(con::CQLConnection, msg::AbstractString)
  id::UInt16, reply = nextReplyID(con);
  sendMessage(con, 0x07, msg, id);
  reply
end

function asyncCommand(con::CQLConnection, msg::AbstractString)
  sendMessage(con, 0x07, msg);
  nothing
end

function getResult(reply::RemoteRef)
  decodeMessage(take!(reply) ...)
end

function sync(con::CQLConnection)
  while 0 < con.pending
    yield();
  end
end

####################################################################
# Helpers
####################################################################

cleanString(str) = replace(str,"'","\\'")
originalString(str) = replace(str,"\\'","'")

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

end


text:
module ConstraintTests

using Test
using WeaklyHard
using WeaklyHard: consword

##############################
### Generate all sequences ###
##############################

const k_max = 5
const types = [Int8, Int16, Int32, Int64, BigInt]

generate_all_seq(k::T) where {T <: Integer} = (k <= 8*sizeof(T) - 1) ? [(T(1) << k) | i for i in T(0):T(2^k-1)] : error("Can't create that long sequences of $T")

const seqs = generate_all_seq(Int32(k_max))
@show length(seqs)

generate_all_valid(lambda::T) where {T <: Constraint} = Set{Int32}(filter(x -> is_satisfied(lambda, x), seqs))

################################################
### Aux functions for individual constraints ###
################################################

function trivial_rowhit(x, k)
    lambda = RowHitConstraint(x, k)

    # Parameter tests
    @test lambda.x == x
    @test lambda.k == k
    z = k - 2*x + 1

    # Function tests
    for T in types
        s = 8*sizeof(T)
        if k + 2*x < s
            # Trivial tests
            @test (x==0) ? is_satisfied(lambda, T(0)) : !is_satisfied(lambda, T(0))
            @test (2*x > k) ? !is_satisfied(lambda, consword(T(k)) << 1) : is_satisfied(lambda, consword(T(k)) << 1)
        end # if
    end # for
end # function

function trivial_rowmiss(x, k)
    lambda = RowMissConstraint(x, k)

    # Parameter tests
    @test lambda.x == x
    @test lambda.k == x+1

    # Function tests
    for T in types
        s = 8*sizeof(T)
        if k < s
            # Trivial tests
            if x == 0
                @test !is_satisfied(lambda, T(1) << 1)
            else
                @test is_satisfied(lambda, T(1) << lambda.x)
                @test !is_satisfied(lambda, T(1) << (lambda.x+1))
            end # if
        end # if
    end # for
end # function

function trivial_anymiss(x, k)
    lambda = AnyMissConstraint(x, k)

    # Parameter tests
    @test lambda.x == x
    @test lambda.k == k

    # Function tests
    for T in types
        s = 8*sizeof(T)
        if k < s
            # Trivial tests
            if x == 0
                @test !is_satisfied(lambda, T(1) << 1)
            elseif x == k
                @test is_satisfied(lambda, T(0))
            else
                @test is_satisfied(lambda, consword(T(lambda.k)) << lambda.x)
                @test !is_satisfied(lambda, consword(T(lambda.k)) << (lambda.x+1))
            end # if
        end #if
    end # for
end # function

function trivial_anyhit(x, k)
    lambda = AnyHitConstraint(x, k)

    # Parameter tests
    @test lambda.x == x
    @test lambda.k == k

    # Function tests
    for T in types
        s = 8*sizeof(T)
        if k < s
            # Trivial tests
            if x == 0
                @test is_satisfied(lambda, T(0))
            elseif x == k
                @test !is_satisfied(lambda, T(1) << 1)
            else
                @test is_satisfied(lambda, consword(T(lambda.k)) << (lambda.k - lambda.x))
                @test !is_satisfied(lambda, consword(T(lambda.k)) << (lambda.k - lambda.x + 1))
            end # if
        end # if
    end # for
end # function

########################################
### Tests for individual constraints ###
########################################

@testset "Trivial Tests - RowHitConstraint" begin
    for k in 1:200
        for x in 0:k
            trivial_rowhit(x, k)
        end # for
    end # for
end # testset

@testset "Trivial Tests - RowMissConstraint" begin
    for k in 1:200
        for x in 0:k
            trivial_rowmiss(x, k)
        end # for
    end # for
end # testset

@testset "Trivial Tests - AnyMissConstraint" begin
    for k in 1:200
        for x in 0:k
            trivial_anymiss(x, k)
        end # for
    end # for
end # testset

@testset "Trivial Tests - AnyHitConstraint" begin
    for k in 1:200
        for x in 0:k
            trivial_anyhit(x, k)
        end # for
    end # for
end # testset

###############################
### Aux dominance functions ###
###############################

function rowmissrowmiss_tests(x1, x2)
    c1 = RowMissConstraint(x1)
    c2 = RowMissConstraint(x2)

    if x1 == x2
        @test is_dominant(c1, c2) && is_dominant(c2, c1)
    elseif x1 < x2
        @test is_dominant(c1, c2)
    else
        @test is_dominant(c2, c1)
    end
end # function

function rowhitrowhit_tests(x1, k1, x2, k2)
    c1 = RowHitConstraint(x1, k1)
    c2 = RowHitConstraint(x2, k2)

    if x2 == 0
        @test is_dominant(c1, c2)
    elseif x1 == 0
        @test !is_dominant(c1, c2)
    elseif k2 >= k1
        @test (x2 <= x1) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
    else
        @test (x2 <= x1 - ceil( (k1-k2)/2 )) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
    end
end # function

function anymissanymiss_tests(x1, k1, x2, k2)
    c1 = AnyMissConstraint(x1, k1)
    c2 = AnyMissConstraint(x2, k2)
    n1 = k1 - x1
    n2 = k2 - x2

    @test (n2 <= max(n1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(n1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function anyhitanyhit_tests(x1, k1, x2, k2)
    c1 = AnyHitConstraint(x1, k1)
    c2 = AnyHitConstraint(x2, k2)

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

#######################
### Dominance tests ###
#######################

@testset "RowMiss-RowMiss - Tests" begin
    for i in 1:10000
        x1 = rand(0:i)
        x2 = rand(0:i)
        rowmissrowmiss_tests(x1, x2)
    end # for
end # testset

@testset "RowHit-RowHit - Tests" begin
    for i in 2:10001
        k1 = rand(1:i)
        k2 = rand(1:i)
        x1 = rand(0:floor(Int, k1/2))
        x2 = rand(0:floor(Int, k2/2))
        rowhitrowhit_tests(x1, k1, x2, k2)
    end # for
end # testset

@testset "AnyMiss-AnyMiss - Tests" begin
    for i in 1:10000
        k1 = rand(1:i)
        k2 = rand(1:i)
        x1 = rand(0:k1)
        x2 = rand(0:k2)
        anymissanymiss_tests(x1, k1, x2, k2)
    end # for
end # testset

@testset "AnyHit-AnyHit - Tests" begin
    for i in 1:10000
        k1 = rand(1:i)
        k2 = rand(1:i)
        x1 = rand(0:k1)
        x2 = rand(0:k2)
        anyhitanyhit_tests(x1, k1, x2, k2)
    end # for
end # testset

#################################
### Aux cross-dominance tests ###
#################################

### RowMiss-X ###
function rowmissrowhit_tests(x1, k1, x2, k2)
    c1 = RowMissConstraint(x1)
    c2 = RowHitConstraint(x2, k2)

    x1 = 1
    k1 = c1.x+1
    z1 = k1-x1

    @test (x2 <= min(floor(k2/(z1+1)), max(floor(k1/(z1+1)), ceil(x1/z1)))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function rowmissanymiss_tests(x1, k1, x2, k2)
    c1 = RowMissConstraint(x1)
    c2 = AnyMissConstraint(x2, k2)

    x1 = 1
    k1 = c1.x+1
    x2 = k2 - x2

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function rowmissanyhit_tests(x1, k1, x2, k2)
    c1 = RowMissConstraint(x1)
    c2 = AnyHitConstraint(x2, k2)

    x1 = 1
    k1 = c1.x+1
    x2 = x2

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

### RowHit-X ###
function rowhitrowmiss_tests(x1, k1, x2, k2)
    c1 = RowHitConstraint(x1, k1)
    c2 = RowMissConstraint(x2)

    z1 = k1-2x1+1
    x2 = 1
    k2 = c2.x+1

    @test (x2 <= max(x1*floor(k2/(z1+x1)), k2 - floor(k2/(z1+x1))*z1 - z1)) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function rowhitanymiss_tests(x1, k1, x2, k2)
    c1 = RowHitConstraint(x1, k1)
    c2 = AnyMissConstraint(x2, k2)

    z1 = k1-2x1+1
    x2 = k2-x2

    @test (x2 <= max(x1*floor(k2/(z1+x1)), k2 - floor(k2/(z1+x1))*z1 - z1)) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function rowhitanyhit_tests(x1, k1, x2, k2)
    c1 = RowHitConstraint(x1, k1)
    c2 = AnyHitConstraint(x2, k2)

    z1 = k1-2x1+1

    return x2 <= max(x1*floor(k2/(z1+x1)), k2 - floor(k2/(z1+x1))*z1 - z1)
end # function

### AnyMiss-X ###
function anymissrowmiss_tests(x1, k1, x2, k2)
    c1 = AnyMissConstraint(x1, k1)
    c2 = RowMissConstraint(x2)

    x1 = k1 - x1
    x2 = 1
    k2 = c2.x+1

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function anymissrowhit_tests(x1, k1, x2, k2)
    c1 = AnyMissConstraint(x1, k1)
    c2 = RowHitConstraint(x2, k2)

    z1 = x1
    x1 = k1-x1

    @test (x2 <= min(floor(k2/(z1+1)), max(floor(k1/(z1+1)), ceil(x1/z1)))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function anymissanyhit_tests(x1, k1, x2, k2)
    c1 = AnyMissConstraint(x1, k1)
    c2 = AnyHitConstraint(x2, k2)

    x1 = k1 - x1

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

### AnyHit-X ###
function anyhitrowmiss_tests(x1, k1, x2, k2)
    c1 = AnyHitConstraint(x1, k1)
    c2 = RowMissConstraint(x2)

    x2 = 1
    k2 = c2.x+1

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function anyhitrowhit_tests(x1, k1, x2, k2)
    c1 = AnyHitConstraint(x1, k1)
    c2 = RowHitConstraint(x2, k2)

    z1 = k1-x1

    @test (x2 <= min(floor(k2/(z1+1)), max(floor(k1/(z1+1)), ceil(x1/z1)))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

function anyhitanymiss_tests(x1, k1, x2, k2)
    c1 = AnyHitConstraint(x1, k1)
    c2 = AnyMissConstraint(x2, k2)

    x2 = k2-x2

    @test (x2 <= max(x1*floor(Int, k2/k1), k2+ceil(Int, k2/k1)*(x1-k1))) ? is_dominant(c1, c2) : !is_dominant(c1, c2)
end # function

#############################
### Cross-dominance tests ###
#############################

@testset "RowMiss-X - Analytic Tests" begin
    for i in 1:10000
        k = rand(1:i, 4)
        x = [rand(0:j) for j in k]

        rowmissrowhit_tests(x[1], x[1]+1, floor(Int, x[2]/2), k[2])
        rowmissanymiss_tests(x[1], x[1]+1, x[3], k[3])
        rowmissanyhit_tests(x[1], x[1]+1, x[4], k[4])
    end # for
end # testset

@testset "RowHit-X - Analytic Tests" begin
    for i in 1:10000
        k = rand(1:i, 4)
        x = [rand(0:j) for j in k]

        rowhitrowmiss_tests(floor(Int, x[2]/2), k[2], x[1], x[1]+1)
        rowhitanymiss_tests(floor(Int, x[2]/2), k[2], x[3], k[3])
        rowhitanyhit_tests(floor(Int, x[2]/2), k[2], x[4], k[4])
    end # for
end # testset

@testset "AnyMiss-X - Analytic Tests" begin
    for i in 1:10000
        k = rand(1:i, 4)
        x = [rand(0:j) for j in k]

        anymissrowmiss_tests(x[3], k[3], x[1], x[1]+1)
        anymissrowhit_tests(x[3], k[3], floor(Int, x[2]/2), k[2])
        anymissanyhit_tests(x[3], k[3], x[4], k[4])
    end # for
end # testset

@testset "AnyHit-X - Analytic Tests" begin
    for i in 1:10000
        k = rand(1:i, 4)
        x = [rand(0:j) for j in k]

        anyhitrowmiss_tests(x[4], k[4], x[1], x[1]+1)
        anyhitrowhit_tests(x[4], k[4], floor(Int, x[2]/2), k[2])
        anyhitanymiss_tests(x[4], k[4], x[3], k[3])
    end # for
end # testset

#######################
### satisfied tests ###
#######################

Threads.@threads for k1 = 1:k_max-1
    for x1 in 0:k1
        l1 = RowMissConstraint(x1)
        S1 = generate_all_valid(l1)

        @testset "BruteForce Tests" begin
            for k2 = 1:k_max
                for x2 = 0:k2
                    l2 = RowHitConstraint(x2, k2)
                    S2 = generate_all_valid(l2)

                    @test is_dominant(l1, l2) ? isempty(setdiff(S1, S2)) : !isempty(setdiff(S1, S2))

                    @test is_dominant(l2, l1) ? isempty(setdiff(S2, S1)) : !isempty(setdiff(S2, S1))

                    for k3 = 1:k_max
                        for x3 = 0:k3
                            l3 = AnyMissConstraint(x3, k3)
                            S3 = generate_all_valid(l3)

                            @test is_dominant(l1, l3) ? isempty(setdiff(S1, S3)) : !isempty(setdiff(S1, S3))

                            @test is_dominant(l2, l3) ? isempty(setdiff(S2, S3)) : !isempty(setdiff(S2, S3))

                            @test is_dominant(l3, l1) ? isempty(setdiff(S3, S1)) : !isempty(setdiff(S3, S1))
                            @test is_dominant(l3, l2) ? isempty(setdiff(S3, S2)) : !isempty(setdiff(S3, S2))

                            for k4 = 1:k_max
                                for x4 = 0:k4
                                    l4 = AnyHitConstraint(x4, k4)
                                    S4 = generate_all_valid(l4)

                                    @test is_dominant(l1, l4) ? isempty(setdiff(S1, S4)) : !isempty(setdiff(S1, S4))

                                    @test is_dominant(l2, l4) ? isempty(setdiff(S2, S4)) : !isempty(setdiff(S2, S4))

                                    @test is_dominant(l3, l4) ? isempty(setdiff(S3, S4)) : !isempty(setdiff(S3, S4))

                                    @test is_dominant(l4, l1) ? isempty(setdiff(S4, S1)) : !isempty(setdiff(S4, S1))
                                    @test is_dominant(l4, l2) ? isempty(setdiff(S4, S2)) : !isempty(setdiff(S4, S2))
                                    @test is_dominant(l4, l3) ? isempty(setdiff(S4, S3)) : !isempty(setdiff(S4, S3))
                                end # for
                            end # for
                        end # for
                    end # for
                end # for
            end # for
        end # testset
    end # for
end # for

end # module

ext: jl


text:
@doc """
Calculates the stationarity KPSS test statistic of Nyblom (1989) and Kwiatkowski et all (1992).

```math
 y_t = \\xi_t + r_t + \\epsilon_t
```
```math
 r_t = r_{t-1} + u_t, \\quad u_t\\sim\\mathcal{N}(0, \\sigma^2_u)
```

The KPSS test has a LM-type statistic given by

```math
  LM = \\frac{\\sum^T_t S_t}{\\epsilon'\\epsilon},
```

where Sₜ is the partial sum of resids up to time t

```math
  S_t = \\sum^t_{j=1} \\epsilon_j
```

Required Input Parameters
----------------------------------
`y::Vector{Float64}`: Observed time series in vector form.

Optional Input Parameters
----------------------------------

`trend::String`: Deterministic trend specification. It determines the asymptotic distribution of the t-statistic. Possible values are `"constant"`, which considers a constant deterministic term; amd `"linear"`, which also considers a linear time trend.

Output
----------------------------------
`kpss_stat::Float64`: KPSS test statistic.

`kpss_pval::Float64`: p-value obtained from 10000 simulations.

```@example 1
  series = simul_arma()
  kpss_test(series)
```

"""->
function kpss_test(y::Vector{Float64}; trend::String="constant")
  T = size(y,1)
  if trend=="constant"
    x = ones(T)
    println("constant deteriministic detrending.")
  elseif trend=="linear"
    x = [ones(T) collect(1:T)]
    println("constant plus linear deteriministic detrending.")
  else
    error(" trend must be either, constant, or linear")
  end
  β, t_stat, resid = kpss_ols(y, x[:,:])
  S = cumsum(resid, 1)
  nominator = sum(S.^2)/(T^2)
  denominator = resid'resid/T
  kpss_stat = nominator./denominator
  kpss_pval = mean(kpss_stat.<kpss_dist(T, trend=trend))
  return kpss_stat[1], kpss_pval
end
#
#
function wiener2(nobs)
    e0 = rand(Normal(), nobs);
    x = ones(nobs,1);
    e1 = e0 - mean(e0);
    b = inv(x'x)*(x'e0);
    e2 = e0-b.*e0;
    y1 = cumsum(e1);
    y2 = cumsum(e2);
    intW2_1 = nobs^(-2.)*sum(y1.^2, 1);
    intW2_2 = nobs^(-2.)*sum(y2.^2, 1);
    return intW2_1, intW2_2
end

function kpss_dist(T; trend="constant")
  nsim = 10000
  kpss = Vector{Float64}(nsim)
  for isim ∈ 1:nsim
    c_tmp, trnd_tmp = wiener2(T)
    (trend=="constant") && (kpss[isim, :] = c_tmp)
    (trend=="linear") && (kpss[isim, :] = trnd_tmp)
  end
  return kpss
end

function kpss_ols(y::Array{Float64,}, x::Array{Float64,})
  T, k = size(x)
  β = inv(x'x)*(x'y)
  ϵ = y - x*β
  yhat = x*β
  sse = ϵ'ϵ
  σ2_y = sse/(T-k)
  Σₓ = inv(x'x)
  Σᵦ = σ2_y.*Σₓ
  t_stat = β./sqrt.(diag(Σᵦ))
  bic = (T-k).*log(sse./(T-k)) .+ k.*log(T-k)
  resid = y - x*β
  return β, t_stat, resid
end

#kpss_test(randn(200))


text:
using Random

""" Energy of the GCM model """
function Energy(X, parameters)
    x, y, px, py = X
    A, B, C, K  = parameters

    b2 = x * x + y * y
    p2 = px * px + py * py

    return 0.5 * p2 / K + A * b2 + B * x * (x * x - 3.0 * y * y) + C * b2 * b2
end

""" Filling initial conditions (first for section (x, px) through plane y = 0)"""
function InitialCondition!(x0, e, parameters)
    numMissingValues = length(x0) - length(collect(skipmissing(x0)))

    x, y, px, py = x0
    A, B, C, K  = parameters

    if numMissingValues == 1 && ismissing(py)
        d = 2.0 * K * (e - Energy([x, y, px, 0], parameters))
        if d < 0
            return false
        end

        sign = rand() >= 0.5 ? 1.0 : -1.0
        py = sign * sqrt(d)
        x0[4] = py

        return true
    end

    return false
end

""" True if given trajectory is within the kinematically accessible domain """
function CheckDomain(x, parameters, t)
    return false
end

function EquationOfMotion!(dx, X, parameters, t)
    x, y, px, py = X[1:4]             # Equations of motion for trajectories
    Φ = reshape(X[5:end], 4, 4)     # Tangent dynamics

    A, B, C, K = parameters[1:4]

    b2 = x * x + y * y
    Hx = 2.0 * A * x + 3.0 * B * (x * x - y * y) + 4.0 * C * x * b2
    Hy = 2.0 * y * (A - 3.0 * B * x + 2.0 * C * b2)

    Hxx = 2.0 * (A + 3.0 * B * x + 2.0 * C * (3.0 * x * x + y * y))
    Hxy = y * (-6.0 * B + 8.0 * C * x)
    Hyy = 2.0 * (A - 3.0 * B * x + 2.0 * C * (x * x + 3.0 * y * y))

    dx[1] = px / K
    dx[2] = py / K
    dx[3] = -Hx
    dx[4] = -Hy

    j = [0 0 1.0/K 0;
         0 0 0 1.0/K;
         -Hxx -Hxy 0 0;
         -Hxy -Hyy 0 0]

    dx[5:end] = j * Φ
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
@testset "PWM source prior setup, PWM source initialisation and manipulation functions" begin
    #test dirichlet prior estimation from wm inputs
    wm_input = [.0 .2 .3 .5; .0 .2 .3 .5]
    est_dirichlet_vec = estimate_dirichlet_prior_on_wm(wm_input)
    @test typeof(est_dirichlet_vec) == Vector{Dirichlet{Float64}}
    for pos in 1:length(est_dirichlet_vec)
        @test isapprox(est_dirichlet_vec[pos].alpha, wm_input[pos,:].*PRIOR_WT)
    end

    bad_input = wm_input .* 2
    @test_throws DomainError estimate_dirichlet_prior_on_wm(bad_input)

    wm_input = [.1 .2 .3 .4; .1 .2 .3 .4]
    est_dirichlet_vec = estimate_dirichlet_prior_on_wm(wm_input)
    @test typeof(est_dirichlet_vec) == Vector{Dirichlet{Float64}}
    for pos in 1:length(est_dirichlet_vec)
        @test est_dirichlet_vec[pos].alpha == wm_input[pos,:].*PRIOR_WT
    end

    length_range = 2:2

    #test informative/uninformative source prior vector assembly
    test_priors = assemble_source_priors(2, [wm_input])
    @test length(test_priors)  == 2
    for pos in 1:length(test_priors[1])
        @test test_priors[1][pos].alpha == wm_input[pos,:].*PRIOR_WT
    end
    @test test_priors[2] == false

    #test source wm initialisation from priors
    test_sources = init_logPWM_sources(test_priors, length_range)
    for source in test_sources
        for pos in 1:size(source[1])[1]
            @test isprobvec(exp.(source[1][pos,:]))
        end
    end

    #test that wm_shift is returning good shifted probvecs
    rando_dist=Dirichlet([.25,.25,.25,.25])
    for i in 1:1000
        wm=rand(rando_dist)
        new_wm=wm_shift(wm,Weibull(1.5,.1))
        @test isprobvec(new_wm)
        @test wm!=new_wm
    end

    #test that legal new sources are generated by permute_source_weights
    permuted_weight_sources=deepcopy(test_sources)
    permuted_weight_sources[1]=permute_source_weights(permuted_weight_sources[1],1.,Weibull(1.5,.1))
    permuted_weight_sources[2]=permute_source_weights(permuted_weight_sources[2],1.,Weibull(1.5,.1))
    @test permuted_weight_sources != test_sources
    for (s,source) in enumerate(permuted_weight_sources)
        for pos in 1:size(source[1],1)
            @test isprobvec(exp.(source[1][pos,:]))
            @test source[1][pos,:] != test_sources[s][1][pos,:]
        end
    end

    #test that get_length_params returns legal length shifts
    lls=1:10
    pr=1:5
    for i in 1:1000
        srcl=rand(1:10)
        sign, permute_length=get_length_params(srcl, lls, pr)
        @test pr[1]<=permute_length<=pr[end]
        @test sign==-1 || sign==1
        @test lls[1]<=(srcl+(sign*permute_length))<=lls[end]
    end

    permuted_length_sources=deepcopy(test_sources)
    permuted_length_sources[1]=permute_source_length(permuted_length_sources[1],test_priors[1],1:3,1:10)
    permuted_length_sources[2]=permute_source_length(permuted_length_sources[2],test_priors[2],1:3,1:10)
    for (s, source) in enumerate(permuted_length_sources)
        @test size(source[1],1) != size(test_sources[1][1],1)
        @test 1<=size(source[1],1)<=3
    end

    info_test_wm=[1. 0. 0. 0.
    .94 .02 .02 .02
    .82 .06 .06 .06
    .7 .1 .1 .1
    .67 .11 .11 .11
    .52 .16 .16 .16
    .4 .2 .2 .2
    .25 .25 .25 .25]

    #test eroding sources by finding most informational position and cutting off when information drops below threshold
    infovec=get_pwm_info(log.(info_test_wm))
    @test infovec==[2.0, 1.5774573308022544, 1.0346297041419121, 0.6432203505529603, 0.5620360019822908, 0.2403724636586433, 0.07807190511263773, 0.0]

    erosion_test_source=(log.([.25 .25 .25 .25
                         .2 .4 .2 .2
                         .7 .1 .1 .1
                         .06 .06 .06 .82
                         .7 .1 .1 .1
                         .25 .25 .25 .25]),1)

    infovec=get_pwm_info(erosion_test_source[1])
    start_idx, end_idx = get_erosion_idxs(infovec, .25, 2:8)
    @test start_idx==3
    @test end_idx==5

    eroded_pwm,eroded_prior_idx=erode_source(erosion_test_source,2:8,.25)
    for pos in 1:size(eroded_pwm,1)
        @test isprobvec(exp.(eroded_pwm[pos,:]))
    end
    @test eroded_prior_idx==3
    @test isapprox(exp.(eroded_pwm),[.7 .1 .1 .1
    .06 .06 .06 .82
    .7 .1 .1 .1])

    #make sure revcomp_pwm is reversing pwms across both dimensions
    revcomp_test_pwm = zeros(2,4)
    revcomp_test_pwm[1,1] = 1
    revcomp_test_pwm[2,3] = 1
    log_revcomp_test_pwm = log.(revcomp_test_pwm)
    @test revcomp_pwm(log_revcomp_test_pwm) == [-Inf 0. -Inf -Inf
                                                        -Inf -Inf -Inf 0.]
end


text:
using BenchmarkTools
using LinearAlgebra

function g(x::Array{Float64, 1})
    y = zeros(length(x))
    z = Diagonal(ones(length(x)))
    q = ones(length(x))
    y .= z * x .+ q
end

function g_notyped(x)
    y = similar(x)
    z = I
    q = ones(eltype(x), length(x))
    y  .= z * x .+ q
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using ClenshawCurtisBessel
using Documenter

DocMeta.setdocmeta!(ClenshawCurtisBessel, :DocTestSetup, :(using ClenshawCurtisBessel); recursive=true)

makedocs(;
    modules=[ClenshawCurtisBessel],
    authors="Zack Li",
    repo="https://github.com/xzackli/ClenshawCurtisBessel.jl/blob/{commit}{path}#{line}",
    sitename="ClenshawCurtisBessel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xzackli.github.io/ClenshawCurtisBessel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/xzackli/ClenshawCurtisBessel.jl",
    devbranch = "main"
)


text:
```
    BitStringAddresses
Module with types and methods pertaining to bitstring addresses.
```
module BitStringAddresses

using LinearAlgebra
using StaticArrays
using Setfield
using Parameters

using Base.Cartesian

export BitString
export num_bits, num_chunks, chunk_type, chunks, chunk_bits, top_chunk_bits

export AbstractFockAddress, SingleComponentFockAddress, BoseFS, BoseFS2C, FermiFS
export CompositeFS, FermiFS2C
export BoseFSIndex, FermiFSIndex
export num_particles, num_modes, num_components
export find_occupied_mode, find_mode, occupied_modes, is_occupied, num_occupied_modes
export excitation, onr, near_uniform, OccupiedModeMap

include("bitstring.jl")
include("fockaddress.jl")
include("bosefs.jl")
include("fermifs.jl")
include("multicomponent.jl")

end


text:
using Turing, Test
import Random

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@stage_testset "dynamichmc" "dynamichmc.jl" begin
    import DynamicHMC
    Random.seed!(100)
    chn = sample(gdemo_default, DynamicNUTS(), 5000);
    check_numerical(chn, [:s, :m], [49/24, 7/6], atol=0.2)
end


text:
using NBT
using DataStructures
# Test
tree = NBT.ReadNBT("bigtest.nbt")
print(tree)


text:
#--------------------------------------------------------------------------
# one-way or two-way analysis of variance
#--------------------------------------------------------------------------
struct ANOVAReturn
    title::String
	colnms::Vector
	array::Vector
	F::Vector
	p::Vector
	bartlett::Float64 # Bartlett's test statistic (chisq with k-1 degrees of freedom)
	bartlett_dof::Int
	bartlett_p::Float64
end

function anova(df::DataFrame,cvar::Symbol,args::Symbol...)
	if length(args) == 1
		return oneway(df,cvar,args[1])
	elseif length(args) == 2
		return towway(df,cvar,args[1],args[2])
	else
		error("More than two categorical variables are not allowed")
	end
end

#--------------------------------------------------------------------------
# oneway - return an one-way ANOVA result in ANOVAReturn
#--------------------------------------------------------------------------
function oneway(df::DataFrame,dep::Symbol,cat::Symbol)

	# establish data
	ba = completecases(df[:,[dep,cat]])
	dep = df[ba,dep]
	cat = df[ba,cat]

	if length(dep) == 0
		error(string(df)," does not contain usable records")
	end

	# title
	title = "One-Way Analysis of Variance"

	# N
	N = length(dep)

	# grand mean
	grandmean = mean(dep)

	# total sum of squares
	sstotal = mapreduce(x->(x-grandmean)^2,+,dep)

	# levels and number of groups
	lev = sort(collect(Set(cat)))
	k = length(lev)

	# group means and between group sum of squares
	groupmean = zeros(Float64,length(lev))
	groupvar = zeros(Float64,length(lev))
	n = zeros(Int64,length(lev))
	ssbetween = 0.
	sswithin = 0.
	for i in 1:k
		depgroup = dep[cat .== lev[i]]
		groupmean[i] = mean(depgroup)
		groupvar[i] = std(depgroup)^2
		n[i] = length(depgroup)
		ssbetween += n[i]*(groupmean[i] - grandmean)^2
		sswithin += (n[i]-1)*groupvar[i]
	end

	dfbetween = k - 1
	dfwithin = N - k
	dftotal = N - 1
	mswithin = sswithin / dfwithin
	msbetween = ssbetween / dfbetween
	mstotal = sstotal / dftotal
	F = msbetween / mswithin

	# Bartlett's chisquare of equal variance
	pooledvar = sum([(n[i] - 1)*groupvar[i] for i in 1:k])/dfwithin
	numer = dfwithin*log(pooledvar) - sum([(n[i] - 1)*log(groupvar[i]) for i in 1:k])
	denom = 1 + (sum([1/(n[i]-1) for i in 1:k]) - 1/dfwithin)/(3*dfbetween)
	bartlett = numer/denom

	return ANOVAReturn(
		title,
		["Source","SS","dof","MS","F","P-Value"],
		[
			["Between groups","Within groups","Total"],
			[ssbetween,sswithin,sstotal],
			[dfbetween,dfwithin,dftotal],
			[msbetween,mswithin,mstotal]
		],
		[F],
		[HypothesisTests.ccdf(Distributions.FDist(dfbetween,dfwithin),F)],
		bartlett,
		k - 1,
		HypothesisTests.ccdf(Distributions.Chisq(k-1),bartlett))
end

#--------------------------------------------------------------------------
# twoway - return a two-way ANOVA result in ANOVAReturn
#--------------------------------------------------------------------------
function crossmap(cat1::AbstractVector,cat2::AbstractVector)

	n = length(cat1)
	if n != length(cat2)
		error("Lengths of both vectors must be identical")
	end

	lev1 = sort(collect(Set(cat1)))
	n1 = length(lev1)
	lev2 = sort(collect(Set(cat2)))
	n2 = length(lev2)

	k=1
	d = Dict()
	for i in lev1
		for j in lev2
			d[i,j] = k
			k += 1
		end
	end

	return [d[cat1[i],cat2[i]] for i = 1:n]
end

function twoway(df::DataFrame,depvar::Symbol,cat1::Symbol,cat2::Symbol)

	# establish data
	ba = completecases(df[[depvar,cat1,cat2]])
	dep = Vector(df[ba,depvar])
	catj = Vector(df[ba,cat1])
	catk = Vector(df[ba,cat2])

	if length(dep) == 0
		error(string(df)," does not contain usable records")
	end

	# title
	title = "Two-Way Analysis of Variance"

	# N
	N = length(dep)

	# grand mean
	ybar = mean(dep)

	# levels and number of groups
	# J
	levj = sort(collect(Set(catj)))
	j = length(levj)
	ybarj = [mean(dep[catj .== i]) for i in levj]
	nj = [sum(catj .== i) for i in levj]

	# K
	levk = sort(collect(Set(catk)))
	k = length(levk)
	ybark = [mean(dep[catk .== i]) for i in levk]
	nk = [sum(catk .== i) for i in levk]

	# J x K groups
	catjk = crossmap(catj,catk)
	levjk = sort(collect(Set(catjk)))
	jk = length(levjk)
	ybarjk = [mean(dep[catjk .== i]) for i in levjk]

	# sums of squares
	ss = Vector(Float64,6)
	ss[2] = sum([nj[i]*(ybarj[i] - ybar)^2 for i in 1:j])
	ss[3] = sum([nk[i]*(ybark[i] - ybar)^2 for i in 1:k])
	ss[5] = sum([(dep[i] - ybarjk[catjk[i]])^2 for i in 1:N])
	ss[6] = sum([(dep[i] - ybar)^2 for i in 1:N])
	ss[1] = ss[6] - ss[5]
	ss[4] = ss[1] - ss[2] - ss[3]

	# degrees of freedom
	df = Vector{Int64}(6)
	df[2] = j - 1 # factor 1
	df[3] = k - 1 # factor 2
	df[4] = (j - 1)*(k - 1)
	df[1] = df[2] + df[3] + df[4]
	df[5] = N - jk
	df[6] = N - 1

	# mean squared error
	ms = ss ./ df

	# F statistics
	F = [ms[i]./ms[5] for i in 1:4]

	# P-values
	pval = [HypothesisTests.ccdf(Distributions.FDist(df[i],df[5]),F[i]) for i = 1:4]

	return ANOVAReturn(
		title,
		["Source","Partial SS","df","MS","F","P-Value"],
		[
			["Model",string(cat1),string(cat2),string(cat1,"*",cat2),"Residual","Total"],
			ss,
			df,
			ms,
		],
		F,
		pval,
		0.,
		0,
		0.)
end
#
# function anova(mout::RegressionModel; type::String = "")
# 	if in(type,["II","III"]) == false
# 		error("`type` must be either `II` or `III`")
# 	end
#
# 	if dof_residual(mout)
# 		error("residual df = 0")
# 	end
#
# 	if deviance(mout) < sqrt(eps(Float64))
# 		error("residual sum of squares is 0 (within rounding error)")
# 	end
#
# 	if typeof(mout) <: LinearModel
# 		if type == "II"
# 			return _anova_lm_II(mout,type=type)
# 		elseif type == "III"
# 			return _anova_lm_III(mout,type=type)
# 		end
# 	end
# end
#
# function _anova_lm_3(mout; type::String = "")
#
# 	# intercept
# 	intercept = mout.mf.terms.intercept
#
# 	# coef and vcov
# 	V = vcov(mout)
# 	b = coef(mout)
#
# 	# identity matrix
# 	blen = length(b)
# 	Ip = eye(Int8,blen)
#
# 	# column names
# 	Source = string.(mout.mf.terms.terms)
#
# 	# number of terms
# 	nterms = length(Source)
#
# 	# column index assignment to variables
# 	assign = mout.mm.assign
#
#
# 	# missing coefficients
# 	notaliased = !isna(coef(mout)) # find out the value returned for missing estimates
#
# 	for term in 1:nterms
# 		subs = find(x->x == term,assign)
#
# 		# hypothesis matrix
# 		hmat = Ip[subs,notaliased]
#
# 		# value
# 		hval = hmat * b
# 		hvcov = hmat * b * hmat'
#
#
#
#
#
#
# 	end
#
# end

import Base.show
Base.show(io::IO,::MIME"text/plain",a::ANOVAReturn) = show(io,a)
function show(io::IO,a::ANOVAReturn)
	width = length.(a.colnms)
	width[1] = max(width[1],maximum(length.(a.array[1])))
	width[2] = max(7,width[2],maximum(length.(Stella.strpr.(a.array[2],4)))) # 7 is minimum width
	width[3] = max(5,width[3],maximum(length.(string.(a.array[3])))) # 5 is minimum
	width[4] = max(7,maximum(length.(Stella.strpr.(a.array[4],4)))) # 7 is minimum width
	width[5] = 7 # F
	width[6] = 7 # P-value

	# title
	println(io,a.title)

	# header row
    print(io,"|  ",lpad(a.colnms[1],width[1]))
    for i=2:length(a.colnms)
        print(io,"  |  ",lpad(a.colnms[i],width[i]))
    end
    print(io,"  |\n")

	# line
	Stella.hline(io,width.+2)

	# array output
	len = length(a.array[1])
	for i = 1:len-1
		print(io,"|  ",lpad(a.array[1][i],width[1]))
		for j=2:4
			print(io,"  |  ",lpad(j == 3 ? a.array[j][i] : strpr(a.array[j][i],4),width[j]))
		end
		if i < len-1
			print(io,"  |  ",lpad(strpr(a.F[i],5),width[5]))
			print(io,"  |  ",lpad(strpr(a.p[i],5),width[6]))
			print(io,"  |")
		else
			print(io,"  |  ",repeat(" ",width[5]),"  |  ",repeat(" ",width[6]),"  |")
		end
		print(io,"\n")
	end

	# line
	Stella.hline(io,width.+2)

	# total
	print(io,"|  ",Stella.lpad(a.array[1][len],width[1]))
	for j=2:4
		print(io,"  |  ",Stella.lpad(j == 3 ? a.array[j][len] : strpr(a.array[j][len],4),width[j]))
	end
	print(io,"  |  ",repeat(" ",width[5]),"  |  ",repeat(" ",width[6]),"  |")
	print(io,"\n\n")
end


text:
using Logging

import Cassette

export
    trace_algorithm,
    time_algorithm,
    trace_runtime

# these are the methods in this file
function toskip(nm)
    startswith(string(nm), "#module_functions") ||
    :trace_algorithm == nm ||
    :time_algorithm == nm ||
    :trace_runtime == nm
end

function module_functions(modname; private=true)
    list = Tuple{Module, Symbol}[]
    for nm in names(modname,all=private)
        toskip(nm) || typeof(eval(:($modname.$nm))) <: Function && push!(list,(modname,nm))
    end
    return list
end

function getparams(sig)
    while isa(sig, UnionAll)
        sig = sig.body
    end
    sig.parameters
end

function runtime_functions(modname; private=true)
    set = Set{Tuple{Module, Symbol}}()
    for nm in names(modname,all=private)
        toskip(nm) && continue
        evald = eval(:($modname.$nm))
        if typeof(evald) <: Function
            for m in methods(evald)
                for tp in getparams(m.sig)
                    tp = isa(tp, UnionAll) ? tp.body : isa(tp, TypeVar) ? tp.ub : tp
                    isa(tp, DataType) && tp <: Runtime && push!(set, (modname, nm))
                end
            end
        end
    end
    return set
end

filteronmethod(filter, m) = filter(string(typeof(m).name.mt.name))
filteronmethod(::Nothing, m) = true

filteronargs(filter, args...) = filter(args...)
filteronargs(::Nothing, args...) = true

function dofilter(ctx, f, args...)
    ctx.metadata === nothing || (filteronmethod(ctx.metadata.methodfilter, f)
        && filteronargs(ctx.metadata.argsfilter, args...))
end

struct Meta
    methodfilter
    argsfilter
end

function meta(fnamefilter, argfilter)
    fnamefilter === nothing && argfilter === nothing ? nothing : Meta(fnamefilter, argfilter)
end

# construct Scruff function list before we start, since all the methods
# in this file will also be in the Scruff module
fns = Tuple{Module, Symbol}[]
fns = append!(fns, module_functions(Scruff))

rtfns = Tuple{Module, Symbol}[]
rtfns = append!(rtfns, runtime_functions(Scruff))

Cassette.@context DoTracing
Cassette.@context DoTiming
Cassette.@context DoTraceRuntime

rtfns = append!(rtfns, runtime_functions(Scruff.Operators))
rtfns = append!(rtfns, runtime_functions(Scruff.Algorithms))
rtfns = append!(rtfns, runtime_functions(Scruff.SFuncs))
rtfns = append!(rtfns, runtime_functions(Scruff.Utils))
rtfns = append!(rtfns, runtime_functions(Scruff.RTUtils))
rtfns = unique(p->p[2], rtfns) # unique function name

# create Cassette prehooks for Tracing Runtime

for (m,op) in rtfns
    @eval begin
        function Cassette.prehook(ctx::DoTraceRuntime, f::typeof($m.$op), args...)
            if dofilter(ctx, f, args...)
                @info "$(f)" called_by=stacktrace()[3] args=args
            end
        end
    end
end

fns = append!(fns, module_functions(Scruff.Operators))
fns = append!(fns, module_functions(Scruff.Algorithms))
fns = append!(fns, module_functions(Scruff.SFuncs))
fns = append!(fns, module_functions(Scruff.RTUtils))
fns = unique(p->p[2], fns) # unique function name

# create Cassette prehooks for Tracing
for (m,op) in fns
    @eval begin
        function Cassette.prehook(ctx::DoTracing, f::typeof($m.$op), args...)
            if dofilter(ctx, f, args...)
                @info "$(f)" called_by=stacktrace()[3] args=args
            end
        end
    end
end

# create Cassette overdub for Timing
for (m,op) in fns
    @eval begin
        function Cassette.overdub(ctx::DoTiming, f::typeof($m.$op), args...)
            val, t, bytes, gctime, memallocs = @timed Cassette.recurse(ctx, f, args...)
            if dofilter(ctx, f, args...)
                @info "$f" time=string(t*1000," ms") bytes=bytes gctime=gctime
            end
            return val
        end
    end
end

"""
    trace_runtime(alg::Function, args...; fnamefilter=nothing, argfilter=nothing)

Wrap `Runtime` methods (called with its `args`) in a trace.  All Scruff functions
that have `Runtime` as a parameter and their calling parameter values will
be logged using the `@info` macro using Julia's logging system.

Filtering can be done on either the name of the method or the args list;
both parameters are functions that take a `x::String` or a `x...`,
respectively, as input and return `true` if the method is to be logged.

The output looks like a series of
```julia
┌ Info: get_variable
│   called_by = get_definition at runtime.jl:93 [inlined]
└   args = (Instance{Int64}(Variable{Int64,DiscreteCPD{1}}(:x5, DiscreteCPD{1}([[0.35, 0.65], [0.45, 0.55]])), 0),)
```
"""
function trace_runtime(alg::Function, args...; fnamefilter=nothing, argfilter=nothing, kwargs...)
    m = meta(fnamefilter, argfilter)
    Cassette.overdub(Scruff.DoTraceRuntime(metadata=m), ()->alg(args...; kwargs...))
end

"""
    trace_algorithm(alg::Function, args...; fnamefilter=nothing, argfilter=nothing)

Wrap an algorithm (called with its `args`) in a trace.  All Scruff functions and
their calling parameter values will be logged using the `@info` macro using
Julia's logging system.

Filtering can be done on either the name of the method or the args list;
both parameters are functions that take a `x::String` or a `x...`,
respectively, as input and return `true` if the method is to be logged.

The output looks like a series of
```julia
┌ Info: outgoing_pis
│   called_by = operate at runtime.jl:387 [inlined]
└   args = (DiscreteCPD{1}([[0.35, 0.65], [0.45, 0.55]]), [0.38, 0.6200000000000001], Any[])
```
"""
function trace_algorithm(alg::Function, args...; fnamefilter=nothing, argfilter=nothing, kwargs...)
    m = meta(fnamefilter, argfilter)
    Cassette.overdub(Scruff.DoTracing(metadata=m), ()->alg(args...; kwargs...))
end

"""
    time_algorithm(alg::Function, args...; fnamefilter=nothing, argfilter=nothing)

Wrap an algorithm (callect with its `args`) such that every Scruff function
is called with the `@timed` macro.  The output is sent to the `@info` macro
using Julia's logging system.

Filtering can be done on either the name of the method or the args list;
both parameters are functions that take a `x::String` or a `x...`,
respectively, as input and return `true` if the method is to be logged.

The output looks like a series of
```julia
┌ Info: outgoing_pis
│   time = "0.0068000000000000005 ms"
│   bytes = 144
└   gctime = 0.0
```
"""
function time_algorithm(alg::Function, args...; fnamefilter=nothing, argfilter=nothing, kwargs...)
    ctx = Cassette.disablehooks(Scruff.DoTiming(metadata=meta(fnamefilter, argfilter)))
    val, t, bytes, gctime, memallocs =
        @timed Cassette.overdub(ctx, alg, args...; kwargs...)
    @info "$alg" time=string(t*1000," ms") bytes=bytes gctime=gctime
    return val
end


text:
using PooksoftOptionsKit
using PyPlot

# create base options struct -
option_parameters = PSOptionKitPricingParameters(349.72, 0.3124, (7.0/365.0),7,360.0,0.17,0.0098);

# what prices to we want to look at?
assetPriceArray = collect(range(340.0,stop=380.0,length=20))

# compute -
optionContractPriceArray = option_contract_price(assetPriceArray, option_parameters; modelTreeType=:ternary)

# compute the payoff at expiration -
result_buyer = compute_call_option_profit_loss_at_expiration(:buy, 360.0, 1.20, 340.0, 380.0)
profit_loss_buyer = result_buyer.value

# make a plot -
plot(profit_loss_buyer[:,1], (1.0/100)*profit_loss_buyer[:,3], "b")       # buyer - payoff
plot(assetPriceArray,optionContractPriceArray.value, color="r")
plot(assetPriceArray,optionContractPriceArray.value, "o", markerfacecolor="w", mec="r")


text:
using DataValues
using Test

@testset "Broadcast" begin

dv1 = DataValue(3)
dv2 = DataValue(5)
dv3 = DataValue{Int}()

@test log.(Ref(dv1)) .+ Ref(dv2) == DataValue(log(3) + 5)
@test log.(Ref(dv1)) .+ Ref(dv2) .- Ref(dv3) == DataValue{Float64}()

@test Ref(dv1) .+ 2 == DataValue(5)
@test 2 .+ Ref(dv1) == DataValue(5)

@test Ref(dv1) .+ (1,3) == (DataValue(4), DataValue(6))
@test (1,3) .+ Ref(dv1) == (DataValue(4), DataValue(6))

end


text:
# --------------
# Initialization
# --------------
using LLLplus
using Test
using LinearAlgebra
using Documenter

# --------------
# tests with small matrices
# --------------
println("tests with small matrices...")

# Matrix from http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=442
H =[1   9   1   2;
    1   8   8   3;
    7   4   5   1;
    2   6   7   1];
B1=[2   3  -2  -4
    3  -1   2   1
    1   1   6  -4
    1   3  -1   3];
T1=[0   0   1  -1
    0   1  -1   0
    0   0   0   1
    1  -3   3  -2];
(B,T) = lll(H);
@test sum(abs.(B-B1))==0
@test sum(abs.(T-T1))==0
println("...done")

# --------------
# Run tests with random matrices
# --------------

println("\nIn all the following tests, the first time includes the "*
        "JIT compilation; \nfor the second execution the "*
        "compilation is already done and the time\n"*
        "should be faster.\n")

# Time LLL decomposition of a 1000x1000 real matrix with randn entries
N = 1000;
H = randn(N,N);
println("Testing LLL on $(N)x$(N) real matrix...")
@time (B,T,Q,R) = lll(H);
@time (B,T,Q,R) = lll(H);

# Time LLL, Seysen, VBLAST decompositions of a 10x10 complex matrix with
# randn entries
N = 10;
println("\nTesting LLL on $(N)x$(N) complex matrix...")
H = randn(N,N) + im*randn(N,N);
@time (B,T,Q,R) = lll(H);
@time (B,T,Q,R) = lll(H);
println("Testing Seysen on same $(N)x$(N) complex matrix...")
@time (B,T) = seysen(H);
@time (B,T) = seysen(H);
println("Testing VBLAST on same $(N)x$(N) complex matrix...")
@time (W,P,B) = vblast(H);
@time (W,P,B) = vblast(H);
println("Testing Brun on real part of same $(N)x$(N) matrix...")
@time (B,T) = brun(real.(H));
@time (B,T) = brun(real.(H));

# Test CVP
N=500;
println("\nTesting cvp on $(N)x$(N) BPSK system...")
H=randn(N,N);
Q,R=qr(H);
u=Int.(rand([-1,1],N));
y=H*u+rand(N)/100;
uhat=cvp(Q'*y,R,Val(false),-1,1);
errRate = sum(abs.(u-uhat))
println("Error Rate is $(errRate). It should be zero or very small.\n")

# --------------
# doctests
# --------------
if VERSION >= v"1.1.0"
    println("\nRunning doctests...")
    DocMeta.setdocmeta!(LLLplus, :DocTestSetup,
                        :(using LLLplus, LinearAlgebra,Random); recursive=true)
    doctest(LLLplus)
end


text:
module GPU

using PyCall
import ..LazyHelp

const qulacs = PyNULL()

const qulacs_class = [
    :QuantumStateGpu,
    :StateVectorGpu,
]

for class in qulacs_class
    @eval begin
        @doc LazyHelp(qulacs, nameof($(class))) struct $class
            pyobj::PyObject
            function $(class)(args...)
                try
                    class4gpu = qulacs.$(class)
                    new(class4gpu(args...))
                catch e
                    if isa(e, KeyError)
                        msg = """
                        Make sure you've installed GPU-version of `qulacs`.
                        Try with the following command:
                        ```console
                        \$ pip3 install qulacs-gpu
                        ```
                        """
                        throw(error(msg))
                    else
                        rethrow(e)
                    end
                end
            end
        end

        PyObject(t::$(class)) = t.pyobj

        function Base.propertynames(t::$(class))
            propertynames(getfield(t, :pyobj))
        end

        function Base.getproperty(t::$(class), s::Symbol)
            if s ∈ fieldnames($(class))
                return getfield(t, s)
            else
                return getproperty(getfield(t, :pyobj), s)
            end
        end

        export $(class)
    end
end

function __init__()
    copy!(qulacs, pyimport("qulacs"))
end

end


text:
# extra matchers that don't do anything exciting, but save typing

Parse(r::Regex, t::Type) = Parse(Pattern(r), t)
Parse(r::Regex, t::Type, base) = Parse(Pattern(r), t, base)
Parse(m::Matcher, t::Type) = m > s -> parse(t, s)
Parse(m::Matcher, t::Type, base) = m > s -> parse(t, s, base)
PUInt() = Parse(p"\d+", UInt)
PUInt8() = Parse(p"\d+", UInt8)
PUInt16() = Parse(p"\d+", UInt16)
PUInt32() = Parse(p"\d+", UInt32)
PUInt64() = Parse(p"\d+", UInt64)
PInt() = Parse(p"-?\d+", Int)
PInt8() = Parse(p"-?\d+", Int8)
PInt16() = Parse(p"-?\d+", Int16)
PInt32() = Parse(p"-?\d+", Int32)
PInt64() = Parse(p"-?\d+", Int64)
PFloat32() = Parse(p"-?(\d*\.?\d+|\d+\.\d*)([eE]\d+)?", Float32)
PFloat64() = Parse(p"-?(\d*\.?\d+|\d+\.\d*)([eE]\d+)?", Float64)

Word() = p"\w+"
Space() = p"\s+"

# see sugar.jl for [] syntax support

Star(m::Matcher; flatten=true) = flatten ? m[0:end] : m[0:end,:&]
Plus(m::Matcher; flatten=true) = flatten ? m[1:end] : m[1:end,:&]
Star!(m::Matcher; flatten=true) = flatten ? m[0:end,:!] : m[0:end,:&,:!]
Plus!(m::Matcher; flatten=true) = flatten ? m[1:end,:!] : m[1:end,:&,:!]

# list with separator

StarList(m::Matcher, s::Matcher) = Alt(Seq(m, Star(Seq(s, m))), Epsilon())
StarList!(m::Matcher, s::Matcher) = Alt!(Seq!(m, Star!(Seq!(s, m))), Epsilon())
PlusList(m::Matcher, s::Matcher) = Seq(m, Star(Seq(s, m)))
PlusList!(m::Matcher, s::Matcher) = Seq!(m, Star!(Seq!(s, m)))


text:
module TestAsInterval

using Test
using ArbFloats

afIvl1 = midpoint_radius(ArbFloat(1.0), ArbFloat(1.0e-8))
afIvl2 = midpoint_radius(ArbFloat(1.0), ArbFloat(1.0e-10))
afIvl3 = nextfloat(afIvl2)

@test contains(afIvl1, afIvl2)
@test !contains(afIvl2, afIvl1)
@test overlap(afIvl2, afIvl3)

@test midpoint(afIvl3) > midpoint(afIvl2)
# this will not be true unless we force the conversion to ArbMag everywhere radius is explicit
# @test radius(afIvl3) == radius(afIvl2)

end # module


text:
"""
    basemap(cmd0::String=""; kwargs...)

Plot base maps and frames.

Full option list at [`psbasemap`]($(GMTdoc)basemap.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- **A** | **polygon** :: [Type => Str | []]

    No plotting is performed. Instead, we determine the geographical coordinates of the polygon
    outline for the (possibly oblique) rectangular map domain.
    ($(GMTdoc)basemap.html#a)
- $(GMT.opt_B)
- **D** | **inset** :: [Type => Str]

    Draw a simple map insert box on the map. Requires -F.
    ($(GMTdoc)basemap.html#d)
- **F** | **box** :: [Type => Str]

    Without further options, draws a rectangular border around any map insert (D), map scale (L)
    or map rose (T)
    ($(GMTdoc)basemap.html#f)
- $(GMT.opt_Jz)
- **L** | **map_scale** :: [Type => Str]

    Draw a map scale.
    ($(GMTdoc)basemap.html#l)
- $(GMT.opt_P)
- **Td** | **rose** :: [Type => Str]

    Draws a map directional rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)basemap.html#t)
- **Tm** | **compass** :: [Type => Str]

    Draws a map magnetic rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)basemap.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_f)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function basemap(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("psbasemap", cmd0, arg1)
	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd = parse_common_opts(d, cmd, [:F :UVXY :JZ :bo :c :f :p :t :params], first)
    cmd = parse_these_opts(cmd, d, [[:A :polygon]])
    #-D[unit]xmin/xmax/ymin/ymax[r][+sfile][+t] | -D[g|j|J|n|x]refpoint+wwidth[/height][+jjustify][+odx[/dy]][+sfile][+t]
	#cmd = parse_TdTmL(d, cmd)
    cmd = parse_type_anchor(d, cmd, [[:Td :rose], [:Tm :compass], [:L :map_scale], [:D :inset]])

	return finish_PS_module(d, "psbasemap " * cmd, "", K, O, true, arg1)
end

#= ---------------------------------------------------------------------------------------------------
function parse_TdTmL(d::Dict, cmd::String)
	cmd = add_opt(cmd, "Td", d, [:Td :rose],
        (map=("g", nothing, 1), inside=("j", nothing, 1), anchor=("", arg2str, 2), width="+w", justify="+j",
         fancy="+f", labels="+l", label="+l", offset="+o"))
	cmd = add_opt(cmd, "Tm", d, [:Tm :compass],
        (map=("g", nothing, 1), inside=("j", nothing, 1), anchor=("", arg2str, 2), width="+w", dec="+d", justify="+j",
         rose_primary=("+i", add_opt_pen), rose_secondary=("+p", add_opt_pen), labels="+l", label="+l", annot="+t", offset="+o"))
	cmd = add_opt(cmd, "L", d, [:L :map_scale],
        (map=("g", nothing, 1), inside=("j", nothing, 1), norm=("n", nothing, 1), paper=("x", nothing, 1),
         anchor=("", arg2str, 2), dec="+d", scale_at_lat="+c", length="+w",
         align="+a1", justify="+j", fancy="_+f", label="+l", offset="+o", units="_+u", vertical="_+v"))
end
=#

# ---------------------------------------------------------------------------------------------------
function parse_type_anchor(d::Dict, cmd::String, symbs)
	# SYMBS: [:Td :rose] | [:Tm :compass] | [:L :map_scale] | [:D :inset] [:D :pos :position]
	# or     [[:Td :rose], [:Tm :compass], ...]
    if (isa(symbs, Array{Symbol,2}))  symbs = [symbs]  end      # So that we can always do a loop
	for k = 1:length(symbs)
		opt = add_opt("", "", d, symbs[k],
			(map=("g", nothing, 1), outside=("J", nothing, 1), inside=("j", nothing, 1), norm=("n", nothing, 1), paper=("x", nothing, 1), anchor=("", arg2str, 2), annot="+t", dec="+d", labels="+l", label="+l", length="+w", width="+w", size="+w", align="+a1", justify="+j", fancy="+f", horizontal="_+h", move_annot="+m", neon="_+mc", nan="+n", offset=("+o", arg2str), rose_primary=("+i", add_opt_pen), rose_secondary=("+p", add_opt_pen), save="+s", scale_at_lat="+c", spacing="+l", translate="_+t", triangles="+e", units="_+u", vertical="_+v"))
		if (opt != "" && opt[1] != 'j' && opt[1] != 'J' && opt[1] != 'g' && opt[1] != 'n' && opt[1] != 'x')
		#if (opt != "" && !occursin(r"[Jjgnx]", opt[1]))
			if (symbs[k][2] == :pos || symbs[k][2] == :position)  opt = 'J' * opt   # Colorbar default is outside
			else                                                  opt = 'j' * opt   # All others are inside
			end
		end
		if (opt != "")  cmd *= " -" * string(symbs[k][1]) * opt  end
	end
	return cmd
end

# ---------------------------------------------------------------------------------------------------
basemap!(cmd0::String="", arg1=nothing; first=false, kw...) = basemap(cmd0, arg1; first=first, kw...)

const psbasemap  = basemap 		# Alias
const psbasemap! = basemap!		# Alias


text:
using YaoBase, YaoArrayRegister, BitBasis
import LegibleLambdas: LegibleLambda, parse_lambda
export MathGate, mathgate, @mathgate

struct MathGate{N, F <: Union{LegibleLambda, Function}, Fv <: Function} <: PrimitiveBlock{N}
    f::F
    v::Fv
end

function MathGate{N}(f::Union{LegibleLambda, Function}; bview::Function=bint) where N
    return MathGate{N, typeof(f), typeof(bview)}(f, bview)
end

"""
    mathgate(f; nbits[, bview=BitBasis.bint])

Create a [`MathGate`](@ref) with a math function `f` and number of bits. You can
select different kinds of view which this `MathGate` will be applied on. Possible
values are [`BitBasis.bint`](@ref), [`BitBasis.bint_r`](@ref),
[`BitBasis.bfloat`](@ref), [`BitBasis.bfloat_r`](@ref).

    mathgate(f; bview=BitBasis.bint) -> f(n)

Lazy curried version of `mathgate`.

# Example

We can make a classical toffoli gate on quantum register.

```julia
julia> r = ArrayReg(bit"110")
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> function toffli(b::BitStr)
           t = @inbounds b[1] ⊻ (b[3] & b[2])
           return @inbounds bit_literal(t, b[2], b[3])
       end
toffli (generic function with 1 method)

julia> g = mathgate(toffli; nbits=3)
mathgate(toffli; nbits=3, bview=bint)

julia> apply!(r, g) == ArrayReg(bit"111")
true

```
"""
function mathgate(f; nbits::Union{Int, Nothing}=nothing, bview::Function=bint)
    if nbits === nothing
        @λ(n->matgate(f; nbits=n, bview=bview))
    else
        return MathGate{nbits}(f; bview=bview)
    end
end

"""
    @mathgate f <nbits> <bview=bint>

Create a [`MathGate`](@ref) with a math function `f` and number of bits `nbits`,
binary view `bview`. Unlike [`mathgate`](@ref), `f` will be automatically
converted to a more legible form.

# Example

```jldoctest
julia> @mathgate x->x + 0b11 nbits=4
mathgate((x -> x + 0x03); nbits=4, bview=bint)
```
"""
:(@mathgate)

macro mathgate(f, nbits, bview)
    if !(nbits.head === :(=) && nbits.args[1] === :nbits)
        return :(error("expect keyword nbits, got $nbits"))
    end

    if !(bview.head === :(=) && bview.args[1] === :bview)
        return :(error("expect keyword bview, got $bview"))
    end

    return quote
        mathgate($(parse_lambda(f)); nbits=$(esc(nbits.args[2])), bview=$(esc(bview.args[2])))
    end
end

macro mathgate(f, nbits)
    if !(nbits.head === :(=) && nbits.args[1] === :nbits)
        return :(error("expect keyword nbits, got $nbits"))
    end

    return quote
        mathgate($(parse_lambda(f)); nbits=$(esc(nbits.args[2])))
    end
end

macro mathgate(f)
    return quote
        mathgate($(parse_lambda(f)))
    end
end

mathop(m::MathGate{N, F, typeof(bint)}, b::Int) where {N, F} = callmath(m)(b)

function callmath(m::MathGate{N}) where N
    @inline _value(x::BitStr) = x.val
    @inline _value(x) = x

    return function (x::T) where T
        if hasmethod(m.f, Tuple{T})
            return btruncate(m.f(x), N)
        elseif hasmethod(m.f, Tuple{T, typeof(N)})
            return m.f(x, N)
        elseif hasmethod(m.f, Tuple{BitStr{T, N}})
            return btruncate(_value(m.f(bit(x; len=N))), N)
        elseif hasmethod(m.f, Tuple{BitStr{T, N}, typeof(N)})
            return _value(m.f(bit(x; len=N), N))
        else
            error("Invalid math function call, math operation should be either f(x) or f(x, N::Int)")
        end
    end
end

function mathop(m::MathGate{N, F, typeof(bint_r)}, b::Int) where {N, F}
    return b |> x->bint_r(x; nbits=N) |> callmath(m) |> x->bint_r(x; nbits=N)
end

function mathop(m::MathGate{N, F, typeof(bfloat)}, b::Int) where {N, F}
    return b |> x->bfloat(x; nbits=N) |> callmath(m) |> x->bint(x; nbits=N)
end

function mathop(m::MathGate{N, F, typeof(bfloat_r)}, b::Int) where {N, F}
    return b |> x->bfloat_r(x; nbits=N) |> callmath(m) |> x->bint_r(x; nbits=N)
end

function apply!(r::ArrayReg, m::MathGate{N, F}) where {N, F}
    nstate = zero(r.state)
    for b in basis(r)
        b2 = mathop(m, b)
        nstate[b2+1, :] = view(r.state, b+1, :)
    end
    r.state .= nstate
    return r
end

# TODO: use trait to correct this
mat(::Type{T}, m::MathGate) where T = Matrix(BlockMap(T, m))


text:
function filterup(f, dir0; kwargs...)
    founds = String[]
    walkup(dir0; kwargs...) do path
        f(path) && push!(founds, path)
        false
    end
    founds
end


text:
load("debug.jl")

module TestDebug
import Base.*
import Debug.*

@debug function f(n)
    x = 0  # this must be line #10
    for k=1:n
        x += k
    end
    x = x*x
    x
end

function debug_hook(line::Int, file, scope::Scope)
    print(line, ":")
    if (line > 10) print("\tx = ", scope[:x]) end
    if (line == 12) print("\tk = ", scope[:k]) end
    println()
end

f(3)

end


text:
using NetworkEpidemics
using LightGraphs
using Statistics
using Plots
using ColorSchemes

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SI}) = Metapopulation(mp.h, mp.D, SI(χ*mp.dynamics.β))

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SIS}) = Metapopulation(mp.h, mp.D, SIS(χ*mp.dynamics.β, mp.dynamics.γ))

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SIR}) = Metapopulation(mp.h, mp.D, SIR(χ*mp.dynamics.β, mp.dynamics.δ))

N = 1000
M = 10

h = complete_graph(M)
g = erdos_renyi(N, 0.3)
g = random_regular_graph(N, 50)
g = static_scale_free(N, Int(20*N), 2)
g = star_graph(N)
g = watts_strogatz(N, 100, 0.1)
g = cycle_graph(N)

is_connected(g)

k = degree(g)

χ₁ = mean(closeness_centrality(g))
χ₂ = mean(k.^2)/(mean(k)*N)
χ₃ = mean(k)/N #ne(g)/(n*(n-1)/2)
χ₄ = global_clustering_coefficient(g)

β = 0.1
γ = 0.2
D = 0.1

mp = Metapopulation(h, D, SIR(β,γ))
mp_1 = CorrectedMetapopulation(χ₁, mp)
mp_2 = CorrectedMetapopulation(χ₂, mp)
mp_3 = CorrectedMetapopulation(χ₃, mp)
mp_4 = CorrectedMetapopulation(χ₄, mp)

mpx = Metaplex(g, h, D, SIR(β,γ))

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]

x0_mp = Array{Int,2}(undef, M,3)
x0_mp[2:M,:] .= [100 0 0]
x0_mp[1, :] .= [90, 10, 0]

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = rand(1:M, N)
x0_mp = zeros(Int, M, 3)
for i in 1:N
    x0_mp[x0_μ[i], x0_i[i]] += 1
end

x0_mp

tmax = 50.0
nmax = 1000000
nsims = 200
nbins = 200

ts, output = gillespie(mp, x0_mp, tmax=tmax, nmax=nmax)
ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax)

plot(ts, output[2], label="")
plot(ts_mpx, output_mpx[2], label="")

ts_av, u_av = average(mp, x0_mp, tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)
ts_av_1, u_av_1 = average(mp_1, x0_mp, tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)
ts_av_2, u_av_2 = average(mp_2, x0_mp, tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)
ts_av_3, u_av_3 = average(mp_3, x0_mp, tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)
ts_av_4, u_av_4 = average(mp_4, x0_mp, tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)
ts_av_mpx, u_av_mpx = average(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax, nbins=nbins, nsims=nsims)

ts_mf, u_mf = meanfield(mp, x0_mp, tmax=tmax, saveat=ts_av)
ts_mf_1, u_mf_1 = meanfield(mp_1, x0_mp, tmax=tmax, saveat=ts_av_1)
ts_mf_2, u_mf_2 = meanfield(mp_2, x0_mp, tmax=tmax, saveat=ts_av_2)
ts_mf_3, u_mf_3 = meanfield(mp_3, x0_mp, tmax=tmax, saveat=ts_av_3)
ts_mf_4, u_mf_4 = meanfield(mp_4, x0_mp, tmax=tmax, saveat=ts_av_4)
ts_mf_mpx, u_mf_mpx = meanfield(mpx, [x0_i,x0_μ], tmax=tmax, saveat=ts_av_mpx)

colors = ColorSchemes.tab10;

plot(ts_av, sum(u_av[2], dims=2)/N,
    label="Metapopulation",
    xlabel="Time",
    ylabel="Fraction of infected individuals",
    linestyle=:dash,
    legend=:bottomright,
    fontfamily="Deja Vu",
    color = colors[1]
    )
    plot!(ts_av_mpx, sum(u_av_mpx[2], dims=2)/N,
    label="Metaplex",
    linestyle=:dash,
    color = colors[2]
    )
    plot!(ts_av_1, sum(u_av_1[2], dims=2)/N,
    label="Closeness Centrality",
    linestyle=:dash,
    color = colors[3]
    )
    plot!(ts_av_2, sum(u_av_2[2], dims=2)/N,
    label="Second Moment",
    linestyle=:dash,
    color = colors[4]
    )
    plot!(ts_av_3, sum(u_av_3[2], dims=2)/N,
    label="Mean Degree",
    linestyle=:dash,
    color = colors[5]
    )
    plot!(ts_av_4, sum(u_av_4[2], dims=2)/N,
    label="Clustering Coefficient",
    linestyle=:dash,
    color = colors[6]
    )
    plot!(ts_mf, sum(u_mf[2], dims=2)/N,
    #label="Meanfield (No correction)",
    label = "",
#    xlabel="Time",
#    ylabel="Fraction of infected individuals",
    #legend=:bottomright,
    color = colors[1]
    )
    plot!(ts_mf_mpx, sum(u_mf_mpx[2], dims=2)/N,
#    label="Meanfield (Metaplex)",
    label = "",
    color = colors[2]
    )
    plot!(ts_mf_1, sum(u_mf_1[2], dims=2)/N,
#    label="Meanfield (Closeness Centrality)",
    label = "",
    color = colors[3]
    )
    plot!(ts_mf_2, sum(u_mf_2[2], dims=2)/N,
#    label="Meanfield (Second Moment)",
    label = "",
    color = colors[4]
    )
    plot!(ts_mf_3, sum(u_mf_3[2], dims=2)/N,
#    label="Meanfield (Mean degree)",
    label = "",
    color = colors[5]
    )
    plot!(ts_mf_4, sum(u_mf_4[2], dims=2)/N,
#    label="Meanfield (Clustering Coefficient)",
    label = "",
    color = colors[6]
    )

cols = map(x->[x.r, x.g, x.b], colors.colors)

using MATLAB

mat"""
figure
xlabel('Time')
ylabel('Fraction of Infected Individuals')
hold on
plot($ts_av, $(sum(u_av[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[1]))
plot($ts_av_mpx, $(sum(u_av_mpx[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[2]))
plot($ts_av_1, $(sum(u_av_1[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[3]))
plot($ts_av_2, $(sum(u_av_2[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[4]))
plot($ts_av_3, $(sum(u_av_3[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[5]))
plot($ts_av_4, $(sum(u_av_4[2], dims=2)/N), 'LineStyle', '--', 'Color', $(cols[6]))

plot($ts_mf, $(sum(u_mf[2], dims=2)/N), 'Color', $(cols[1]))
plot($ts_mf_mpx, $(sum(u_mf_mpx[2], dims=2)/N), 'Color', $(cols[2]))
plot($ts_mf_1, $(sum(u_mf_1[2], dims=2)/N), 'Color', $(cols[3]))
plot($ts_mf_2, $(sum(u_mf_2[2], dims=2)/N), 'Color', $(cols[4]))
plot($ts_mf_3, $(sum(u_mf_3[2], dims=2)/N), 'Color', $(cols[5]))
plot($ts_mf_4, $(sum(u_mf_4[2], dims=2)/N), 'Color', $(cols[6]))
hold off
legend('Metapopulation', 'Metaplex', 'Closeness Centrality', 'Second Moment', 'Mean Degree', 'Clustering Coefficient')
"""

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
import Downloads, JSON

@testset "PokePastes" begin
    expected = """蕾冠王（結凍頭） @ 吃剩的東西
Ability: As One (Glastrier)
Level: 50
EVs: 252 HP / 4 Atk / 252 Spe
爽朗 Nature
- 守住
- 雪矛
- 鐵壁
- 撲擊

鳳王 @ 弱點保險
Ability: 再生力
Level: 50
EVs: 252 HP / 4 Def / 164 SpA / 4 SpD / 84 Spe
內斂 Nature
IVs: 0 Atk
- 燃盡
- 空氣斬
- 大地之力
- 原始之力

雷丘 @ 氣勢披帶
Ability: 避雷針
Level: 50
EVs: 252 HP / 4 SpD / 252 Spe
膽小 Nature
- 擊掌奇襲
- 蹭蹭臉頰
- 伏特替換
- 怪異電波

熾焰咆哮虎 @ 腰木果
Ability: 威嚇
Level: 50
EVs: 244 HP / 12 SpD / 252 Spe
爽朗 Nature
- 擊掌奇襲
- 閃焰衝鋒
- 地獄突刺
- 拋下狠話

雷電雲 (M) @ 突擊背心
Ability: 不服輸
Level: 50
EVs: 76 HP / 140 Atk / 4 SpA / 172 Spe
急躁 Nature
- 瘋狂伏特
- 飛翔
- 蠻力
- 電網

多邊獸Ⅱ @ 進化奇石
Ability: 複製
Level: 50
EVs: 244 HP / 164 Def / 100 SpD
大膽 Nature
IVs: 0 Atk
- 怪異電波
- 自我再生
- 戲法空間
- 欺詐"""
    input = open(Downloads.download("https://pokepast.es/95c4a2af6731969b/json")) do fi
        JSON.parse(fi)
    end
    team = import_team(
        input["paste"],
        title = input["title"],
        author = input["author"],
        notes = input["notes"]
    )
    output = export_team(team, language = "zht")
    @test strip(output) == expected
end


text:
instructions = readlines("day12/input.txt")

mutable struct Ship
    x::Int
    y::Int
    angle::Int
end

ship = Ship(0, 0, 0)

function step_instruction(instruction)
    global ship
    command, amount = instruction[1], parse(Int, instruction[2:end])
    if command == 'N'
        ship.y += amount
    elseif command == 'S'
        ship.y -= amount
    elseif command == 'E'
        ship.x += amount
    elseif command == 'W'
        ship.x -= amount
    elseif command == 'L'
        ship.angle = mod((ship.angle + amount), 360)
    elseif command == 'R'
        ship.angle = mod((ship.angle - amount), 360)
    elseif command == 'F'
        directions = Dict(90 => 'N', 180 => 'W', 270 => 'S', 0 => 'E')
        step_instruction(directions[ship.angle] * string(amount))
    end
end

for ins in instructions
    step_instruction(ins)
end

println(abs(ship.x)+abs(ship.y))


text:
@doc raw"""
    SymmetricMatrices{n,𝔽} <: AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding}

The [`Manifold`](@ref) $ \operatorname{Sym}(n)$ consisting of the real- or complex-valued
symmetric matrices of size $n × n$, i.e. the set

````math
\operatorname{Sym}(n) = \bigl\{p  ∈ 𝔽^{n × n}\ \big|\ p^{\mathrm{H}} = p \bigr\},
````
where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transpose,
and the field $𝔽 ∈ \{ ℝ, ℂ\}$.

Though it is slightly redundant, usually the matrices are stored as $n × n$ arrays.

Note that in this representation, the complex valued case has to have a real-valued diagonal,
which is also reflected in the [`manifold_dimension`](@ref manifold_dimension(::SymmetricMatrices{N,𝔽}) where {N,𝔽}).

# Constructor

    SymmetricMatrices(n::Int, field::AbstractNumbers=ℝ)

Generate the manifold of $n × n$ symmetric matrices.
"""
struct SymmetricMatrices{n,𝔽} <: AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding} end

function SymmetricMatrices(n::Int, field::AbstractNumbers = ℝ)
    return SymmetricMatrices{n,field}()
end

function allocation_promotion_function(
    M::SymmetricMatrices{<:Any,ℂ},
    ::typeof(get_vector),
    args::Tuple,
)
    return complex
end

@doc raw"""
    check_manifold_point(M::SymmetricMatrices{n,𝔽}, p; kwargs...)

Check whether `p` is a valid manifold point on the [`SymmetricMatrices`](@ref) `M`, i.e.
whether `p` is a symmetric matrix of size `(n,n)` with values from the corresponding
[`AbstractNumbers`](@ref) `𝔽`.

The tolerance for the symmetry of `p` can be set using `kwargs...`.
"""
function check_manifold_point(M::SymmetricMatrices{n,𝔽}, p; kwargs...) where {n,𝔽}
    mpv =
        invoke(check_manifold_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    if !isapprox(norm(p - p'), 0.0; kwargs...)
        return DomainError(
            norm(p - p'),
            "The point $(p) does not lie on $M, since it is not symmetric.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::SymmetricMatrices{n,𝔽}, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to manifold point `p` on the
[`SymmetricMatrices`](@ref) `M`, i.e. `X` has to be a symmetric matrix of size `(n,n)`
and its values have to be from the correct [`AbstractNumbers`](@ref).
The optional parameter `check_base_point` indicates, whether to call
 [`check_manifold_point`](@ref)  for `p`.
The tolerance for the symmetry of `p` and `X` can be set using `kwargs...`.
"""
function check_tangent_vector(
    M::SymmetricMatrices{n,𝔽},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {n,𝔽}
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)),typeof(p),typeof(X)},
        M,
        p,
        X;
        check_base_point = false, # already checked above
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(norm(X - X'), 0.0; kwargs...)
        return DomainError(
            norm(X - X'),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not symmetric.",
        )
    end
    return nothing
end

decorated_manifold(M::SymmetricMatrices{N,𝔽}) where {N,𝔽} = Euclidean(N, N; field = 𝔽)

function get_basis(M::SymmetricMatrices, p, B::DiagonalizingOrthonormalBasis)
    Ξ = get_basis(M, p, DefaultOrthonormalBasis()).data
    κ = zeros(real(eltype(p)), manifold_dimension(M))
    return CachedBasis(B, κ, Ξ)
end

function get_coordinates!(
    M::SymmetricMatrices{N,ℝ},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    @assert size(Y) == (dim,)
    @assert size(X) == (N, N)
    @assert dim == div(N * (N + 1), 2)
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds Y[k] = X[i, j] * scale
        k += 1
    end
    return Y
end
function get_coordinates!(
    M::SymmetricMatrices{N,ℂ},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis{ℂ},
) where {N}
    dim = manifold_dimension(M)
    @assert size(Y) == (dim,)
    @assert size(X) == (N, N)
    @assert dim == N * N
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds Y[k] = real(X[i, j]) * scale
        k += 1
        if i != j # imag zero on the diagonal
            @inbounds Y[k] = imag(X[i, j]) * scale
            k += 1
        end
    end
    return Y
end

function get_vector!(
    M::SymmetricMatrices{N,ℝ},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    @assert size(X) == (dim,)
    @assert size(Y) == (N, N)
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds Y[i, j] = X[k] * scale
        @inbounds Y[j, i] = X[k] * scale
        k += 1
    end
    return Y
end
function get_vector!(
    M::SymmetricMatrices{N,ℂ},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis{ℂ},
) where {N}
    dim = manifold_dimension(M)
    @assert size(X) == (dim,)
    @assert size(Y) == (N, N)
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds Y[i, j] = (X[k] + (i == j ? 0 : X[k + 1] * 1im)) * scale
        @inbounds Y[j, i] = Y[i, j]
        k += (i == j ? 1 : 2)
    end
    return Y
end
## unify within bases later.

@doc raw"""
    manifold_dimension(M::SymmetricMatrices{n,𝔽})

Return the dimension of the [`SymmetricMatrices`](@ref) matrix `M` over the number system
`𝔽`, i.e.

````math
\begin{aligned}
\dim \mathrm{Sym}(n,ℝ) &= \frac{n(n+1)}{2},\\
\dim \mathrm{Sym}(n,ℂ) &= 2\frac{n(n+1)}{2} - n = n^2,
\end{aligned}
````

where the last $-n$ is due to the zero imaginary part for Hermitian matrices
"""
function manifold_dimension(::SymmetricMatrices{N,𝔽}) where {N,𝔽}
    return div(N * (N + 1), 2) * real_dimension(𝔽) - (𝔽 === ℂ ? N : 0)
end

@doc raw"""
    project(M::SymmetricMatrices, p)

Projects `p` from the embedding onto the [`SymmetricMatrices`](@ref) `M`, i.e.

````math
\operatorname{proj}_{\operatorname{Sym}(n)}(p) = \frac{1}{2} \bigl( p + p^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
"""
project(::SymmetricMatrices, ::Any)

project!(M::SymmetricMatrices, q, p) = copyto!(q, (p + p') ./ 2)

@doc raw"""
    project(M::SymmetricMatrices, p, X)

Project the matrix `X` onto the tangent space at `p` on the [`SymmetricMatrices`](@ref) `M`,

````math
\operatorname{proj}_p(X) = \frac{1}{2} \bigl( X + X^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
"""
project(::SymmetricMatrices, ::Any, ::Any)

project!(M::SymmetricMatrices, Y, p, X) = (Y .= (X .+ transpose(X)) ./ 2)

function Base.show(io::IO, ::SymmetricMatrices{n,F}) where {n,F}
    return print(io, "SymmetricMatrices($(n), $(F))")
end


text:
# This file is a part of TextSearch.jl

export TextModel, VectorModel,
    TfWeighting, IdfWeighting, TpWeighting,
    FreqWeighting, BinaryLocalWeighting, BinaryGlobalWeighting,
    LocalWeighting, GlobalWeighting,
    fit, vectorize, vectorize_corpus, prune, prune_select_top

#####
##
## LocalWeighting
##
#####
"""
    LocalWeighting

Abstract type for local weighting
"""
abstract type LocalWeighting end

"""
    TfWeighting()

Term frequency weighting
"""
struct TfWeighting <: LocalWeighting end

"""
    TpWeighting()

Term probability weighting
"""
struct TpWeighting <: LocalWeighting end

"""
    FreqWeighting()

Frequency weighting
"""
struct FreqWeighting <: LocalWeighting end

"""
    BinaryLocalWeighting()

The weight is 1 for known tokens, 0 for out of vocabulary tokens
"""
struct BinaryLocalWeighting <: LocalWeighting end

#####
##
## GlobalWeighting
##
#####
"""
    GlobalWeighting

Abstract type for global weighting
"""
abstract type GlobalWeighting end

"""
    IdfWeighting()

Inverse document frequency weighting
"""
struct IdfWeighting <: GlobalWeighting end

"""
    BinaryGlobalWeighting()

The weight is 1 for known tokens, 0 for out of vocabulary tokens
"""
struct BinaryGlobalWeighting <: GlobalWeighting end

#####
##
## TextModels
##
#####
"""
    Model

An abstract type that represents a weighting model
"""
abstract type TextModel end

"""
    TokenStats(occs, ndocs)

Stores useful information for a token; i.e., the number of occurrences in the collection, the number of documents having that token
"""
struct TokenStats
    occs::Int32
    ndocs::Int32
    weight::Float32
end

const UnknownTokenStats = TokenStats(0, 0, 0f0)

const Vocabulary = Dict{UInt64, TokenStats}

mutable struct VectorModel{_G<:GlobalWeighting, _L<:LocalWeighting} <: TextModel
    global_weighting::_G
    local_weighting::_L
    tokens::Vocabulary
    maxfreq::Int
    m::Int  # vocabulary size
    n::Int  # training collection size
end

function Base.copy(
        e::VectorModel;
        local_weighting=e.local_weighting,
        global_weighting=e.global_weighting,
        tokens=e.tokens,
        maxfreq=e.maxfreq,
        m=e.m,
        n=e.n
    )
    VectorModel(global_weighting, local_weighting, tokens, maxfreq, m, n)
end

function create_vocabulary(corpus)
    V = Vocabulary()
    for tokenlist in corpus
        for (token, occs) in tokenlist
            s = get(V, token, UnknownTokenStats)
            V[token] = TokenStats(s.occs + occs, s.ndocs + 1, 0f0)
        end
    end

    V
end

"""
    VectorModel(global_weighting::GlobalWeighting, local_weighting::LocalWeighting, corpus; mindocs=1)

Creates a vector model using the input corpus.
"""
function VectorModel(global_weighting::GlobalWeighting, local_weighting::LocalWeighting, corpus; mindocs=1)
    tokens = Vocabulary()

    tokens = Vocabulary()
    maxfreq = 0
    for (t, s) in create_vocabulary(corpus)
        s.ndocs < mindocs && continue
        tokens[t] = TokenStats(s.occs, s.ndocs, 0.0f0)
        maxfreq = max(maxfreq, s.occs)
    end

    VectorModel(global_weighting, local_weighting, tokens, maxfreq, length(tokens), length(corpus))
end

Base.show(io::IO, model::VectorModel) = print(io, "{VectorModel global_weighting=$(model.global_weighting), local_weighting=$(model.local_weighting), train-voc=$(model.m), train-n=$(model.n), maxfreq=$(model.maxfreq)}")

"""
    prune(model::VectorModel, lowerweight::AbstractFloat)

Creates a new vector model without terms with smaller global weights than `lowerweight`.
"""
function prune(model::VectorModel, lowerweight::AbstractFloat)
    tokens = Vocabulary()
    for (token, s) in model.tokens
        if prune_global_weighting(model, s) >= lowerweight
            tokens[token] = TokenStats(s.occs, s.ndocs, s.weight)
        end
    end

    VectorModel(model.global_weighting, model.local_weighting, tokens, model.maxfreq, model.m, model.n)
end

"""
    prune_select_top(model::VectorModel, k::Integer)
    prune_select_top(model::VectorModel, ratio::AbstractFloat)

Creates a new vector model with the best `k` tokens from `model` based on global weighting.
`ratio` is a floating point between 0 and 1 indicating the ratio of the vocabulary to be kept
"""
function prune_select_top(model::VectorModel, k::Integer)
    w = [prune_global_weighting(model, s) for (token, s) in model.tokens]
    sort!(w, rev=true)
    prune(model, Float64(w[k]))
end

prune_select_top(model::VectorModel, ratio::AbstractFloat) = prune_select_top(model, floor(Int, length(model.tokens) * ratio))

"""
    vectorize(model::VectorModel, bow::BOW, maxfreq::Integer=maximum(bow); normalize=true) where Tv<:Real

Computes a weighted vector using the given bag of words and the specified weighting scheme.
"""
function vectorize(model::VectorModel{_G, _L}, bow::BOW; normalize=true) where {_L,_G}
    numtokens = 0.0
    if _L === TpWeighting
        for (token, freq) in bow
            numtokens += freq
        end
    end

    maxfreq = (_L === TfWeighting) ? maximum(bow) : 0.0
    vec = SVEC()

    for (token, freq) in bow
        s = get(model.tokens, token, nothing)
        s === nothing && continue

        w = local_weighting(model.local_weighting, freq, maxfreq, numtokens) * global_weighting(model, s)
        if w > 1e-9
            vec[token] = w
        end
    end

    if length(vec) == 0
        # error("empty final vector: the original vector has $(length(bow)) features")
        vec[0] = 1e-9
    end

    normalize && normalize!(vec)
    vec
end

function vectorize(tok::Tokenizer, model::VectorModel, text; bow=BOW(), normalize=true)
    compute_bow(tok, text, bow)
    v = vectorize(model, bow; normalize=normalize)
    # if length(v) == 1 && haskey(v, 0)
    #     @info "empty vector"
    # end
    v
end

function vectorize_corpus(tok::Tokenizer, model::VectorModel, corpus; bow=BOW(), normalize=true)
    V = Vector{SVEC}(undef, length(corpus))

    for (i, text) in enumerate(corpus)
        empty!(bow)
        V[i] = vectorize(tok, model, text; bow=bow, normalize=normalize)
    end

    V
end

function broadcastable(model::VectorModel)
    (model,)
end

# local weightings: TfWeighting, TpWeighting, FreqWeighting, BinaryLocalWeighting
# global weightings: IdfWeighting, BinaryGlobalWeighting

local_weighting(::TfWeighting, occs, maxfreq, numtokens) = occs / maxfreq
local_weighting(::FreqWeighting, occs, maxfreq, numtokens) = occs
local_weighting(::TpWeighting, occs, maxfreq, numtokens) = occs / numtokens
local_weighting(::BinaryLocalWeighting, occs, maxfreq, numtokens) = 1.0
global_weighting(m::VectorModel{IdfWeighting}, s::TokenStats) = log2(1 + m.n / s.ndocs)
global_weighting(m::VectorModel{BinaryGlobalWeighting}, s::TokenStats) = 1.0
prune_global_weighting(m::VectorModel, s) = global_weighting(m, s)
prune_global_weighting(m::VectorModel{BinaryGlobalWeighting}, s) = -s.ndocs


text:
using Flux
using Flux: onehotbatch, onecold, flatten, Zygote
using Flux.Losses: logitcrossentropy
using Flux.Data: DataLoader
using Parameters: @with_kw
using Statistics: mean
using CUDA
using ProgressMeter: @showprogress
using Formatting

using DeepUncertainty
include("utils.jl")
include("models/bayesian_nn.jl")

if CUDA.has_cuda()
    @info "CUDA is on"
    CUDA.allowscalar(false)
end

@with_kw mutable struct Args
    batchsize::Int = 64
    lr::Float64 = 0.1
    epochs::Int = 50
    valsplit::Float64 = 0.1
    sample_size = 10
    complexity_constant = 1e-8
end

function accuracy(preds, labels)
    acc = sum(onecold(preds |> cpu) .== onecold(labels |> cpu))
    return acc
end

function test(args, loader, model)
    loss = 0
    acc = 0
    ece = 0
    entropy = 0
    ntot = 0
    for (x, y) in loader
        predictions = []
        x, y = x |> gpu, y |> gpu
        logits = model(x)

        # Loop through each model's predictions
        for ensemble = 1:args.sample_size
            logits = model(x)
            push!(predictions, logits)
        end

        # Get the mean predictions
        mean_preds = Flux.batch(predictions)
        mean_predictions = mean(mean_preds, dims = ndims(mean_preds))
        logits = dropdims(mean_predictions, dims = ndims(mean_predictions))

        n = size(logits)[end]
        loss += logitcrossentropy(logits, y) * n
        acc += accuracy(logits, y)
        ece += expected_calibration_error(logits |> cpu, onecold(y |> cpu)) * n
        entropy += mean(calculate_entropy(logits |> cpu)) * n
        ntot += n
    end

    mean_loss = loss / ntot |> round4
    mean_acc = acc / ntot * 100 |> round4
    mean_ece = ece / ntot |> round4
    mean_entropy = entropy / ntot |> round4

    @info (format(
        "Mean Loss: {} Mean Accuracy: {} Mean ECE: {} Mean Entropy: {}",
        mean_loss,
        mean_acc,
        mean_ece,
        mean_entropy,
    ))
    @info "==========================================================="
    return nothing
end

function train(; kws...)
    # Initialize the hyperparameters
    args = Args(; kws...)

    # Load the train, validation data
    train_loader, test_loader = get_data(args)

    @info("Constructing Model")
    m = VariationalResNet18(nclasses = 10) |> gpu

    ## Training
    # Defining the optimizer
    opt = Nesterov(args.lr)
    ps = Flux.params(m)

    test(args, test_loader, m)

    @info("Training....")
    # Starting to train models
    for epoch = 1:args.epochs
        @info "Epoch $epoch"

        function loss_fn(x, y)
            l = logitcrossentropy(m(x), y)
            layers = Zygote.@ignore Flux.modules(m)
            l += args.complexity_constant * sum(normal_kl_divergence.(layers))
        end

        @showprogress for (x, y) in train_loader
            x, y = x |> gpu, y |> gpu
            total_loss = 0
            gs = Flux.gradient(ps) do
                for sample in args.sample_size
                    total_loss += loss_fn(x, y)
                end
                total_loss /= args.sample_size
            end
            Flux.update!(opt, ps, gs)
        end
        test(args, test_loader, m)
    end
    test(args, test_loader, m)

    return m
end

m = train()


text:
using SelfPropelledModel
using Test
using PyPlot
include(dirname(@__FILE__)*"/mplstyle.jl")

N=160
sv=sample_particles(N,()->500.0.*(rand(2).-0.5),()->(rand(2).-0.5))

frames=[]
for i in 1:170
    X=foldl(hcat,[[p.x...] for p in sv])
    push!(frames,X)
    spm_step!(sv,0.4,0.2,0.1,1.5)
end

pts=PermutedDimsArray(cat(frames...,dims=3),(2,3,1))|>collect
pygui(false)
figure(figsize=(5,5))
for i in 1:N
    plot(pts[i,:,1],pts[i,:,2],color="black",lw=0.5)
end
xlim((-250,250))
ylim((-250,250))
tight_layout()
title("N=$(N) particles")
savefig("N$(N).pdf")


text:
"Create an empty blueprint for `n` copies of a chromosome of `len` base pairs."
plan_chrom(len::Int, n::Int) = ChromosomeBlueprint(len, n)

function plan_chrom(; len::Int = 500, ploidy::Int = 1, ngenome::Int = 1, het = 0.0, div::Float64 = 0.0)
    bp = ChromosomeBlueprint(len, ploidy * ngenome)
    @info "- ✔ Created blank blueprint"
    if div > 0.0
        if ngenome > 1
            nsites = Int(ceil(len * div))
            bp = plan_het(bp, nsites, repeat(1:ngenome, inner = ploidy))
            @info "- ✔ Planned genome diverging polymorphisms"
        else
            @error "- ✘ Could not plan genome diverging polymorphisms: only one genome is planned"
        end
    end

    bp = _try_plan_het(bp, ploidy, ngenome, het)

    return bp
end

function _try_plan_het(bp::ChromosomeBlueprint, ploidy::Int, ngenome::Int, het::Float64)
    @info "A single value $het has been provided to plan heterozygosity"
    if het > 0.0 && ngenome === 1 && ploidy === 2
        nsites = Int(floor(length(bp) * het))
        bp = plan_het(bp, nsites, 2)
        @info "- ✔ Planned heterozygosity"
    else
        @error """- ✘ Could not plan heterozygosity:
                      Don't know how to plan heterozygosity when:
                      Ploidy is $ploidy
                      The number of genomes is $ngenome
                      het has been set as $het"""
    end
    return bp
end


text:
"""
    default_options()

Sets McCormickParameters to default.
* `env_max_int = 100`: number of iterations used for envelope calculations
* `env_tol = 1E-10`: tolerance used for envelope calculations
* `mu = 0`: nonsmooth McCormick relaxations used
* `valid_check = true`: relaxtions and intervals are checked for validity
* `subgrad_refine = false`: don't use interval refinement by subgradient propagation
* `multivar_refine = false`: don't use multivariant relaxations
* `mv_tol = 1E-15`: tolerance for multivariant relaxations.
* `outer_rnding = false`: outer rounding of interval and relaxation disabled
* `outer_param = 0.0`: amount of outer rounding
"""
function default_options()
  MC_param = McCormickParamters()
  println("Options set to default values")
end

"""
    set_iterations(val::Integer)

Set iterations for used in envelope calculations.
"""
function set_iterations(val::Integer)
  MC_param.env_max_int = val
  println("Maximum number of iteration for root-finding algorithms used to determing convex/concave envelopes set to $(MC_param.env_max_int)")
end

"""
    set_tolerance(val::Float64)

Set tolerance for used in envelope calculations.
"""
function set_tolerance(val::Float64)
  MC_param.env_tol = val
end

"""
    set_diff_relax(val::Integer)

Set differentiability of relaxations used.
"""
function set_diff_relax(val::Integer)
  diff_relax = val>0
  if (diff_relax>0)
    MC_param.mu = val+1
  else
    MC_param.mu = 0
  end
end

"""
    set_validated(val)

Sets flag for using validity check for McCormick relaxations (cv<=cc).
"""
function set_validated(val)
  MC_param.valid_check = valid_check
end

"""
    set_subgrad_refine(val)

Set flag for using subgradient refinement of interval bounds.
"""
function set_subgrad_refine(val)
  MC_param.subgrad_refine = val
end

"""
    set_multivar_refine(bool,tol)

Set flag for using using multivariant MC relaxations and their tolerance.
"""
function set_multivar_refine(bool,tol)
  MC_param.multivar_refine = bool
  MC_param.mv_tol = tol
end

"""
    set_outer_rnd(bool,tol)

Set flag `bool` for outer round (`cc`/`cv`) and do so by `tol`.
"""
function set_outer_rnd(bool,tol)
  MC_param.outer_rnding = bool
  MC_param.outer_param = tol
end


text:
module AutomotiveDrivingModels

using Reexport
using Parameters
using Printf
using LinearAlgebra
using SparseArrays

@reexport using Records
@reexport using Distributions
@reexport using Vec

include("interface/main.jl") # core interface
include("1d/main.jl") # 1D interface
include("2d/main.jl") # 2D interface

end # module

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
export  EmptyQuad,
        Quad

struct Quad{T, V<:AbstractVector{<:T}} <: AbstractQuad{T}
    name::String
    Nquad::Int # number of qudrature points
    nodes::V
    weights::V

    function Quad(name::String, N::Int, nodes, weights)
        N <= 0 && throw(DomainError(N,"number of qudrature points has to be positive"))
        !(length(nodes) == length(weights)) && throw(InconsistencyError("inconsistent numbers of nodes and weights"))
        new{promote_type(eltype(nodes), eltype(weights)), promote_type(typeof(nodes), typeof(weights))}(lowercase(name), N, nodes, weights)
    end
end

struct EmptyQuad{T} <: AbstractQuad{T}
    EmptyQuad() = new{Float64}()
end

# general constructor
function Quad(N::Int, α::AbstractVector{<:Real}, β::AbstractVector{<:Real})
    !(length(α) == length(β)) && throw(InconsistencyError("inconsistent numbers of recurrence coefficients"))
    !(N <= length(α) - 1) && throw(DomainError(N,"requested number of quadrature points $N cannot be provided with $(length(α)) recurrence coefficients"))
    nodes, weights = gauss(N,α,β)
    Quad("golubwelsch", N, nodes, weights)
end

# quadrature rules for orthoPolys
Quad(N::Int, op::AbstractOrthoPoly) = Quad(N, op.α, op.β)
Quad(op::AbstractOrthoPoly) = Quad(op.deg, op)

#####################################################
#####################################################
#####################################################
# the constructor below is probably not relevant
#####################################################
#####################################################
#####################################################
# function Quad(N::Int, weight::Function, α::Vector{<:Real}, β::Vector{<:Real}, supp::Tuple{<:Real,<:Real}, symm::Bool, d::Dict=Dict())
#     m = Measure("fun_"*String(nameof(weight)),weight,supp,symm,d)
#     Quad(N,α,β,m)
# end
#####################################################
#####################################################
#####################################################

# all-purpose constructor (last resort!)
function Quad(N::Int, w::Function, dom::Tuple{<:Real,<:Real}; quadrature::Function=clenshaw_curtis)
    N <= 0 && throw(DomainError(N, "number of quadrature points has to be positive"))
    nodes, weights = quadgp(w, dom[1], dom[2], N; quadrature=quadrature)
    Quad("quadgp", N, nodes, weights)
end

function Quad(N::Int, measure::AbstractMeasure; quadrature::Function=clenshaw_curtis)
    typeof(measure) != Measure && @warn "For measures of type $(typeof(measure)) the quadrature rule should be based on the recurrence coefficients."
    Quad(N, measure.w, (measure.dom[1], measure.dom[2]); quadrature=quadrature)
end


text:
# Accessing a.L directly might involve an extra copy();
# instead, always use the stored Cholesky factor:
chol_lower(a::Cholesky) = a.uplo === 'L' ? a.L : a.U'
chol_upper(a::Cholesky) = a.uplo === 'U' ? a.U : a.L'

# For a dense Matrix, the following allows us to avoid the Adjoint wrapper:
chol_lower(a::Matrix) = cholesky(Symmetric(a, :L)).L
# NOTE: Formally, the line above should use Hermitian() instead of Symmetric(),
# but this currently has an AutoDiff issue in Zygote.jl, and PDMat is
# type-restricted to be Real, so they are equivalent.

if HAVE_CHOLMOD
    CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf.L
end