text:
using IndexableBitVectors

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

srand(12345)

function test_access{T}(::Type{T})
    b = convert(T, Bool[])
    @test_throws BoundsError b[0]
    @test_throws BoundsError b[1]

    b = convert(T, [false])
    @test_throws BoundsError b[0]
    @test b[1] == 0
    @test_throws BoundsError b[2]

    b = convert(T, [true])
    @test_throws BoundsError b[0]
    @test b[1] == 1
    @test_throws BoundsError b[2]

    b = convert(T, [false, false, true, true])
    @test_throws BoundsError b[0]
    @test b[1] == 0
    @test b[2] == 0
    @test b[3] == 1
    @test b[4] == 1
    @test_throws BoundsError b[5]

    b = convert(T, falses(1024))
    @test b[1] == false
    @test b[end] == false
    @test all([b[i] for i in 1:1024]) == false

    b = convert(T, trues(1024))
    @test b[1] == true
    @test b[end] == true
    @test all([b[i] for i in 1:1024]) == true

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @test b[i] == bitv[i]; end

    b = convert(T, rand(10) .> 0.5)
    for i in 1:10; @test isa(b[i], Bool); end
end

function test_rank{T}(::Type{T})
    b = convert(T, Bool[])
    @test rank0(b, 0) == 0
    @test rank0(b, 1) == 0
    @test rank1(b, 0) == 0
    @test rank1(b, 1) == 0

    b = convert(T, [false])
    @test rank0(b, 0) == 0
    @test rank0(b, 1) == 1
    @test rank0(b, 2) == 1
    @test rank1(b, 0) == 0
    @test rank1(b, 1) == 0
    @test rank1(b, 2) == 0

    b = convert(T, [true])
    @test rank0(b, 0) == 0
    @test rank0(b, 1) == 0
    @test rank0(b, 2) == 0
    @test rank1(b, 0) == 0
    @test rank1(b, 1) == 1
    @test rank1(b, 2) == 1

    b = convert(T, [false, false, true, true])
    # rank0
    @test rank0(b, 0) == 0
    @test rank0(b, 1) == 1
    @test rank0(b, 2) == 2
    @test rank0(b, 3) == 2
    @test rank0(b, 4) == 2
    @test rank(0, b, 0) == 0
    @test rank(0, b, 4) == 2
    @test rank(0, b, 5) == 2
    # rank1
    @test rank1(b, 0) == 0
    @test rank1(b, 1) == 0
    @test rank1(b, 2) == 0
    @test rank1(b, 3) == 1
    @test rank1(b, 4) == 2
    @test rank(1, b, 0) == 0
    @test rank(1, b, 4) == 2
    @test rank(1, b, 5) == 2

    b = convert(T, falses(1024))
    for i in 1:1024; @test rank0(b, i) == i; end
    for i in 1:1024; @test rank1(b, i) == 0; end

    b = convert(T, trues(1024))
    for i in 1:1024; @test rank0(b, i) == 0; end
    for i in 1:1024; @test rank1(b, i) == i; end

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @test rank0(b, i) == rank0(bitv, i); end
    for i in 1:1024; @test rank1(b, i) == rank1(bitv, i); end

    b = convert(T, rand(10) .> 0.5)
    for i in 0:11; @test isa(rank0(b, i), Int); end
    for i in 0:11; @test isa(rank1(b, i), Int); end
end

function test_select{T}(::Type{T})
    b = convert(T, Bool[])
    @test select0(b, -1) == 0
    @test select0(b, 0) == 0
    @test select0(b, 1) == 0
    @test select1(b, -1) == 0
    @test select1(b, 0) == 0
    @test select1(b, 1) == 0

    b = convert(T, [false])
    @test select0(b, -1) == 0
    @test select0(b, 0) == 0
    @test select0(b, 1) == 1
    @test select0(b, 2) == 0
    @test select1(b, -1) == 0
    @test select1(b, 0) == 0
    @test select1(b, 1) == 0
    @test select1(b, 2) == 0

    b = convert(T, [true])
    @test select0(b, -1) == 0
    @test select0(b, 0) == 0
    @test select0(b, 1) == 0
    @test select0(b, 2) == 0
    @test select1(b, -1) == 0
    @test select1(b, 0) == 0
    @test select1(b, 1) == 1
    @test select1(b, 2) == 0

    b = convert(T, [false, false, true, true])
    # select0
    @test select0(b, -1) == 0
    @test select0(b, 0) == 0
    @test select0(b, 1) == 1
    @test select0(b, 2) == 2
    @test select0(b, 3) == 0
    # select1
    @test select1(b, -1) == 0
    @test select1(b, 0) == 0
    @test select1(b, 1) == 3
    @test select1(b, 2) == 4
    @test select1(b, 3) == 0

    b = convert(T, falses(1024))
    for i in 0:1024; @test select0(b, i) == i; end
    @test select0(b, 1025) == 0
    for i in 0:1025; @test select1(b, i) == 0; end

    b = convert(T, trues(1024))
    for i in 0:1025; @test select0(b, i) == 0; end
    for i in 0:1024; @test select1(b, i) == i; end
    @test select1(b, 1025) == 0

    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024; @test select0(b, i) == select0(bitv, i); end
    for i in 1:1024; @test select1(b, i) == select1(bitv, i); end

    b = convert(T, rand(10) .> 0.5)
    for i in 0:11; @test isa(select0(b, i), Int); end
    for i in 0:11; @test isa(select1(b, i), Int); end
end

function test_search{T}(::Type{T})
    b = convert(T, Bool[])
    @test search0(b, 0) == 0
    @test search0(b, 1) == 0
    @test search1(b, 0) == 0
    @test search1(b, 1) == 0

    b = convert(T, [false])
    @test search0(b, 0) == 1
    @test search0(b, 1) == 1
    @test search0(b, 2) == 0
    @test search1(b, 0) == 0
    @test search1(b, 1) == 0
    @test search1(b, 2) == 0

    b = convert(T, [true])
    @test search0(b, 0) == 0
    @test search0(b, 1) == 0
    @test search0(b, 2) == 0
    @test search1(b, 0) == 1
    @test search1(b, 1) == 1
    @test search1(b, 2) == 0

    b = convert(T, [false, false, true, true])
    # search0
    @test search0(b, 0) == 1
    @test search0(b, 1) == 1
    @test search0(b, 2) == 2
    @test search0(b, 3) == 0
    @test search0(b, 4) == 0
    # search1
    @test search1(b, 0) == 3
    @test search1(b, 1) == 3
    @test search1(b, 2) == 3
    @test search1(b, 3) == 3
    @test search1(b, 4) == 4

    b = convert(T, falses(1024))
    for i in 1:1024; @test search0(b, i) == i; end
    for i in 1:1024; @test search1(b, i) == 0; end

    b = convert(T, trues(1024))
    for i in 1:1024; @test search0(b, i) == 0; end
    for i in 1:1024; @test search1(b, i) == i; end

    function linsearch(x, bv, i)
        while i ≤ endof(bv) && bv[i] != x
            i += 1
        end
        return i > endof(bv) ? 0 : i
    end
    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024
        @test search0(b, i) == linsearch(0, bitv, i)
        @test search1(b, i) == linsearch(1, bitv, i)
    end

    b = convert(T, rand(10) .> 0.5)
    for i in 0:11; @test isa(search0(b, i), Int); end
    for i in 0:11; @test isa(search1(b, i), Int); end
end

function test_rsearch{T}(::Type{T})
    b = convert(T, Bool[])
    @test rsearch0(b, 0) == 0
    @test rsearch0(b, 1) == 0
    @test rsearch1(b, 0) == 0
    @test rsearch1(b, 1) == 0

    b = convert(T, [false])
    @test rsearch0(b, 0) == 0
    @test rsearch0(b, 1) == 1
    @test rsearch0(b, 2) == 1
    @test rsearch1(b, 0) == 0
    @test rsearch1(b, 1) == 0
    @test rsearch1(b, 2) == 0

    b = convert(T, [true])
    @test rsearch0(b, 0) == 0
    @test rsearch0(b, 1) == 0
    @test rsearch0(b, 2) == 0
    @test rsearch1(b, 0) == 0
    @test rsearch1(b, 1) == 1
    @test rsearch1(b, 2) == 1

    b = convert(T, [false, false, true, true])
    # rsearch0
    @test rsearch0(b, 0) == 0
    @test rsearch0(b, 1) == 1
    @test rsearch0(b, 2) == 2
    @test rsearch0(b, 3) == 2
    @test rsearch0(b, 4) == 2
    # rsearch1
    @test rsearch1(b, 0) == 0
    @test rsearch1(b, 1) == 0
    @test rsearch1(b, 2) == 0
    @test rsearch1(b, 3) == 3
    @test rsearch1(b, 4) == 4

    b = convert(T, falses(1024))
    for i in 1:1024; @test rsearch0(b, i) == i; end
    for i in 1:1024; @test rsearch1(b, i) == 0; end

    b = convert(T, trues(1024))
    for i in 1:1024; @test rsearch0(b, i) == 0; end
    for i in 1:1024; @test rsearch1(b, i) == i; end

    function linrsearch(x, bv, i)
        while i ≥ 1 && bv[i] != x
            i -= 1
        end
        return i < 1 ? 0 : i
    end
    bitv = rand(1024) .> 0.5
    b = convert(T, bitv)
    for i in 1:1024
        @test rsearch0(b, i) == linrsearch(0, bitv, i)
        @test rsearch1(b, i) == linrsearch(1, bitv, i)
    end

    b = convert(T, rand(10) .> 0.5)
    for i in 0:11; @test isa(rsearch0(b, i), Int); end
    for i in 0:11; @test isa(rsearch1(b, i), Int); end
end

function test_long{T}(::Type{T})
    len = 2^33 + 1000
    b  = bitrand(len)
    while sum(b) ≤ typemax(UInt32)
        b  = bitrand(len)
    end
    b′ = T(b)
    for i in [1, 2, 2^32-1, 2^32, 2^32+1, len-1, len]
        @test b′[i] == b[i]
        @test rank1(b′, i) == sum(b[1:i])
    end
end

function test_copy{T}(::Type{T})
    b = convert(T, Bool[])
    @test copy(b) == b
    @test copy(b) !== b

    b = convert(T, rand(100) .> 0.5)
    @test copy(b) == b
    @test copy(b) !== b

    b = convert(T, rand(10000) .> 0.5)
    @test copy(b) == b
    @test copy(b) !== b
end

@testset "BitVector" begin
    @testset "access" begin
        test_access(BitVector)
    end
    @testset "rank" begin
        test_rank(BitVector)
    end
    @testset "select" begin
        test_select(BitVector)
    end
end

@testset "SucVector" begin
    @testset "access" begin
        test_access(SucVector)
    end
    @testset "rank" begin
        test_rank(SucVector)
    end
    @testset "select" begin
        test_select(SucVector)
    end
    @testset "search" begin
        test_search(SucVector)
    end
    @testset "rsearch" begin
        test_rsearch(SucVector)
    end
    @testset "long" begin
        test_long(SucVector)
    end
    @testset "copy" begin
        test_copy(SucVector)
    end
end

#=
facts("CompactBitVector") do
    context("access") do
        test_access(IndexableBitVectors.CompactBitVector)
    end
    context("rank") do
        test_rank(IndexableBitVectors.CompactBitVector)
    end
    context("select") do
        test_select(IndexableBitVectors.CompactBitVector)
    end
end
=#

@testset "RRR" begin
    @testset "access" begin
        test_access(RRR)
    end
    @testset "rank" begin
        test_rank(RRR)
    end
    @testset "select" begin
        test_select(RRR)
    end
    @testset "search" begin
        test_search(RRR)
    end
    @testset "rsearch" begin
        test_rsearch(RRR)
    end
    @testset "copy" begin
        test_copy(RRR)
    end
end

@testset "LargeRRR" begin
    @testset "access" begin
        test_access(IndexableBitVectors.LargeRRR)
    end
    @testset "rank" begin
        test_rank(IndexableBitVectors.LargeRRR)
    end
    @testset "select" begin
        test_select(IndexableBitVectors.LargeRRR)
    end
    @testset "copy" begin
        test_copy(IndexableBitVectors.LargeRRR)
    end
end

@testset "Size" begin
    @testset "compressible" begin
        bv = falses(10_000)
        @test sizeof(RRR(bv)) < sizeof(bv) < sizeof(SucVector(bv))
        bv = trues(10_000)
        @test sizeof(RRR(bv)) < sizeof(bv) < sizeof(SucVector(bv))
    end
    @testset "incompressible" begin
        bv = bitrand(10_000)
        @test sizeof(bv) < sizeof(SucVector(bv)) < sizeof(RRR(bv))
    end
end

ext: jl


text:
module GenericRound2

using Hecke
import AbstractAlgebra, Nemo
import Base: +, -, *, gcd, lcm, divrem, div, rem, mod, ^, ==
export integral_closure

mutable struct Order <: AbstractAlgebra.Ring
  F::AbstractAlgebra.Field
  R::AbstractAlgebra.Ring
  trans::MatElem #both matrices are over base_ring(F)
  itrans::MatElem

  function Order(R::AbstractAlgebra.Ring, F::AbstractAlgebra.Field, empty::Bool = false)
    r = new()
    r.F = F
    r.R = R
    empty && return r
    Qt = base_ring(F)
    d = lcm(map(x->denominator(x, R), coefficients(defining_polynomial(F))))
    f = map_coefficients(x->Qt(d*numerator(x, R)),defining_polynomial(F))
    if !ismonic(f) #need Lenstra Order
      d = degree(F)
      M = zero_matrix(Qt, d, d)
      M[1, 1] = one(Qt)
      for i in 2:d
        for j in i:-1:2
          M[i, j] = coeff(f, d - (i - j))
        end
      end
      O = Order(r, M, one(Qt))
      return O
    end
    return r
  end

  function Order(O::Order, T::MatElem, d::RingElem)
    F = base_ring(O.F)
    T = map_entries(F, T)
    T = divexact(T, base_ring(T)(d))
#    T = map(x->F(numerator(x)//denominator(x)), T)  #Rat{fmpq} is not simplifies #fixed in Bill's pull request in AA:canonicalisation
    Ti = inv(T)
#    Ti = map(x->F(numerator(x)//denominator(x)), Ti)#Rat{fmpq} is not simplifies
    r = Order(O.R, O.F, true)
    if isdefined(O, :trans)
      r.trans = T*O.trans
      r.itrans = O.itrans*Ti
    else
      r.trans = T
      r.itrans = Ti
    end
    return r
  end
end

function Base.show(io::IO, O::Order)
  println(io, "generic order over $(O.R) in $(O.F)")
end

Hecke.base_ring(O::Order) = O.R

mutable struct OrderElem <: RingElem
  parent :: Order
  data :: FieldElem
  coord :: Vector{RingElem}

  function OrderElem(O::Order, f::FieldElem)
    @assert parent(f) == O.F
    r = new()
    r.parent = O
    r.data = f
    return r
  end
  function OrderElem(O::Order, f::fmpz)
    return OrderElem(O, O.F(f))
  end

end

function Base.show(io::IO, a::OrderElem)
  print(io, a.data)
end

Nemo.elem_type(::Order) = OrderElem
Nemo.parent_type(::OrderElem) = Order
Nemo.parent_type(::Type{OrderElem}) = Order
Nemo.isdomain_type(::Type{OrderElem}) = true

Base.parent(a::OrderElem) = a.parent

(R::Order)(a::FieldElem) = OrderElem(R, a)
(R::Order)(a::fmpz) = OrderElem(R, a)
(R::Order)(a::Integer) = OrderElem(R, fmpz(a))
(R::Order)(a::OrderElem) = OrderElem(R, a.data)
(R::Order)() = R(0)

Nemo.iszero(a::OrderElem) = iszero(a.c)
Nemo.isone(a::OrderElem) = isone(a.c) && isone(a.f) && isone(a.g)

Nemo.zero(R::Order) = R(0)
Nemo.one(R::Order) = R(1)
Nemo.canonical_unit(a::OrderElem) = OrderElem(parent(a), fmpz(1))

Base.deepcopy_internal(a::OrderElem, dict::IdDict) = OrderElem(parent(a), Base.deepcopy_internal(a.data, dict))

+(a::OrderElem, b::OrderElem) = OrderElem(parent(a), a.data + b.data)
-(a::OrderElem, b::OrderElem) = OrderElem(parent(a), a.data - b.data)
-(a::OrderElem) = OrderElem(parent(a), -a.data)
*(a::OrderElem, b::OrderElem) = OrderElem(parent(a), a.data * b.data)

==(a::OrderElem, b::OrderElem) = parent(a) == parent(b) && a.data == b.data

function Hecke.mul!(a::OrderElem, b::OrderElem, c::OrderElem)
  a.data = Hecke.mul!(a.data, b.data, c.data)
  return a
end

function Hecke.add!(a::OrderElem, b::OrderElem, c::OrderElem)
  a.data = Hecke.add!(a.data, b.data, c.data)
  return a
end

function Hecke.addeq!(a::OrderElem, b::OrderElem)
  a.data = Hecke.addeq!(a.data, b.data)
  return a
end

function Hecke.tr(a::OrderElem)
  return parent(a).R(trace(a.data))
end

function Hecke.coordinates(a::OrderElem)
  c = coordinates(a.data)
  O = parent(a)
  if isdefined(O, :itrans)
    d = matrix(c)'*O.itrans
  else
    d = matrix(c)'
  end
  return d
end

function Hecke.coordinates(a::Generic.FunctionFieldElem)
  return [coeff(a, i) for i=0:degree(parent(a))-1]
end

Hecke.degree(O::Order) = degree(O.F)

function Hecke.basis(O::Order)
  if isdefined(O, :itrans)
    return [O(O.F(vec(collect(O.trans[i,:])))) for i=1:degree(O)]
  else
    return map(O, basis(O.F))
  end
end

function (O::Order)(c::Vector)
  if isdefined(O, :itrans)
    return O(O.F(vec(collect(matrix(map(base_ring(O.trans), c))'*O.trans))))
  else
    return O(O.F(c))
  end
end

struct FFElemCoeffs{T <: RingElement}
   f::T
end

function Nemo.coefficients(f::Generic.FunctionFieldElem)
   return FFElemCoeffs(f)
end

function Base.iterate(PC::FFElemCoeffs, st::Int = -1)
   st += 1
   if st > degree(parent(PC.f))
       return nothing
   else
       return coeff(PC.f, st), st
   end
end
Base.length(PC::FFElemCoeffs) = degree(parent(PC.f))

function Hecke.mod(a::OrderElem, p::RingElem)
  O = parent(a)
  R = parent(p)
  S = base_ring(O.F)

  if isdefined(O, :itrans)
    a = map(x->S(R(x) % p), coordinates(a))
    b = a*O.trans
    return O(O.F(vec(collect(b'))))
  else
    return O(O.F([O.R(x) % p for x = coefficients(a.data)]))
  end
end

Base.:^(a::OrderElem, n::Integer) = parent(a)(a.data^n)

function Hecke.powermod(a::OrderElem, n::fmpz, p::RingElem)
  c = parent(a)(1)
  for i = Hecke.BitsMod.bits(n)
    c *= c
    if i
      c *= a
    end
    c = Hecke.mod(c, p)
  end
  return c
end

function radical_basis_power(O::Order, p::RingElem)
  t = ResidueField(parent(p), p)
  if isa(t, Tuple)
    F, mF = t
  else
    F = t
    mF = MapFromFunc(x->F(x), y->lift(y), parent(p), F)
  end
#  @assert characteristic(F) == 0 || (isfinite(F) && characteristic(F) > degree(O))
  q = characteristic(F)
  @assert q > 0
  while q < degree(O)
    q *= characteristic(F)
  end

  b = basis(O)
  m = zero_matrix(F, degree(O), degree(O))
  for i=1:degree(O)
    c = coordinates(powermod(b[i], q, p))
    for j=1:degree(O)
      m[j,i] = mF(O.R(c[j]))
    end
  end
  k, B = kernel(m)

  M2 = B[:, 1:k]'
  M2 = map_entries(x->preimage(mF, x), M2)
  M3 = Hecke.hnf(vcat(M2, p*identity_matrix(parent(p), degree(O))))[1:degree(O), :]
  return return M3 #[O(vec(collect((M3[i, :])))) for i=1:degree(O)]
end

function Hecke.gcdinv(a::KInftyElem{fmpq}, b::KInftyElem{fmpq})
  g, q, w = gcdx(a, b)
  @assert isunit(g)
  return one(parent(a)), q*inv(g)
end

function Hecke.mul!(a::KInftyElem{fmpq}, b::KInftyElem{fmpq}, c::KInftyElem{fmpq})
  return b*c
end
function Hecke.addeq!(a::KInftyElem{fmpq}, b::KInftyElem{fmpq})
  return a+b
end

Hecke.isdomain_type(::Type{KInftyElem{T}}) where {T} = true

function radical_basis_trace(O::Order, p::RingElem)
  T = trace_matrix(O)

  t = ResidueField(parent(p), p)
  if isa(t, Tuple)
    R, mR = t
  else
    R = t
    mR = MapFromFunc(x->R(x), y->lift(y), parent(p), R)
  end

  TT = map_entries(mR, T)
  k, B = kernel(TT)
  M2 = map_entries(x->preimage(mR, x), B[:, 1:k])'
  M3 = Hecke.hnf(vcat(M2, p*identity_matrix(parent(p), degree(O))))[1:degree(O), :]
  return return M3 #[O(vec(collect((M3[i, :])))) for i=1:degree(O)]
end

function radical_basis_power_non_perfect(O::Order, p::RingElem)
  t = ResidueField(parent(p), p)
  if isa(t, Tuple)
    F, mF = t
  else
    F = t
    mF = MapFromFunc(x->F(x), y->lift(y), parent(p), F)
  end
  @assert isa(F, Generic.RationalFunctionField) && characteristic(F) != 0
#  @assert characteristic(F) == 0 || (isfinite(F) && characteristic(F) > degree(O))
  q = characteristic(F)
  @assert q > 0
  while q < degree(O)
    q *= characteristic(F)
  end

  q = Int(q)
  b = basis(O)
  m = zero_matrix(F, q*degree(O), degree(O))
  for i=1:degree(O)
    c = coordinates(powermod(b[i], fmpz(q), p))
    for j=1:degree(O)
      d = mF(O.R(c[j]))
      @assert isone(denominator(d))
      d = numerator(d)
      for k=0:degree(d)
        iszero(coeff(d, k)) && continue
        m[(j-1)*q+(k%q)+1,i] += gen(parent(d))^div(k, q)*coeff(d, k)
      end
    end
  end
  k, B = kernel(m)

  M2 = B[:, 1:k]'
  M2 = map_entries(x->preimage(mF, x), M2)
#  M2 = map_entries(x->preimage(mF, _root(numerator(x), q)//_root(denominator(x), q)), B[:, 1:k])'
  M3 = Hecke.hnf(vcat(M2, p*identity_matrix(parent(p), degree(O))))[1:degree(O), :]
  return return M3 #[O(vec(collect((M3[i, :])))) for i=1:degree(O)]

end

function Hecke.representation_matrix(a::OrderElem)
  O = parent(a)
  b = basis(O)
  m = zero_matrix(O.R, degree(O), degree(O))
  for i=1:degree(O)
    c = coordinates(b[i]*a)
    for j=1:degree(O)
      m[i,j] = c[j]
    end
  end
  return m
end

function Hecke.representation_matrix(a::Generic.FunctionFieldElem)
  F = parent(a)
  g = gen(F)
  m = zero_matrix(base_ring(F), degree(F), degree(F))
  b = a
  for i=1:degree(F)
    for j=1:degree(F)
      m[i,j] = coeff(b, j-1)
    end
    if i<degree(F)
      b *= g
    end
  end
  return m
end

function ring_of_multipliers(O::Order, I::MatElem)
  #TODO: modular big hnf, peu-a-peu, not all in one
  II, d = pseudo_inv(I)
#  return II, d, [representation_matrix(O(vec(collect(I[i, :])))) for i=1:nrows(I)]
  m = hcat([divexact(representation_matrix(O(vec(collect(I[i, :]))))*II, d) for i=1:nrows(I)]...)
  H = hnf(m')[1:degree(O), 1:degree(O)]
  Hi, d = pseudo_inv(H)

  O = Order(O, Hi', d)
  return O
end

function Hecke.trace_matrix(O::Order)
  return trace_matrix(basis(O))
end

function Hecke.trace_matrix(b::Vector{OrderElem})
  O = parent(b[1])
  m = zero_matrix(O.R, length(b), length(b))
  for i=1:length(b)
    m[i,i] = trace(b[i]*b[i])
    for j=i+1:length(b)
      m[i,j] = m[j,i] = trace(b[i]*b[j])
    end
  end
  return m
end

function Hecke.trace_matrix(b::Vector{OrderElem}, c::Vector{OrderElem}, exp::fmpz = fmpz(1))
  O = parent(b[1])
  m = zero_matrix(O.R, length(b), length(c))
  for i=1:length(b)
    for j=1:length(c)
      m[i,j] = trace((b[i]*c[j])^exp)
    end
  end
  return m
end

function Hecke.pmaximal_overorder(O::Order, p::RingElem)
  t = ResidueField(parent(p), p)

  if isa(t, Tuple)
    R, mR = t
  else
    R = t
    mR = MapFromFunc(x->R(x), y->lift(y), parent(p), R)
  end
#  @assert characteristic(F) == 0 || (isfinite(F) && characteristic(F) > degree(O))
  if characteristic(R) == 0 || characteristic(R) > degree(O)
    rad = radical_basis_trace
  elseif isa(R, Generic.RationalFunctionField)
    rad = radical_basis_power_non_perfect
  else
    rad = radical_basis_power
  end
  while true
    I = rad(O, p)
    S = ring_of_multipliers(O, I)
    if discriminant(O) == discriminant(S)
      return O
    end
    O = S
  end
end

function integral_closure(S::AbstractAlgebra.Ring, F::Generic.FunctionField)
  O = Order(S, F)
  d = discriminant(O)
  ld = factor(d)
  local Op
  first = true
  for (p,k) = ld.fac
    if k<2
      continue
    end
    T = integral_split(pmaximal_overorder(O, p).trans, S)
    isone(T[2]) && continue
    if first
      Op = T
      first = false
    else
      Op = (Hecke._hnf(vcat(Op[1]*T[2], T[1]*Op[2]), :lowerleft)[degree(O)+1:end, :], T[2]*Op[2])
    end
  end
  if first
    return O
  else
    return Order(O, Op[1], Op[2])
  end
end

#goal IntCls(Z[x], F) as the intersection of IntCls(Z<x>) and IntCls(Q[x])
#ala Hess

function Hecke.discriminant(O::Order)
  d = discriminant(O.F)
  if isdefined(O, :trans)
    d *= det(O.trans)^2
  end
  return O.R(d)
end

function Hecke.basis(O::Order, F::Generic.FunctionField)
  return map(F, basis(O))
end
function Hecke.basis(F::Generic.FunctionField)
  a = gen(F)
  bas = [a^0, a]
  while length(bas) < degree(F)
    push!(bas, bas[end]*a)
  end
  return bas
end

(R::PolyRing{T})(a::Generic.Rat{T}) where {T} = R(numerator(a))

function Hecke.ResidueField(R::FmpqPolyRing, p::fmpq_poly)
  K, _ = number_field(p)
  return K, MapFromFunc(x->K(x), y->R(y), R, K)
end

function (F::Generic.FunctionField)(p::PolyElem)
  @assert parent(p) == parent(F.pol)
  @assert degree(p) < degree(F) # the reduction is not implemented
  R = parent(gen(F).num)
  S = base_ring(R)
  d = reduce(lcm, map(denominator, coefficients(p)), init = one(S))
  return F(map_coefficients(S, d*p, parent = R), d)
end

function (F::Generic.FunctionField)(c::Vector{<:RingElem})
  return F(parent(F.pol)(map(base_ring(F), c)))
end

(F::Generic.FunctionField)(a::OrderElem) = a.data

function Hecke.charpoly(a::OrderElem)
  return charpoly(representation_matrix(a))
end

function Hecke.minpoly(a::OrderElem)
  return minpoly(representation_matrix(a))
end

function Hecke.charpoly(a::Generic.FunctionFieldElem)
  return charpoly(representation_matrix(a))
end

function Hecke.minpoly(a::Generic.FunctionFieldElem)
  return minpoly(representation_matrix(a))
end

function Hecke.discriminant(F::Generic.FunctionField)
  return discriminant(defining_polynomial(F))
end

function Hecke.factor(a::Generic.Rat{T}, R::Generic.PolyRing{T}) where {T}
  @assert parent(numerator(a)) == R
  f1 = factor(numerator(a))
  f2 = factor(denominator(a))
  for (p,e) = f2.fac
    @assert !haskey(f1.fac, p)
    f1.fac[p] = -e
  end
  f1.unit = divexact(f1.unit, f2.unit)
  return f1
end

Hecke.characteristic(R::KInftyRing) = characteristic(R.K)

function Hecke.factor(d::KInftyElem)
  t = gen(parent(d))
  a = degree(d)
  return Fac(d*t^a, Dict(t=>-a))
end

function Hecke.numerator(a::Generic.Rat, S::FmpqPolyRing)
  return numerator(a)
end

function Hecke.denominator(a::Generic.Rat, S::FmpqPolyRing)
  return denominator(a)
end

function Hecke.integral_split(a::Generic.Rat, S::FmpqPolyRing)
  return numerator(a), denominator(a)
end

Base.denominator(x::AbstractAlgebra.Generic.Rat{fmpq}, R::KInftyRing{fmpq}) = Hecke.integral_split(x, R)[2]
Base.numerator(x::AbstractAlgebra.Generic.Rat{fmpq}, R::KInftyRing{fmpq}) = Hecke.integral_split(x, R)[1]

function Hecke.integral_split(x::AbstractAlgebra.Generic.Rat{fmpq}, R::KInftyRing{fmpq})
  if iszero(x)
    return zero(R), one(R)
  end
  a = degree(numerator(x))
  b = degree(denominator(x))
  if b >= a
    return R(x), one(R)
  end
  t = gen(parent(x))
  return R(x*t^(b-a)), R(t^(b-a))
end

function Hecke.lcm(a::KInftyElem{fmpq}, b::KInftyElem{fmpq})
  return gen(parent(a))^-min(degree(a), degree(b))
end
function Hecke.gcd(a::KInftyElem{fmpq}, b::KInftyElem{fmpq})
  return gen(parent(a))^-max(degree(a), degree(b))
end

Hecke.gen(R::KInftyRing) = R(inv(gen(R.K)))

(R::Generic.RationalFunctionField{fmpq})(x::KInftyElem{fmpq}) = x.d

function Hecke.integral_split(M::Generic.MatElem{<:Generic.Rat}, S::Generic.Ring)
  m = zero_matrix(S, nrows(M), ncols(M))
  den = one(S)
  for i=1:nrows(M)
    for j=1:ncols(M)
      n, d = integral_split(M[i,j], S)
      q, r = divrem(den, d)
      if iszero(r)
        m[i,j] = n*q
      else
        t = lcm(den, d)
        m *= divexact(t, den)
        m[i,j] = n*divexact(t, d)
        den = t
      end
    end
  end
  return m, den
end

end  # ModuleRound2

module HessQRModule
using Hecke
import AbstractAlgebra, Nemo
import Base: +, -, *, gcd, lcm, divrem, div, rem, mod, ^, ==
export HessQR

struct HessQR <: AbstractAlgebra.Ring
  R::FmpzPolyRing
  Qt::Generic.RationalFunctionField{fmpq}

  function HessQR(R::FmpzPolyRing, Qt::Generic.RationalFunctionField)
    new(R, Qt)
  end
end

mutable struct HessQRElem <: RingElem
  parent::HessQR
  c::fmpz
  f::fmpz_poly
  g::fmpz_poly

  function HessQRElem(P::HessQR, c::fmpz, f::fmpz_poly, g::fmpz_poly)
    if iszero(c) || iszero(f)
      r = new(P, fmpz(0), zero(P.R), one(P.R))
      @assert parent(r.f) == P.R
      @assert parent(r.g) == P.R
      return r
    end
    if parent(f) != P.R
      f = map_coefficients(FlintZZ, f, parent = P.R)
    end
    if parent(g) != P.R
      g = map_coefficients(FlintZZ, g, parent = P.R)
    end
    gc = gcd(f, g)
    f = divexact(f, gc)
    g = divexact(g, gc)
    cf = content(f)
    cg = content(g)
    @assert (c*cf) % cg == 0
    cu = canonical_unit(g)
    r = new(P, divexact(c*cf, cg), cu*divexact(f, cf), cu*divexact(g, cg))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
  function HessQRElem(P::HessQR, q::Generic.Rat{fmpq})
    f = numerator(q)
    g = denominator(q)
    return HessQRElem(P, f, g)
  end
  function HessQRElem(P::HessQR, f::fmpq_poly)
    return HessQRElem(P, f, one(parent(f)))
  end
  function HessQRElem(P::HessQR, f::fmpq_poly, g::fmpq_poly)
    df = reduce(lcm, map(denominator, coefficients(f)), init = fmpz(1))
    dg = reduce(lcm, map(denominator, coefficients(g)), init = fmpz(1))
    ff = map_coefficients(FlintZZ, df*f, parent = P.R)
    gg = map_coefficients(FlintZZ, dg*g, parent = P.R)
    #ff/df//gg/dg = dg/df * ff/gg
    return HessQRElem(P, divexact(dg, df), ff, gg)
  end
  function HessQRElem(P::HessQR, c::fmpz)
    r = new(P, c, one(P.R), one(P.R))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
  function HessQRElem(P::HessQR, c::fmpz_poly)
    d = content(c)
    r = new(P, d, divexact(c, d), one(P.R))
    @assert parent(r.f) == P.R
    @assert parent(r.g) == P.R
    return r
  end
end

function Base.show(io::IO, a::HessQRElem)
  print(io, "$(a.c) * ($(a.f))//($(a.g))")
end

function Hecke.integral_split(a::Generic.Rat{fmpq}, S::HessQR)
  if iszero(a)
    return zero(S), one(S)
  end
  n = numerator(a)
  d = denominator(a)
  dn = reduce(lcm, map(denominator, coefficients(n)), init = fmpz(1))
  dd = reduce(lcm, map(denominator, coefficients(d)), init = fmpz(1))
  zn = S.R(n*dn)
  zd = S.R(d*dd)
  cn = content(zn)
  cd = content(zd)
  zn = divexact(zn, cn)
  zd = divexact(zd, cd)
  cn *= dd
  cd *= dn
  g = gcd(cn, cd)
  cn = divexact(cn, g)
  cd = divexact(cd, g)
  return HessQRElem(S, cn, zn, zd), S(cd)
end

function Hecke.numerator(a::Generic.Rat, S::HessQR)
  return integral_split(a, S)[1]
end

function Hecke.denominator(a::Generic.Rat, S::HessQR)
  return integral_split(a, S)[2]
end

Nemo.elem_type(::HessQR) = HessQRElem
Nemo.elem_type(::Type{HessQR}) = HessQRElem
Nemo.parent_type(::HessQRElem) = HessQR
Nemo.parent_type(::Type{HessQRElem}) = HessQR
Nemo.isdomain_type(::Type{HessQRElem}) = true

Base.parent(a::HessQRElem) = a.parent

(R::HessQR)(a::Generic.Rat{fmpq}) = HessQRElem(R, a)
(R::HessQR)(a::fmpz) = HessQRElem(R, a)
(R::HessQR)(a::Integer) = HessQRElem(R, fmpz(a))
(R::HessQR)(a::fmpz_poly) = HessQRElem(R, a)
(R::HessQR)(a::fmpq_poly) = HessQRElem(R, a)
(R::HessQR)(a::HessQRElem) = a
(R::HessQR)() = R(0)

(F::Generic.RationalFunctionField)(a::HessQRElem) = a.c*F(a.f)//F(a.g)

Nemo.iszero(a::HessQRElem) = iszero(a.c)
Nemo.isone(a::HessQRElem) = isone(a.c) && isone(a.f) && isone(a.g)

Nemo.zero(R::HessQR) = R(0)
Nemo.one(R::HessQR) = R(1)
Nemo.canonical_unit(a::HessQRElem) = HessQRElem(parent(a), fmpz(1), a.f, a.g)

Base.deepcopy_internal(a::HessQRElem, dict::IdDict) = HessQRElem(parent(a), Base.deepcopy_internal(a.c, dict), Base.deepcopy_internal(a.f, dict), Base.deepcopy_internal(a.g, dict))

Base.hash(a::HessQRElem, u::UInt=UInt(12376599)) = hash(a.g, hash(a.f, hash(a.c, u)))

+(a::HessQRElem, b::HessQRElem) = HessQRElem(parent(a), fmpz(1), a.c*a.f*b.g+b.c*b.f*a.g, a.g*b.g)
-(a::HessQRElem, b::HessQRElem) = HessQRElem(parent(a), fmpz(1), a.c*a.f*b.g-b.c*b.f*a.g, a.g*b.g)
-(a::HessQRElem) = HessQRElem(parent(a), -a.c, a.f, a.g)
*(a::HessQRElem, b::HessQRElem) = HessQRElem(parent(a), a.c*b.c, a.f*b.f, a.g*b.g)

==(a::HessQRElem, b::HessQRElem) = parent(a) == parent(b) && a.c == b.c && a.f == b.f && a.g == b.g

Base.:^(a::HessQRElem, n::Int) = HessQRElem(parent(a), a.c^n, a.f^n, a.g^n)

function Hecke.mul!(a::HessQRElem, b::HessQRElem, c::HessQRElem)
  d = b*c
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function Hecke.add!(a::HessQRElem, b::HessQRElem, c::HessQRElem)
  d = b+c
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function Hecke.addeq!(a::HessQRElem, b::HessQRElem)
  d = a+b
  @assert parent(a.f) == parent(d.f)
  @assert parent(a.g) == parent(d.g)
  a.c = d.c
  a.f = d.f
  a.g = d.g
  return a
end

function divrem(a::HessQRElem, b::HessQRElem)
  if iszero(b)
    error("div by 0")
  end
  if iszero(a)
    return a, a
  end
  #= af/g mod b:
    it is enough to figure this out for b in Z>0, the rest is units
    that will disappaear into the quotient q

    - c = cont(g mod d), then since cont(g) is 1, c is a unit mod d
    - cd = 1 mod d, then
    - cont(g(cx^?+1) mod d) = 1 if ? = deg(g)+1 (or larger):
      g(cx?+1) = cg x^? + g and cont(cg mod d = 1)...
    - af/g -d/(cx^?+1) now has a denominator with cont( mod d) = 1
    - af/g - (af - du)/(g - dv) =
      (af(g-dv) -g(af-du))/(g(g-dv)) = d*..., so once cont(g %d) =1,
      we can replace f and g mod d (and figure out the quotient afterwards)

    - for d = p a prime, the rep is unqiue, thus F_p(t)
  =#
  r = rem(a,b)
  return divexact(a-r, b), r

  return HessQRElem(parent(a), q, a.f*b.g, a.g*b.f), HessQRElem(parent(a), r, a.f, a.g)
end

function div(a::HessQRElem, b::HessQRElem)
  if iszero(a)
    return a
  end
  return divexact(a-rem(a, b), b)
  q = div(a.c, b.c)
  return HessQRElem(parent(a), q, a.f*b.g, a.g*b.f)
end

function rem(a::HessQRElem, b::HessQRElem)
  if iszero(a)
    return a
  end
  #see above for reasoning
  d = abs(b.c)
  if isone(d)
    return zero(parent(a))
  end
  R = parent(a).R
  gd = mod(a.g, d)
  c = content(gd)
  ci = invmod(c, d)
  e = gen(parent(gd))^(degree(a.g)+1)+1
  f = a.c*a.f*e
  g = a.g*e
  gd = mod(g, d)
  @assert content(gd) == 1

  fd = mod(f, d)

  return HessQRElem(parent(a), fmpz(1), fd, gd)
end

function Nemo.divexact(a::HessQRElem, b::HessQRElem)
  @assert parent(a) == parent(b)
  @assert parent(a.f) == parent(a).R
  @assert parent(a.g) == parent(a).R
  @assert parent(b.f) == parent(a).R
  @assert parent(b.g) == parent(a).R
  return HessQRElem(parent(a), divexact(a.c, b.c), a.f*b.g, a.g*b.f)
end

function gcd(a::HessQRElem, b::HessQRElem)
  return HessQRElem(parent(a), gcd(a.c, b.c))
end

function Nemo.gcdx(a::HessQRElem, b::HessQRElem)
  R = parent(a)
  g, u, v = Nemo.gcdx(a.c, b.c)
  return R(g), HessQRElem(R, u, a.g, a.f), HessQRElem(R, v, b.g, b.f)
end

function lcm(a::HessQRElem, b::HessQRElem)
  return HessQRElem(parent(a), lcm(a.c, b.c))
end

Hecke.isunit(a::HessQRElem) = isunit(a.c)

Nemo.dense_poly_type(::Type{gfp_fmpz_elem}) = gfp_fmpz_poly

function Nemo.ResidueField(a::HessQR, b::HessQRElem)
  @assert isprime(b.c)
  F = GF(b.c)
  Ft, t = RationalFunctionField(F, String(var(a.R)), cached = false)
  R = parent(numerator(t))
  return Ft, MapFromFunc(x->F(x.c)*Ft(map_coefficients(F, x.f, parent = R))//Ft(map_coefficients(F, x.g, parent = R)),
                         y->HessQRElem(a, fmpz(1), map_coefficients(lift, numerator(y)), map_coefficients(lift, denominator(y))), a, Ft)
end

function Nemo.ResidueRing(a::HessQR, b::HessQRElem)
  error("wrong wrong wrong")
  F = ResidueRing(FlintZZ, b.c)
  return F, MapFromFunc(x->F(x.c), y->a(lift(y)), a, F)
end

function Hecke.factor(a::HessQRElem)
  f = factor(a.c)
  R = parent(a)
  return Fac(R(a.f), Dict((R(p),k) for (p,k) = f.fac))
end

function Hecke.factor(a::Generic.Rat, R::HessQR)
  d1 = reduce(lcm, map(denominator, coefficients(numerator(a))), init = fmpz(1))
  f1 = factor(R(d1*numerator(a)))
  d2 = reduce(lcm, map(denominator, coefficients(denominator(a))), init = fmpz(1))
  f2 = factor(R(d1*denominator(a)))

  for (p,k) = f2.fac
    if haskey(f1.fac, p)
      f1.fac[p] -= k
    else
      f1.fac[p] = k
    end
  end
  f1.unit = divexact(f1.unit, f2.unit)
  return f1
end

end

using .HessQRModule
using .GenericRound2

export HessQR, integral_closure

module HessMain
using Hecke
using ..HessQRModule, ..GenericRound2

function GenericRound2.integral_closure(Zx::FmpzPolyRing, F::Generic.FunctionField)
  Qt = base_ring(F)
  t = gen(Qt)
  S = HessQR(Zx, Qt)
  o1 = integral_closure(S, F)
  o2 = integral_closure(parent(denominator(t)), F)
  T = o1.trans * o2.itrans
  q, w = integral_split(T, S)
  h, T1 = hnf_with_transform(q)
  return basis(GenericRound2.Order(o1, T1, one(S)), F)
  #the matrix H below has to be the identity and the 2 order thus then have the
  #same basis. Hence no reason to compute...
  T = map_entries(x->Qt(x)//Qt(w), h)
  qq, ww = integral_split(T', base_ring(o2))
  H, T2 = hnf_with_transform(qq)
  return H, GenericRound2.Order(o1, T1, one(S)), GenericRound2.Order(o2, inv(T2'), one(base_ring(T1)))
end

end

using .HessMain


text:
using DLSMethod
using FITSIO
using HDF5
using LinearAlgebra
using Statistics

include("BayesianSZNet.jl")
using .BayesianSZNet

ε = 0.00005

@time begin
    fid = h5open("data/dr16q_superset.hdf5", "r+")
    id = read(fid["id"])
    n = size(id, 2)

    fluxes = zeros(3752, n)

    Threads.@threads for i = 1:n
        loglam, flux = get_spectrum("DR16Q_Superset_v3", id[:, i]...)

        # standardise flux
        flux_mean = mean(flux)
        standard_flux = Vector{Float64}((flux .- flux_mean) ./ std(flux, mean=flux_mean))

        X = polynomial_features(loglam, 3)
        n_best, a_best = dls_fit(X, standard_flux, 2, 0.9)

        # padding
        idx = (LOGLAMMIN - ε .< loglam) .& (loglam .< LOGLAMMAX + ε)
        pixels = sum(idx)
        Δ_loglam = minimum(loglam) - LOGLAMMIN
        offset = Δ_loglam <= 0 ? 1 : round(Int, Δ_loglam * 1e4) + 1
        fluxes[offset:offset + pixels - 1, i] = (standard_flux - X * a_best)[idx]
    end

    write_dataset(fid, "X", fluxes)
    close(fid)
end

max_stars_repo_name: podondra/bdl-redshift

max_issues_repo_name: podondra/bdl-redshift

max_forks_repo_name: podondra/bdl-redshift


text:
using Test
using MIPVerify: get_example_network_params, read_datasets, frac_correct
@isdefined(TestHelpers) || include("../TestHelpers.jl")

TestHelpers.@timed_testset "import_example_nets.jl" begin
    TestHelpers.@timed_testset "get_example_network_params" begin
        TestHelpers.@timed_testset "MNIST.n1" begin
            nn = get_example_network_params("MNIST.n1")
            mnist = read_datasets("mnist")
            @test frac_correct(nn, mnist.test, 10000) == 0.9695
        end
        TestHelpers.@timed_testset "MNIST.WK17a_linf0.1_authors" begin
            nn = get_example_network_params("MNIST.WK17a_linf0.1_authors")
            mnist = read_datasets("mnist")
            @test frac_correct(nn, mnist.test, 10000) == 0.9811
        end
        TestHelpers.@timed_testset "MNIST.RSL18a_linf0.1_authors" begin
            nn = get_example_network_params("MNIST.RSL18a_linf0.1_authors")
            mnist = read_datasets("mnist")
            @test frac_correct(nn, mnist.test, 10000) == 0.9582
        end
    end

    @testset "unrecognized example network" begin
        @test_throws DomainError get_example_network_params("the social network")
    end
end


text:
using DataValues

a = [DataValue(3.), DataValue(2.), DataValue(5.)]

a[a .> 2.]


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

import Base: ∘
"""
 `∘(A::WFST,B::WFST; filter=EpsSeq, connect=true)`

Perform composition of the transucers `A` and `B`.

```julia
julia> A = linearfst(["a","b","c"],["α","β","γ"],ones(3))
WFST #states: 4, #arcs: 3, #isym: 3, #osym: 3
|1/0.0f0|
a:α/1.0f0 → (2)
(2)
b:β/1.0f0 → (3)
(3)
c:γ/1.0f0 → (4)
((4/0.0f0))

julia> B = matrix2wfst(["α","β","γ"],[:w,:x,:z],ones(3,3))
WFST #states: 4, #arcs: 9, #isym: 3, #osym: 3
|1/0.0f0|
α:w/1.0f0 → (2)
β:x/1.0f0 → (2)
γ:z/1.0f0 → (2)
(2)
α:w/1.0f0 → (3)
β:x/1.0f0 → (3)
γ:z/1.0f0 → (3)
(3)
α:w/1.0f0 → (4)
β:x/1.0f0 → (4)
γ:z/1.0f0 → (4)
((4/0.0f0))

julia> A∘B
WFST #states: 4, #arcs: 3, #isym: 3, #osym: 3
|1/0.0f0|
a:w/2.0f0 → (2)
(2)
b:x/2.0f0 → (3)
(3)
c:z/2.0f0 → (4)
((4/0.0f0))

```

The keyword `filter` can specify the composition filter to be used, which makes it possible to handle epsilon-transitions.
See [Allauzen et al. "Filters for Efficient Composition of Weighted Finite-State Transducers"](https://storage.googleapis.com/pub-tools-public-publication-data/pdf/36838.pdf).

If `connect` is set to `true` after completing the composition the [`connect`](@ref) algorithm is applied.

"""
function ∘(A::WFST{W1,A1,D,I1,O1}, B::WFST{W2,A2,D,I2,O2};
           filter=EpsSeq, connect=true
          ) where {W1, W2, D, I1, I2, O1, O2,
                   A1 <: Arc{W1,D},
                   A2 <: Arc{W2,D}}
  if O1 != I2
    error("Output label of A must be of the same type of input labels of B.")
  end
  if W1 != W2
    error("Weights of A must be of the same type of weights of B.")
  end
  if !(A.osym === B.isym)
    if A.osym != B.isym # construncting A and B using the same object avoids this check
      throw(ErrorException("Output symbol table of A not consistent with input symbol table of B."))
    end
  end

  function updateC!(it::ComposeIterator,
                    i::D,q1::D,q2::D,q3::D,w3::W,
                    e1::Arc{W,D},e2::Arc{W,D},
                    src::D,dest::D) where {W,A,D}
    wc = get_weight(e1)*get_weight(e2)
    if max(src,dest) - length(it.C) > 0
      add_states!(it.C, max(src,dest) - length(it.C) )  # eventually expand states
    end
    push_arc!(it.C,src,Arc(get_ilabel(e1),get_olabel(e2),wc,dest))
  end

  function updateCfinal!(it::ComposeIterator,
                         i::D,q1::D,q2::D,q3::D,w3::W,
                         src::D) where {W,A,D}
    if isinitial(it.A,q1) && isinitial(it.B,q2)
      w1, w2 = get_initial_weight(it.A,q1), get_initial_weight(it.B,q2)
      initial!(it.C, src, w1*w2)
    end
    if isfinal(it.A,q1) && isfinal(it.B,q2) &&
      q3 <= it.filter.q_max && w3 != zero(typeof(w3))
      w1, w2 = get_final_weight(it.A,q1), get_final_weight(it.B,q2)
      final!(it.C, src, w1*w2*w3)
    end
  end

  ϵ  = 0
  ϵL = max(length(get_osym(A)),length(get_isym(B)))+1
  f = filter{W1,D}(ϵL)
  # add self loops
  for i in 1:size(A,1)
    push_arc!(A,i,Arc{W1,D}(ϵ,ϵL,one(W1),i))
  end
  for i in 1:size(B,1)
    push_arc!(B,i,Arc{W1,D}(ϵL,ϵ,one(W1),i))
  end
  C = WFST(get_isym(A),get_osym(B); W = W1)
  z = ComposeIterator(A,B,C,f,updateC!,updateCfinal!)
  for zi in z
    nothing
  end
  # delete self loops
  for s in get_states(A) deleteat!(s,length(s)) end
  for s in get_states(B) deleteat!(s,length(s)) end
  if connect
    connect!(C)
  end
  return C
end


text:
using DataFrames
using Unitful

import Unitful:
    nm, μm, mm, cm, ns, μs, ms, ps, s

"""
    setunits(df::DataFrame)

Return a dataframe for lor computation where quantities are defined with their units
"""

function setunits(df::DataFrame)

    DataFrame(
        nsipm1 = df.nsipm1,     # number of sipms in cluster
        nsipm2 = df.nsipm2,
        phistd1 = df.phistd1,   # although these have units they are never used as
        phistd2 = df.phistd2,   # quantities, but as an estimator, so leave it implicit.
        zstd1 = df.zstd1,       # ditto
        zstd2 = df.zstd2,

        q1 = df.q1,             # q in pes
        q2 = df.q2,

        r1  = df.r1  * mm,      # best radius
        r1x = df.r1x * mm,      # radius from best estimator
        r2  = df.r2  * mm,
        r2x = df.r2x * mm,

        t1  = df.t1  * ns,      # time true
        t2  = df.t2  * ns,
        ta1 = df.ta1 * ns,      # average (over 5 mimimum in t)
        ta2 = df.ta2 * ns,
        tr1 = df.tr1 * ns,      # reco (smeared take minimum)
        tr2 = df.tr2 * ns,

        ux = df.ux,             # unit direction vector of gammas
        uy = df.uy,
        uz = df.uz,

        x1  = df.x1  * mm,      # best reco position
        x2  = df.x2  * mm,
        xb1 = df.xb1 * mm,      # position of sipm that gives time stamp
        xb2 = df.xb2 * mm,
        xr1 = df.xr1 * mm,      # reco position
        xr2 = df.xr2 * mm,
        xs  = df.xs  * mm,      # position of source
        xt1 = df.xt1 * mm,      # true position
        xt2 = df.xt2 * mm,

        y1  = df.y1  * mm,
        y2  = df.y2  * mm,
        yb1 = df.yb1 * mm,
        yb2 = df.yb2 * mm,
        yr1 = df.yr1 * mm,
        yr2 = df.yr2 * mm,
        ys  = df.ys  * mm,
        yt1 = df.yt1 * mm,
        yt2 = df.yt2 * mm,

        z1  = df.z1  * mm,
        z2  = df.z2  * mm,
        zb1 = df.zb1 * mm,
        zb2 = df.zb2 * mm,
        zr1 = df.zr1 * mm,
        zr2 = df.zr2 * mm,
        zs  = df.zs  * mm,
        zt1 = df.zt1 * mm,
        zt2 = df.zt2 * mm)
end

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
# Julia script for benchmarking binary signal operators

using BenchmarkTools
using ProgressBars
using DataFrames
using JLD2

import Romeo

# Benchmark Convolution
# ---------------------
# func = Romeo.LTI.Operators.convolve
# filename = "benchmarks/runs/conv_N1000.jld2"

# Benchmark Convolution
# ---------------------
func = Romeo.LTI.Operators.deconvolve
filename = "benchmarks/runs/deconv_N1000.jld2"

types = [Int32, Int64, Float32, Float64, ComplexF32, ComplexF64]
results = DataFrame(type=Type[], trial=BenchmarkTools.Trial[])
N = 1000
types_iter = ProgressBar(types)
for T in types_iter
    set_postfix(types_iter, type="$T")
    if T <: Integer
        g = rand(-10 * one(T):10 * one(T), N)
        u = rand(-10 * one(T):10 * one(T), N)
    else
        g = rand(T, N)
        u = rand(T, N)
    end
    r = @benchmark func($g, $u)
    push!(results, (T, copy(r)))
end

show(results)

jldsave(filename; results)


text:
include("solver_preamble.jl")
using MosekTools
optimizer_constructor = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true)
config = MOI.Test.Config(atol=1e-3, rtol=1e-3)
@testset "Mosek" begin
    @testset "Square" begin
        Tests.square_test(optimizer_constructor, config)
    end
    @testset "Invariant in Square" begin
        Tests.invariant_test(optimizer_constructor, config)
    end
    @testset "Controlled Invariant in Square" begin
        Tests.ci_test(optimizer_constructor, config)
    end
    @testset "Switched System" begin
        Tests.switched_system_test(optimizer_constructor, config)
    end
end


text:
using BitInformation
using StatsBase

function get_bitinformation(X::AbstractArray{T}; kwargs...) where {T<:Base.IEEEFloat}

    BitInformation.signed_exponent!(X)
    IC = bitinformation(X; kwargs...)

    return IC

end

function get_bitinformation(X::AbstractArray{T}; kwargs...) where {T<:Union{Int16,Int32,Int64}}

    IC = bitinformation(X; kwargs...)

    return IC

end

function get_keepbits(bitinfo_dict)
    # filter insignificant information via binomial distribution as explained in methods
    # for some variables the last bits do not have a free entropy of 1 bit, in that case
    # adjust the threshold by the information in the very last bits (which should always be 0)
    bitinfo = bitinfo_dict["bitinfo"]
    maskinfo = bitinfo_dict["maskinfo"]
    inflevel = bitinfo_dict["inflevel"]

    ic = values(bitinfo)
    p = BitInformation.binom_confidence(maskinfo,0.99)  # get chance p for 1 (or 0) from binom distr
    M₀ = 1 - entropy([p,1-p],2)                            # free entropy of random 50/50 at trial size
    threshold = max(M₀,1.5*maximum(ic[end-3:end]))         # in case the information never drops to zero
                                                               # use something a bit bigger than maximum
                                                               # of the last 4 bits
    insigni = (ic .<= threshold) .& (collect(1:length(ic)) .> 9)
    ic[insigni] .= floatmin(Float64)

    # find bits with 99/100% of information
    ICcsum = cumsum(ic)
    ICcsum_norm = copy(ICcsum)
    ICcsum_norm[:] ./= ICcsum_norm[end]

    infbits = argmax(ICcsum_norm[:] .> inflevel)
    #infbits100 = [argmax(ICcsum_norm[i,:] .> 0.999999999) for i in 1:nvars];

    return infbits

end


text:
module PandorasBoxPropellerBox

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/propellerBox" PropellerBox(x::Integer, y::Integer, texture::String="default", flashUseColor::String="3F437C", flashChargeColor::String="5A1C1C", maxCharges::Integer=3, rechargeOnGround::Bool=true, glideMode::String="AfterUse")

const placements = Ahorn.PlacementDict(
    "Propeller Box (Pandora's Box)" => Ahorn.EntityPlacement(
        PropellerBox
    )
)

function getSprite(entity::PropellerBox)
    texture = get(entity, "texture", "default")

    return "objects/pandorasBox/propellerBox/$texture/default_charges00"
end

Ahorn.editingOptions(entity::PropellerBox) = Dict{String, Any}(
    "glideMode" => Dict{String, String}(
        "Always" => "Always",
        "After Using Charge" => "AfterUse",
        "Never" => "Never"
    )
)

function Ahorn.selection(entity::PropellerBox)
    x, y = Ahorn.position(entity)
    sprite = getSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PropellerBox, room::Maple.Room) = Ahorn.drawSprite(ctx, getSprite(entity), 0, 0, jx=0.5, jy=1.0)

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
# Holds a copy of the problem passed to EAGO
# Needs separate storage for the MOI compilant model and the nlp_data
Base.@kwdef mutable struct InputModel{T<:Real}
    _input_model::MOIU.Model{T} = MOIU.Model{T}()
    _nlp_data::Union{MOI.NLPBlockData, Nothing} = nothing
end

MOI.is_empty(m::InputModel) = (m._nlp_data == nothing) && MOI.is_empty(m._input_model)

@enum(ObjectiveType, UNSET, SCALAR_AFFINE, NONLINEAR)

"""
$(TYPEDEF)

A structure used to expressions and problem descriptions EAGO uses to formulate
relaxed problems.
"""
Base.@kwdef mutable struct ParsedProblem{N,T} <: MOI.ModelLike

    # Problem classification (set in parse_classify_problem!)
    _problem_type::ProblemType = UNCLASSIFIED

    # objectives (set in initial_parse)
    _objective::SAF = SAF(SAT{T}[], zero(T))
    _objective_parsed::AffineFunction{T} = AffineFunction{T}()
    _optimization_sense::MOI.OptimizationSense = MOI.MIN_SENSE

    # non-single variable constraints (set in initial_parse)
    _saf_leq::Vector{AffineFunction{T}} = AffineFunction{T}[]
    _saf_eq::Vector{AffineFunction{T}} = AffineFunction{T}[]
    _sqf_leq::Vector{BufferedQuadraticIneq{T}} = BufferedQuadraticIneq{T}[]
    _sqf_eq::Vector{BufferedQuadraticEq{T}} = BufferedQuadraticEq{T}[]
    _conic_second_order::Vector{BufferedSOC{T}} = BufferedSOC{T}[]

    # nlp constraints (set in initial_parse)
    _nlp_data::Union{MOI.NLPBlockData,Nothing} = nothing

    # storage for nonlinear functions
    _nonlinear_constr::Vector{BufferedNonlinearFunction{T}} = BufferedNonlinearFunction{T}[]

    # nonlinear evaluator
    _relaxed_evaluator = Evaluator{N,T}()

    # variables (set in initial_parse)
    _variable_info::Vector{VariableInfo} = VariableInfo[]
    _variable_num::Int = 0
end

function Base.isempty(x::ParsedProblem)

    is_empty_flag = true

    new_input_problem = ParsedProblem()
    for field in fieldnames(ParsedProblem)

        field_value = getfield(x, field)
        if field_value isa Array
            if !isempty(field_value)
                is_empty_flag = false
                break
            end

        elseif field_value isa Number
            if getfield(new_input_problem, field) !== field_value
                is_empty_flag = false
                break
            end
        end
    end

    is_empty_flag &= x._nlp_data.evaluator isa EmptyNLPEvaluator
    is_empty_flag &= !x._nlp_data.has_objective
    is_empty_flag &= isempty(x._nlp_data.constraint_bounds)

    is_empty_flag &= isempty(x._objective.terms)
    is_empty_flag &= x._objective.constant == 0.0

    is_empty_flag &= isempty(x._objective.terms)
    is_empty_flag &= x._objective.constant == 0.0
    is_empty_flag &= isempty(x._objective_parsed.func.terms)
    is_empty_flag &= x._objective_parsed.func.constant == 0.0
    is_empty_flag &= x._objective_parsed.len == 0

    return is_empty_flag
end

#
# `GeometricMeanCone` is `SecondOrderCone` representable
# `RootDetConeTriangle` is representable by a `PositiveSemidefiniteConeTriangle` and an `GeometricMeanCone`
# `LogDetConeTriangle` is representable by a `PositiveSemidefiniteConeTriangle` and `ExponentialCone`
# Transforms a `G`-in-`GreaterThan{T}` constraint into an `F`-in-`LessThan{T}` drop support for GT
# IndicatorSOS1Bridge not need for
# The `NormSpectralCone` is representable with a PSD constraint, since
# ``t \\ge \\sigma_1(X)`` if and only if ``[tI X^\\top; X tI] \\succ 0``.
# The `RelativeEntropyCone` is representable with exponential cone and LP constraints
# The `RotatedSecondOrderCone` is `SecondOrderCone` representable
#=
The `geometric mean cone` is representable with a relative entropy constraint and a
nonnegative auxiliary variable
=#

#=
The `SemiToBinaryBridge` replaces an Semicontinuous constraint:
``x \\in \\mathsf{Semicontinuous}(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.
The `SemiToBinaryBridge` replaces an Semiinteger constraint:
``x \\in Semiinteger(l, u)``
is replaced by:
``z \\in \\{0, 1\\}``,
``x \\in \\mathbb{Z}``,
``x \\leq z \\cdot u ``,
``x \\geq z \\cdot l ``.
=#

#=
No need for Semicontinuous or Semiinteger wp representation... but may speed
up mixed integer solves.
=#

# LP relax bridges ->
# The `NormInfinityCone` is representable with LP constraints
# The `NormOneCone` is representable with LP constraint

# SlackBridge{T, F, G} to perform epigraph rearrangement -> FunctionizeBridge{T}
# only support scalar_affine_objectives


text:
using Test, CUDA, Knet, AutoGrad
using CUDA.CUDNN:
    cudnnSoftmaxForward,
    cudnnSoftmaxForward!,
    cudnnSoftmaxBackward,
    cudnnSoftmaxAlgorithm_t,
        CUDNN_SOFTMAX_FAST,     # 0, /* straightforward implementation */
        CUDNN_SOFTMAX_ACCURATE, # 1, /* subtract max from every point to avoid overflow */
        CUDNN_SOFTMAX_LOG,      # 2
    cudnnSoftmaxMode_t,
        CUDNN_SOFTMAX_MODE_INSTANCE, # 0, /* compute the softmax over all C, H, W for each N */
        CUDNN_SOFTMAX_MODE_CHANNEL,  # 1  /* compute the softmax over all C for each H, W, N */
    handle

@testset "cudnn/softmax" begin

    x = CUDA.randn(Float64,10,10)
    x2 = Param(x)
    x3 = Param(reshape(x, 10, 10, 1))
    @test @gcheck cudnnSoftmaxForward(x2, alpha=1, algo=CUDNN_SOFTMAX_FAST,     mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x2, alpha=2, algo=CUDNN_SOFTMAX_FAST,     mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x2, alpha=1, algo=CUDNN_SOFTMAX_ACCURATE, mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x2, alpha=2, algo=CUDNN_SOFTMAX_ACCURATE, mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x2, alpha=1, algo=CUDNN_SOFTMAX_LOG,      mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x2, alpha=2, algo=CUDNN_SOFTMAX_LOG,      mode=CUDNN_SOFTMAX_MODE_INSTANCE)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=1, algo=CUDNN_SOFTMAX_FAST,     mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=2, algo=CUDNN_SOFTMAX_FAST,     mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=1, algo=CUDNN_SOFTMAX_ACCURATE, mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=2, algo=CUDNN_SOFTMAX_ACCURATE, mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=1, algo=CUDNN_SOFTMAX_LOG,      mode=CUDNN_SOFTMAX_MODE_CHANNEL)
    @test @gcheck cudnnSoftmaxForward(x3, alpha=2, algo=CUDNN_SOFTMAX_LOG,      mode=CUDNN_SOFTMAX_MODE_CHANNEL)

end


text:
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Program: pi_monte_carlo.jl
#          Monte-Carlo calculation of PI
#
# Usage: julia pi_monte_carlo.jl
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function montepi(n::Int)
   R = 1.0
   s = 0
   for i = 1: n
      x = R * rand()
      y = R * rand()
      if x^2 + y^2 <= R^2
         s = s + 1
      end
   end
   return 4.0*s/n
end

# Main program
for i in 3: 8
    n = 10^i
    p = montepi(n)
    println("N = $n: PI = $p")
end

max_stars_repo_name: nikhilklath/julia-test

max_issues_repo_name: nikhilklath/julia-test

max_forks_repo_name: nikhilklath/julia-test


text:
export SeqFilter, IdFilter, idfilter

"""
    SeqFilter

An abstract type for the ensemble filtering algorithms.

"""
abstract type SeqFilter end

"""
    IdFilter <: SeqFilter

An immutable structure for the identity filter.

## Fields

- `Δtdyn::Float64`: time-step of the dynamical model
- `Δtobs::Float64`: time-step between two observations

## Constructors

- `IdFilter(Δtdyn, Δtobs)`: set up an identity filter
"""
struct IdFilter<:SeqFilter
	"Time step dynamic"
    Δtdyn::Float64

    "Time step observation"
    Δtobs::Float64
end

"""
	(filter::IdFilter)(X, ystar, t)

	Applies the identity transformation to the ensemble `X`.
"""
function (filter::IdFilter)(X, ystar, t)
	return X
end


text:
#-----------------------------------------------------------------------------------
# Helper;
#-----------------------------------------------------------------------------------
import ..Helper: Looping, current, next!, restart!,
         unzip, remove, is_line_number, linesof

"""
`Looping(iterable)` creates an objects that permits
iterating an infinite repetition of `iterable`.

`current(L::Looping)` gets the current element.

`next!(L)` makes `current(L)` return the next object.
if the current object is the last of the iterable,
start from the beginning.

`restart!(L)` makes `current(L)` return the first
object.

"""
Looping;

"""
`current(L::Looping)` returns the current element,
which can be changed through `next!` and `restart!`.

"""
current;

"""
Makes `current(L::Looping)` return the next object.
if the current object is the last of the iterable,
start from the beginning.

"""
next!;

"""
Makes `current(L::Looping)` return the first object.

"""
restart!;

"""
`unzip(Iterable{Tuple{F,S}}) -> Tuple{Vector{F}, Vector{S}}`

Takes a iterable of the form `[(x₁,y₁), (x₂,y₂), ..., (xₙ, yₙ)]` and
transforms it
in `([x₁, x₂, ..., xₙ], [y₁, y₂, ..., yₙ])`

"""
unzip;

"""
`remove((T -> Bool), Iterable{T}) -> Vector{T}`

Takes a function `f` and a collection and return that
collection with all elements `x` for which `f(x)` is
true removed.

"""
remove;

"""
`is_line_number(ex) -> Bool`

Returns true if `ex` represents a line number.
(`LineNumberNode` or `Expr(:line)`)

"""
is_line_number;

"""
`linesof(Expr) -> Expr`

Takes a block expression and return all non-linenumber
elements from it.

"""
linesof;


text:
using FileIO
using CSVFiles
using DataFrames

function save_datatable(data::DataFrame, file_path::AbstractString)
    return save_DataFrame(data, file_path;
        delimiter = space::Delimiter, header = true)
end

function save_DataFrame(
    data::DataFrame,
    file_path::AbstractString;
    delimiter::Delimiter = space::Delimiter,
    header::Bool = true)
    @assert(!isdir(file_path), ArgumentError("The file '$file_path' is a directory."))
    file_object = query(file_path)
    @assert(!unknown(file_object), ArgumentError("The file '$file_path' has unknown extension and/or format."))
    file_format = typeof(file_object).parameters[1]
    @info "Detected the format $file_format for the file '$file_path'. Saving to DataFrame..."

    data |> save(file_path;
        delim=DELIMITER[delimiter],
        nastring="NA", header=header)
    @info "... Data saved from a DataFrame to the file $file_path."
end


text:
export gui

function gui()
    ui = Builder(filename=joinpath(@__DIR__, "ui.glade"))
    showall(ui["win"])
    setproperty!(ui["imagebox"], :file, joinpath(@__DIR__, "Julia_prog_language_logo.svg"))
    visible(ui["forestfire_gif_button"], false)

    outputdir = joinpath(tempdir(), "JuliaPercolation")
    mkpath(outputdir)
    filename = ""

    function getparameters()
        linsize = parse(getproperty(ui["linent"], :text, String))
        nsample = parse(getproperty(ui["nsament"], :text, String))
        ps = parse(getproperty(ui["psent"], :text, String))
        pinc = parse(getproperty(ui["pincent"], :text, String))
        pf = parse(getproperty(ui["pfent"], :text, String))

        return linsize, nsample, ps, pinc, pf
    end

    function heatmap_button()
        linsize = parse(getproperty(ui["linent"], :text, String))
        p = parse(getproperty(ui["probent"], :text, String))
        latticetype = getproperty(ui["latticetype_cb"], "active", Int)
        if latticetype == 0
            l = Square(linsize, p)
        elseif latticetype == 1
            l = Triangular(linsize, p)
        elseif latticetype == 2
            l = Honeycomb(linsize, p)
        end
        heatmap!(l)
        filename = joinpath(outputdir, "cluster.png")
        savefig(filename)

        setproperty!(ui["imagebox"], :file, filename)
        println("Visualize $(typeof(l))")

        return nothing
    end

    function phase_transition_button()
        latticetype = getproperty(ui["latticetype_cb"], "active", Int)
        if latticetype in [0,1,2]
            plot_percolation_prob_button()
        elseif latticetype == 3
            plot_lifetime_button()
        end

        return nothing
    end

    function plot_percolation_prob_button()
        latticetypes = Dict(0=>"square", 1=>"triangular", 2=>"honeycomb")
        latticetype = latticetypes[getproperty(ui["latticetype_cb"], "active", Int)]
        linsize, nsample, ps, pinc, pf = getparameters()

        plot_percolation_prob(latticetype, linsize, ps, pinc, pf, nsample)
        filename = joinpath(outputdir, "percoprob.png")
        savefig(filename)

        setproperty!(ui["imagebox"], :file, filename)
        println("Plot Percolation Probability of $(latticetype) lattice")

        return nothing
    end

    function plot_lifetime_button()
        linsize, nsample, ps, pinc, pf = getparameters()

        plot_lifetime(linsize, ps, pinc, pf, nsample)
        filename =  joinpath(outputdir, "lifetime.png")
        savefig(filename)

        setproperty!(ui["imagebox"], :file, filename)
        println("Plot Lifetime of forest")
        return nothing
    end

    function forestfire_gif_button()
        linsize = parse(getproperty(ui["linent"], :text, String))
        p = parse(getproperty(ui["probent"], :text, String))

        forest = Forest(linsize, p)
        filename = forestfire!(forest, gifanime=true)[2]

        setproperty!(ui["imagebox"], :file, filename)
        println("Genarate GIF animation")
    end

    function savefig_button()
        dstname = save_dialog("Save figure")
        if contains(filename, ".png")
            if !contains(dstname, ".png")
                dstname *= ".png"
            end
        elseif contains(filename, ".gif")
            if !contains(dstname, ".gif")
                dstname *= ".gif"
            end
        end
        cp(filename, dstname, remove_destination=true)

        return nothing
    end

    signal_connect(x -> savefig_button(), ui["savefig_button"], "clicked")
    signal_connect(x -> heatmap_button(), ui["visualize_button"], "clicked")
    signal_connect(x -> phase_transition_button(), ui["phase_transition_button"], "clicked")
    signal_connect(x -> forestfire_gif_button(), ui["forestfire_gif_button"], "clicked")
    signal_connect(ui["latticetype_cb"], "changed") do widget, others...
        idx = getproperty(ui["latticetype_cb"], "active", Int)
        if idx in [0, 1, 2]
            setproperty!(ui["phase_transition_button"], :label, "Plot percolation probability")
            visible(ui["visualize_button"], true)
            visible(ui["forestfire_gif_button"], false)
        elseif idx == 3
            setproperty!(ui["phase_transition_button"], :label, "Plot Lifetime")
            visible(ui["visualize_button"], false)
            visible(ui["forestfire_gif_button"], true)
        end
    end

    if !isinteractive()
        c = Condition()
        signal_connect(ui["win"], :destroy) do widget
            notify(c)
        end
        wait(c)
    end

    return nothing
end


text:
using PyCall
using FlightSims
using FSimROS
using Plots
theme(:lime)
using FSimPlots
using UnPack
using OnlineStats

rclpy = pyimport("rclpy")
rosNode = pyimport("rclpy.node")
fsim_msg = pyimport("fsim_interfaces.msg")
std_msg = pyimport("std_msgs.msg")

@pydef mutable struct StateNode <: rosNode.Node
    function __init__(self)
        rosNode.Node.__init__(self, "viz_node")
        self.i = 0
        self.t = nothing
        self.time_received = false
        self.fig_multicopter = plot()  # it might not be a good idea to run plotting function in StateNode for latency
        self.fig_stat = plot()
        self.stat_E = MovingWindow(Float64, 500)
        self.stat_N = MovingWindow(Float64, 500)
        self.stat_U = MovingWindow(Float64, 500)
        # env
        self.multicopter = LeeHexacopter()
        # subscriber
        self.control_received = false
        function listener_callback(self, msg_state)
            # if self.t != nothing
                x = msg_to_state(self.multicopter, msg_state)
                self.i = self.i + 1
                # if false
                # plotting multicopter
                self.fig_multicopter = plot()
                plot!(self.fig_multicopter, self.multicopter, x; xlim=(-2, 2), ylim=(-2, 2), zlim=(-1, 10))
                # plotting stats
                fit!(self.stat_E, x.p[2])
                fit!(self.stat_N, x.p[1])
                fit!(self.stat_U, -x.p[3])
                fig_stat_E = plot(1:minimum([self.stat_E.b, self.stat_E.n]), self.stat_E.value; title="E [m]", label=nothing)
                fig_stat_N = plot(1:minimum([self.stat_N.b, self.stat_N.n]), self.stat_N.value; title="N [m]", label=nothing)
                fig_stat_U = plot(1:minimum([self.stat_U.b, self.stat_U.n]), self.stat_U.value; title="U [m]", label=nothing)
                self.fig_stat = plot(fig_stat_E, fig_stat_N, fig_stat_U; layout=(3, 1))
            # end
            display(plot(self.fig_multicopter, self.fig_stat, layout=(1, 2)))
        end
        function listener_callback_time(self, msg)
            t = msg.data
            if !self.time_received
                self.time_received = true
            end
            self.t = t
        end
        self.subscription = self.create_subscription(fsim_msg.PoseTwist, "state", msg -> listener_callback(self, msg), 10)
        self.subscription_time = self.create_subscription(std_msg.Float64, "time", msg -> listener_callback_time(self, msg), 10)
    end
end

function main(args=nothing)
    rclpy.init(args=args)
    println("viz node started")
    state_node = StateNode()
    rclpy.spin(state_node)
    state_node.destroy_node()
    rclpy.shutdown()
end

main()


text:
module NewPackage

using StaticArrays

greet() = print("Hello World!")

include("foo.jl")

end # module


text:
import Base: rand,
             ∈,
             isempty

export EmptySet, ∅,
       an_element,
       linear_map

"""
    EmptySet{N<:Real} <: LazySet{N}

Type that represents the empty set, i.e., the set with no elements.
"""
struct EmptySet{N<:Real} <: LazySet{N} end

# default constructor of type Float64
EmptySet() = EmptySet{Float64}()

"""
    ∅

An `EmptySet` instance of type `Float64`.
"""
const ∅ = EmptySet{Float64}()

# --- LazySet interface functions ---

"""
    dim(∅::EmptySet)

Return the dimension of the empty set, which is -1 by convention.

### Input

- `∅` -- an empty set

### Output

`-1` by convention.
"""
function dim(∅::EmptySet)::Int
    return -1
end

"""
    σ(d::AbstractVector{N}, ∅::EmptySet{N}) where {N<:Real}

Return the support vector of an empty set.

### Input

- `d` -- direction
- `∅` -- an empty set

### Output

An error.
"""
function σ(d::AbstractVector{N}, ∅::EmptySet{N}) where {N<:Real}
    error("the support vector of an empty set does not exist")
end

"""
    isbounded(∅::EmptySet)::Bool

Determine whether an empty set is bounded.

### Input

- `∅` -- empty set

### Output

`true`.
"""
function isbounded(::EmptySet)::Bool
    return true
end

"""
    ∈(x::AbstractVector{N}, ∅::EmptySet{N})::Bool where {N<:Real}

Check whether a given point is contained in an empty set.

### Input

- `x` -- point/vector
- `∅` -- empty set

### Output

The output is always `false`.

### Examples

```jldoctest
julia> ∈([1.0, 0.0], ∅)
false
```
"""
function ∈(x::AbstractVector{N}, ∅::EmptySet{N})::Bool where {N<:Real}
    return false
end

"""
    an_element(∅::EmptySet)

Return some element of an empty set.

### Input

- `∅` -- empty set

### Output

An error.
"""
function an_element(∅::EmptySet)
    error("an empty set does not have any element")
end

"""
    rand(::Type{EmptySet}; [N]::Type{<:Real}=Float64, [dim]::Int=0,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::EmptySet{N}

Create an empty set (note that there is nothing to randomize).

### Input

- `EmptySet` -- type for dispatch
- `N`        -- (optional, default: `Float64`) numeric type
- `dim`      -- (optional, default: 0) dimension (is ignored)
- `rng`      -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`     -- (optional, default: `nothing`) seed for reseeding

### Output

The (only) empty set of the given numeric type.
"""
function rand(::Type{EmptySet};
              N::Type{<:Real}=Float64,
              dim::Int=0,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing
             )::EmptySet{N}
    rng = reseed(rng, seed)
    return EmptySet{N}()
end

"""
    isempty(∅::EmptySet)::Bool

Return if the empty set is empty or not.

### Input

- `∅` -- empty set

### Output

`true`.
"""
function isempty(∅::EmptySet)::Bool
    return true
end

"""
    norm(S::EmptySet, [p]::Real=Inf)

Return the norm of an empty set.
It is the norm of the enclosing ball (of the given ``p``-norm) of minimal volume
that is centered in the origin.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function norm(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a norm")
end

"""
    radius(S::EmptySet, [p]::Real=Inf)

Return the radius of an empty set.
It is the radius of the enclosing ball (of the given ``p``-norm) of minimal
volume with the same center.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function radius(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a radius")
end

"""
    diameter(S::EmptySet, [p]::Real=Inf)

Return the diameter of an empty set.
It is the maximum distance between any two elements of the set, or,
equivalently, the diameter of the enclosing ball (of the given ``p``-norm) of
minimal volume with the same center.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function diameter(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a diameter")
end

"""
    linear_map(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N}

Return the linear map of an empty set.

### Input

- `M` -- matrix
- `∅` -- empty set

### Output

The empty set.
"""
linear_map(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N} = ∅

"""
    translate(∅::EmptySet{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) an empty set by a given vector.

### Input

- `∅` -- empty set
- `v` -- translation vector

### Output

The empty set.
"""
function translate(∅::EmptySet{N}, v::AbstractVector{N}) where {N<:Real}
    return ∅
end

"""
    plot_recipe(∅::EmptySet{N}, [ε]::N=zero(N)) where {N<:Real}

Convert an empty set to a sequence of points for plotting.
In the special case of an empty set, we define the sequence as `nothing`.

### Input

- `∅` -- empty set
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

`nothing`.
"""
function plot_recipe(∅::EmptySet{N}, ε::N=zero(N)) where {N<:Real}
    return nothing
end


text:
using Documenter
using SymbolicNumericIntegration

makedocs(
    sitename="SymbolicNumericIntegration",
    format=Documenter.HTML(),
    modules=[SymbolicNumericIntegration]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo="https://github.com/SciML/SymbolicNumericIntegration.jl";
    push_preview=true
)


text:
using Clang
using Test

@testset "struct_enum_field" begin
    # parse file
    trans_unit = parse_header(joinpath(@__DIR__, "c", "struct_enum_field.h"))
    GC.@preserve trans_unit begin
        # get root cursor
        root_cursor = getcursor(trans_unit)
        cursors = children(root_cursor)
        # enum
        enum_cursor = cursors[1]
        @test name(enum_cursor) == "Fruit"
        # struct
        struct_cursor = cursors[2]
        @test name(struct_cursor) == "Foo"
        # field
        struct_fields = children(struct_cursor)
        enum_field_cursor = struct_fields[1]
        int_field_cursor = struct_fields[2]
        array_field_cursor = struct_fields[3]
        @test clang2julia(enum_field_cursor) == :Fruit
        @test clang2julia(int_field_cursor) == :Cint
        @test clang2julia(array_field_cursor) == :(NTuple{3, Ptr{Cvoid}})
    end
end

@testset "struct_nested" begin
    # parse file
    trans_unit = parse_header(joinpath(@__DIR__, "c", "struct_nested.h"))
    ctx = DefaultContext()
    push!(ctx.trans_units, trans_unit)
    # get root cursor
    root_cursor = getcursor(trans_unit)
    cursors = children(root_cursor)
    # case 1
    wrap!(ctx, cursors[1])
    expr = :(struct ANONYMOUS1_union_in_struct
                x::Cstring
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:ANONYMOUS1_union_in_struct].items[1] == expr
    expr = :(struct nested_struct
                union_in_struct::ANONYMOUS1_union_in_struct
                struct_in_struct::ANONYMOUS2_struct_in_struct
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:nested_struct].items[1] == expr
    # case 2
    ctx.anonymous_counter = 0
    wrap!(ctx, cursors[3])
    expr = :(struct ANONYMOUS2_nest2
                x::Cint
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:ANONYMOUS2_nest2].items[1] == expr
    expr = :(struct ANONYMOUS1_nest1
                nest2::ANONYMOUS2_nest2
                nest3::ANONYMOUS2_nest2
                s::UInt8
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:ANONYMOUS1_nest1].items[1] == expr
    expr = :(struct nested_struct2
                nest1::ANONYMOUS1_nest1
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:nested_struct2].items[1] == expr
end


text:
function set_indices_to!(mutee::AbstractArray{Int,2},indices::AbstractArray{Int,2},set_val::Integer)
    for y in 1:size(indices,1)
        mutee[indices[y,1], indices[y,2]] = set_val
    end
end

function towards(origin::Integer, destination::Integer)
    origin > destination && return origin - 1
    origin < destination && return origin + 1
    return origin
end

max_stars_repo_name: RatanRSur/robots-from-scratch

max_issues_repo_name: RatanRSur/robots-from-scratch

max_forks_repo_name: RatanRSur/robots-from-scratch


text:
"""
This function is exactly like `eval_element_derivative`, but it evaluates
the antierivative of the element instead.
"""
function eval_element_antiderivative(dict::Dictionary, idx, x)
    idxn = native_index(dict, idx)
    @boundscheck checkbounds(dict, idxn)
    unsafe_eval_element_antiderivative1(dict, idxn, x)
end

function eval_element_antiderivative(dict::Dictionary, idx::LinearIndex, x)
    @boundscheck checkbounds(dict, idx)
    unsafe_eval_element_antiderivative1(dict, native_index(dict, idx), x)
end

function eval_element_extension_antiderivative(dict::Dictionary, idx, x)
    @boundscheck checkbounds(dict, idx)
    unsafe_eval_element_antiderivative(dict, native_index(dict, idx), x)
end

function unsafe_eval_element_antiderivative1(dict::Dictionary{S,T}, idx, x) where {S,T}
    in_support(dict, idx, x) ? unsafe_eval_element_antiderivative(dict, idx, x) : zero(T)
end

unsafe_eval_element_antiderivative(dict::Dictionary, idx, x) =
    unsafe_eval_element_antiderivative(dict, native_index(dict, idx), x)


text:
import Plots.plot

function plot(s1::Schumaker, interval::Tuple{R,R} = (s1.IntStarts_[1], s1.IntStarts_[length(s1.IntStarts_)]); derivs::Bool = false, grid_len::Integer = 200, plot_options::NamedTuple = (label = "Spline",), deriv_plot_options::NamedTuple = (label = "Spline - 1st deriv",),
              deriv2_plot_options::NamedTuple = (label = "Spline - 2nd deriv",), plt = missing) where R<:Real
    grid = collect(range(interval[1], interval[2], length=grid_len))
    return plot(s1, grid;  derivs = derivs, plot_options = plot_options, deriv_plot_options = deriv_plot_options, deriv2_plot_options = deriv2_plot_options, plt = plt)
end

function plot(s1::Schumaker, grid::AbstractArray{R,1}; derivs::Bool = false, plot_options::NamedTuple = (label = "Spline",), deriv_plot_options::NamedTuple = (label = "Spline - 1st deriv",),
              deriv2_plot_options::NamedTuple = (label = "Spline - 2nd deriv",), plt = missing) where R<:Real
    evals = s1.(grid)
    plt = ismissing(plt) ? plot(grid, evals; plot_options...) : plot!(plt, grid, evals; plot_options...)
    if derivs
        evals_1 = evaluate.(s1,grid, 1)
        plt = plot!(plt, grid, evals_1; deriv_plot_options...)
        evals_2 = evaluate.(s1,grid, 2)
        plt = plot!(plt, grid, evals_2; deriv2_plot_options...)
        return plt
    else
        return plt
    end
end

function plot(ss::AbstractArray{Schumaker,1}, interval::Tuple{R,R} = (ss[1].IntStarts_[1], ss[1].IntStarts_[length(ss[1].IntStarts_)]); derivs::Bool = false, grid_len::Integer = 200, plot_options::Union{AbstractArray{Tuple,1},Missing} = missing,
              deriv_plot_options::Union{AbstractArray{Tuple,1},Missing} = missing, deriv2_plot_options::Union{AbstractArray{Tuple,1},Missing} = missing, plt = missing)  where R<:Real
    grid = collect(range(interval[1], interval[2], length=grid_len))
    return plot(ss, grid; derivs = derivs, plot_options = plot_options, deriv_plot_options = deriv_plot_options, deriv2_plot_options = deriv2_plot_options, plt = plt)
end

function plot(ss::AbstractArray{Schumaker,1}, grid::AbstractArray{R,1}; derivs::Bool = false, plot_options::Union{AbstractArray{Tuple,1},Missing} = missing, deriv_plot_options::Union{AbstractArray{Tuple,1},Missing} = missing,
    deriv2_plot_options::Union{AbstractArray{Tuple,1},Missing} = missing, plt = missing) where R<:Real
    plt = ismissing(plt) ? plot() : plt
    for i in 1:length(ss)
        plot_options_i = ismissing(plot_options) ? (label = string("Spline ", i),) : plot_options[i]
        deriv_plot_options_i = ismissing(deriv_plot_options) ? (label = string("Spline ", i, " - 1st deriv"),) : deriv_plot_options[i]
        deriv2_plot_options_i = ismissing(deriv2_plot_options) ? (label = string("Spline ", i, " - 2nd deriv"),) : deriv2_plot_options[i]
        plt = plot(ss[i], grid; derivs = derivs, plot_options = plot_options_i, deriv_plot_options = deriv_plot_options_i, deriv2_plot_options = deriv2_plot_options_i, plt = plt)
    end
    return plt
end


text:
# This file is a part of EncodedArrays.jl, licensed under the MIT License (MIT).

using EncodedArrays
using Test

@testset "encoded_array" begin
    data = rand(Int16(-1000):Int16(2000), 21)
    codec = VarlenDiffArrayCodec()

    data_enc = data |> codec

    @testset "ctors and conversions" begin
        @test @inferred(Array(data_enc)) == data
        @test typeof(Array(data_enc)) == Array{eltype(data),1}
        @test @inferred(Array{Int16}(data_enc)) == data
        @test @inferred(Array{Int16,1}(data_enc)) == data
        @test @inferred(Vector(data_enc)) == data
        @test @inferred(Vector{Int16}(data_enc)) == data
        @test typeof(Vector{Int16}(data_enc)) == Vector{Int16}
        @test @inferred(Vector{Int32}(data_enc)) == data
        @test typeof(Vector{Int32}(data_enc)) == Vector{Int32}

        @test @inferred(convert(Array, data_enc)) == data
        @test typeof(convert(Array, data_enc)) == Array{eltype(data),1}
        @test @inferred(convert(Array{Int16}, data_enc)) == data
        @test @inferred(convert(Array{Int16,1}, data_enc)) == data
        @test @inferred(convert(Vector, data_enc)) == data
        @test @inferred(convert(Vector{Int16}, data_enc)) == data
        @test typeof(convert(Vector{Int16}, data_enc)) == Vector{Int16}
        @test @inferred(convert(Vector{Int32}, data_enc)) == data
        @test typeof(convert(Vector{Int32}, data_enc)) == Vector{Int32}
    end

    @testset "collect" begin
        @test @inferred(collect(data_enc)) == data
        @test typeof(collect(data_enc)) == typeof(data)
    end

    @testset "getindex" begin
        @test @inferred(data_enc[:]) == data
        @test typeof(data_enc[:]) == Array{eltype(data),1}

        @test @inferred(data_enc[:]) == data
        @test @inferred(data_enc[1:21]) == data
        @test @inferred(data_enc[5:15]) == data[5:15]
        @test @inferred(data_enc[7]) == data[7]
    end

    @testset "setindex!" begin
        tmp = zero.(data)
        @test (tmp[:] = data_enc) == data

        tmp = vcat(zero.(data), zero.(data))
        tmp2 = copy(tmp)
        tmp[10:30] = data_enc
        tmp2[10:30] = data
        @test tmp == tmp2
    end

    @testset "equality" begin
        @test @inferred data == data_enc
        @test @inferred data_enc == data
        @test @inferred data_enc == data_enc
    end

    @testset "append!" begin
        A = similar(data, 0)
        @test @inferred(append!(A, data_enc)) === A
        @test A == data

        A = data[1:4]
        @test @inferred(append!(A, data_enc)) === A
        @test A == vcat(data[1:4], data)
    end

    @testset "append!" begin
        A = similar(data)
        @test @inferred(copyto!(A, data_enc)) === A
        @test A == data

        @test_throws BoundsError @inferred(copyto!(similar(data, 5), data_enc))
    end
end # testset


text:
"""
    JumpDispersal <: SetCellRule

    JumpDispersal(; prob_threshold, spotrange)
    JumpDispersal{R}(; prob_threshold, spotrange)
    JumpDispersal{R,W}(; prob_threshold, spotrange)

Jump dispersal simulates a random long distance dispersal event, which occurs at a given
probablity. A random cell within the `spotrange` is invaded by duplicating the entire
population of the source cell at the destination cell.

# Keyword Arguments

- `prob_threshold`: probability of dispersal, between one and zero.
    May be a `Number`, an `Aux` array or another `Grid`:
- `spotrange`: number of cells in range of jumps, in any direction.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct JumpDispersal{R,W,PT,SR} <: SetCellRule{R,W}
    prob_threshold::PT
    spotrange::SR
end
function JumpDispersal{R,W}(;
    prob_threshold=Param(0.1, bounds=(0.0, 1.0)),
    spotrange=Param(30.0, bounds=(0.0, 100.0)),
) where {R,W}
    JumpDispersal{R,W}(prob_threshold, spotrange)
end

@inline function applyrule!(data, rule::JumpDispersal{R,W}, N, I) where {R,W}
    # Ignore empty cells
    N > zero(N) || return nothing
    # Random dispersal events
    p = get(data, rule.prob_threshold, I...)
    rand() < p || return nothing

    # Randomly select spotting distance
    intspot = round(Int, rule.spotrange)
    rnge = -intspot:intspot
    dest, is_inbounds = inbounds(data, (rand(rnge), rand(rnge)) .+ I)

    # Update spotted cell if it's on the grid
    is_inbounds && @inbounds add!(data[W], N, dest...)

    return nothing
end


text:
using Documenter, PDESolver, ODLCommonTools

# need to use the physics modules here too
using AdvectionEquationMod
using EulerEquationMod
using SimpleODEMod
include("../test/TestSystem.jl")
using TestSystem

# some abbreviations
advec = "solver/advection"
euler = "solver/euler"
ode = "solver/simpleODE"

makedocs(
  format = :html,
  sitename = "PDESolver.jl",
  pages = Any["PDESolver Introduction" => "index.md"
              "PDESolver Concepts" => Any[
                    "Intro" => "concepts/intro.md"
                    "PUMI" => "concepts/pumi.md"
                    "SBP" => "concepts/sbp.md"
                   ]
              "Building PDESolver" => Any[
                    "build.md",
                    "Build Options" => "deps_readme.md"
                   ]
              "DOC To be broken up or organized" => Any[
                    "Code Interfaces" => "interfaces.md"
                    "Code Parallelization" => "parallel.md"
                   ]
              "Examples" => Any[
                    "Isentropic Vortex" => "examples/isentropic.md"
                    "Unsteady Vortex" => "examples/unsteady.md"
                   ]
              "Frontend" => Any[
                    "Introduction" => "pdesolver.md"
                    "PDESolver User Interface" => "pdesolver_user.md"
                    "PDESolver Physics Interface" => "pdesolver_physics.md"
                    "PDESolver Structure" => "pdesolver_structure.md"
                   ]
              "Invocation" => Any[
                    "Calling PDESolver" => "invocation/calling.md"
                    "Interactive Session (experimental)" => "invocation/interactive.md"
                   ]
              "Solver" => Any[
                  "solver/Readme.md"
                  "solver/misc.md"
                  "solver/SolverCommon.md"
                  "Advection" => Any[
                      "Introduction" => "$advec/advection.md"
                      "Datatypes" => "$advec/types.md"
                      "Volume Integrals" => "$advec/volume.md"
                      "Face Integrals" => "$advec/flux.md"
                      "Boundary Integrals" => "$advec/bc.md"
                      "Initial Condition" => "$advec/ic.md"
                      "Source Term" => "$advec/source.md"
                      "Common Functions" => "$advec/common.md"
                      "Adjoint" => "$advec/adjoint.md"
                      "Boundary Functional" => "$advec/boundary_functional.md"
                     ]
                  "Euler" => Any[
                      "Introduction" => "$euler/euler.md"
                      "Datatypes" => "$euler/types.md"
                      "Volume Integrals" => "$euler/volume.md"
                      "Volume Integrals Jacobian" => "$euler/volume_diff.md"
                      "Face Integrals" => "$euler/flux.md"
                      "Face Integrals Jacobian" => "$euler/flux_diff.md"
                      "Face Element Integrals" => "$euler/faceElementIntegrals.md"
                      "Boundary Integrals" => "$euler/bc.md"
                      "Boundary Integrals Jacobian" => "$euler/bc_diff.md"
                      "Initial Conditions" => "$euler/ic.md"
                      "Source Term" => "$euler/source.md"
                      "Common Functions" => "$euler/common.md"
                      "Conversion" => "$euler/conversion.md"
                      "Numerical Flux Functions" => "$euler/flux_functions.md"
                      "Numerical Flux Functions Jacobian" => "$euler/flux_functions_diff.md"
                      "Stabilization" => "$euler/stabilization.md"
                      "Adjoint" => "$euler/adjoint.md"
                      "Boundary Functional" => "$euler/boundary_functional.md"
                      "Misc" => "$euler/misc.md"
                      "Homotopy" => "$euler/homotopy.md"
                      "Homotopy Jacobian" => "$euler/homotopy_diff.md"
                      "Eigensystem" => "$euler/eigensystem.md"
                      "Startup" => "$euler/startup.md"
                     ]
                  "Simple ODE" => Any[
                      "Main" => "$ode/simpleODE.md"
                     ]

                 ]  # end Solver
              "Input" => Any[
                  "Introduction" => "input/input.md"
                  "Important Keys" => "input/keys.md"
                 ]
              "LinearSolvers" => Any[
                "Introduction" => "linearsolvers/linearsolvers.md"
                "Preconditioners" => "linearsolvers/pc.md"
                "Linear Operators" => "linearsolvers/lo.md"
                "Linear Solvers" => "linearsolvers/ls.md"
                ]  # end LinearSolvers
              "NonlinearSolvers" => Any[
                  "Introduction" => "NonlinearSolvers/nonlinearsolvers.md"
                  "Steady" => "NonlinearSolvers/steady.md"
                  "Unsteady" => Any[
                      "Intro" => "NonlinearSolvers/unsteady/intro.md"
                      "Runge-Kutta" => "NonlinearSolvers/unsteady/rk4.md"
                      "LSERK" => "NonlinearSolvers/unsteady/lserk.md"
                      "Crank-Nicolson" => "NonlinearSolvers/unsteady/cn.md"
                      "Crank-Nicolson: Unsteady Adjoint" => "NonlinearSolvers/unsteady/cn_uadj.md"
                    ]
                  "Newton's Method" => "NonlinearSolvers/newton.md"
                  "Jacobian Calculation" => "NonlinearSolvers/jacobian.md"
                  "Jacobian Freezing" => "NonlinearSolvers/jac_recalc.md"
                  "Residual Evalution" => "NonlinearSolvers/residual_evaluation.md"
                  "Matrix Interface" => "NonlinearSolvers/matrix.md"
                  "Newton Inner" => "NonlinearSolvers/newton_inner.md"
                 ]
              "Utils" => Any[
                  "Main" => "Utils/Utils.md"
                  "Parallel Constructs" => "Utils/parallel.md"
                  "Projections" => "Utils/projections.md"
                  "Logging" => "Utils/logging.md"
                  "Input/Output" => "Utils/io.md"
                  "Checkpointing" => "Utils/checkpoint.md"
                  "Misccellaneous" => "Utils/misc.md"
                 ]
                "Testing" => Any[
                    "Introduction" => "test/Testing.md"
                    "Local Testing" => "test/Readme.md"
                    "CI Testing" => "test/Travis.md"
                    "Test API" => "test/TestSystem.md"
                  ]
             ] # end Home
)  # end mkdocs


text:
# Note: To support testing against the latest Julia nightly we'll assume the serialization
# format is the same as the last version of Julia we've tested against.
const SERIALIZED_HEADER = if VERSION >= v"1.5"
    "7JL\n\x04\0\0"
elseif VERSION >= v"1.4"
    "7JL\t\x04\0\0"
elseif VERSION >= v"1.2"
    "7JL\b\x04\0\0"
elseif VERSION >= v"1.0"
    "7JL\a\x04\0\0"
else
    error("Julia versions earlier than 1.0 are unsupported")
end

# Declare a new `isnan` and `isinf` functions to avoid type piracy. These new function work
# with `Char` and `Period` types as well as other types.
# Note: A generic `isfinite` is declared in Intervals.
isnan(x) = (x != x)::Bool
isnan(x::Real) = Base.isnan(x)

isinf(x) = !isnan(x) && !isfinite(x)
isinf(x::Real) = Base.isinf(x)


text:
using ErlangTerm
using Test

@testset "ErlangTerm" begin
    @testset "Streams" begin
        buffer = IOBuffer()
        serialize(buffer, 1)
        @test take!(buffer) == UInt8[131, 97, 1]
    end
    @testset "Small Int" begin
        @test serialize(1) == UInt8[131, 97, 1]
        @test deserialize(serialize(1)) == 1
    end
    @testset "Large Int" begin
        @test serialize(256) == UInt8[131, 98, 0, 0, 1, 0]
        @test deserialize(serialize(256)) == 256
    end
    @testset "Small Big, Int64" begin
        i = Int64(0x100000000)
        @test serialize(i) == UInt8[131, 110, 5, 0, 0, 0, 0, 0, 1]
        @test deserialize(serialize(i)) == i
        @test typeof(deserialize(serialize(i))) == Int64
        @test serialize(-i) == UInt8[131, 110, 5, 1, 0, 0, 0, 0, 1]
        @test deserialize(serialize(-i)) == -i
        @test typeof(deserialize(serialize(-i))) == Int64
    end
    @testset "Small Big, Int128" begin
        i = Int128(0x10000000000000000)
    @test serialize(i) == UInt8[131, 110, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
        @test deserialize(serialize(i)) == i
        @test typeof(deserialize(serialize(i))) == Int128
        @test serialize(-i) == UInt8[131, 110, 9, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1]
        @test deserialize(serialize(-i)) == -i
        @test typeof(deserialize(serialize(-i))) == Int128
    end
    @testset "Small Big, BigInt" begin
        i = BigInt(16)^99     # 100 byte integer
        @test serialize(i)[1:5] == UInt8[131, 110, 50, 0, 0]
        @test deserialize(serialize(i)) == i
        @test typeof(deserialize(serialize(i))) == BigInt
        @test serialize(-i)[1:5] == UInt8[131, 110, 50, 1, 0]
        @test deserialize(serialize(-i)) == -i
        @test typeof(deserialize(serialize(-i))) == BigInt
    end
    @testset "Large Big, BigInt" begin
        i = BigInt(16)^599    # 600 byte integer
        @test serialize(i)[1:8] == UInt8[131, 111, 0, 0, 1, 44, 0, 0]
        @test deserialize(serialize(i)) == i
        @test typeof(deserialize(serialize(i))) == BigInt
        @test serialize(-i)[1:8] == UInt8[131, 111, 0, 0, 1, 44, 1, 0]
        @test deserialize(serialize(-i)) == -i
        @test typeof(deserialize(serialize(-i))) == BigInt
    end
    @testset "Float" begin
        @test serialize(1.0) == UInt8[131, 70, 63, 240, 0, 0, 0, 0, 0, 0]
        @test deserialize(serialize(1.0)) == 1.0
    end
    @testset "Atom" begin
        @test deserialize(UInt8[131, 100, 0, 4, 97, 116, 111, 109]) == :atom
    end
    @testset "String" begin
        @test deserialize(UInt8[131, 107, 0, 2, 1, 2]) == [1, 2]
    end
    @testset "Atom UTF-8" begin
        @test serialize(:μ) == UInt8[131, 119, 2, 206, 188]
        @test deserialize(serialize(:μ)) == :μ
    end
    @testset "Large Atom UTF-8" begin
        sym = Symbol(repeat("μ", 256))
        @test deserialize(serialize(sym)) == sym
    end
    @testset "Binary" begin
        @test serialize("erlang") == UInt8[131, 109, 0, 0, 0, 6, 101, 114, 108, 97, 110, 103]
        @test deserialize(serialize("erlang")) == "erlang"
        @test serialize("αω") == UInt8[131, 109, 0, 0, 0, 4, 206, 177, 207, 137]
        @test deserialize(serialize("αω")) == "αω"
    end
    @testset "List" begin
        @test serialize([]) == UInt8[131, 106]
        @test deserialize(serialize([])) == []

        list = [1024, :μ, "erlang"]
        @test serialize(list) == UInt8[131, 108, 0, 0, 0, 3, 98, 0, 0, 4, 0, 119, 2, 206, 188, 109, 0,
                                       0, 0, 6, 101, 114, 108, 97, 110, 103, 106]
        @test deserialize(serialize(list)) == list
    end
    @testset "Map" begin
        dict = Dict(:α => 1, :β => 2)
        exp1 = UInt8[131, 116, 0, 0, 0, 2, 119, 2, 206, 177, 97, 1, 119, 2, 206, 178, 97, 2]
        exp2 = UInt8[131, 116, 0, 0, 0, 2, 119, 2, 206, 178, 97, 2, 119, 2, 206, 177, 97, 1]
        act = serialize(dict)
        @test act == exp1 || act == exp2
        @test deserialize(serialize(dict)) == dict
    end
    @testset "Tuple" begin
        small_tuple = (1, 2, 3)
        @test serialize(small_tuple) == UInt8[131, 104, 3, 97, 1, 97, 2, 97, 3]
        @test deserialize(serialize(small_tuple)) == small_tuple

        large_tuple = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
                       23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
                       43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62,
                       63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82,
                       83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101,
                       102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
                       118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133,
                       134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
                       150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165,
                       166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,
                       182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197,
                       198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213,
                       214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229,
                       230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245,
                       246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256)
        large_result = UInt8[131, 105, 0, 0, 1, 0, 97, 1, 97, 2, 97, 3, 97, 4, 97, 5, 97, 6, 97, 7, 97,
                             8, 97, 9, 97, 10, 97, 11, 97, 12, 97, 13, 97, 14, 97, 15, 97, 16, 97, 17,
                             97, 18, 97, 19, 97, 20, 97, 21, 97, 22, 97, 23, 97, 24, 97, 25, 97, 26, 97,
                             27, 97, 28, 97, 29, 97, 30, 97, 31, 97, 32, 97, 33, 97, 34, 97, 35, 97, 36,
                             97, 37, 97, 38, 97, 39, 97, 40, 97, 41, 97, 42, 97, 43, 97, 44, 97, 45, 97,
                             46, 97, 47, 97, 48, 97, 49, 97, 50, 97, 51, 97, 52, 97, 53, 97, 54, 97, 55,
                             97, 56, 97, 57, 97, 58, 97, 59, 97, 60, 97, 61, 97, 62, 97, 63, 97, 64, 97,
                             65, 97, 66, 97, 67, 97, 68, 97, 69, 97, 70, 97, 71, 97, 72, 97, 73, 97, 74,
                             97, 75, 97, 76, 97, 77, 97, 78, 97, 79, 97, 80, 97, 81, 97, 82, 97, 83, 97,
                             84, 97, 85, 97, 86, 97, 87, 97, 88, 97, 89, 97, 90, 97, 91, 97, 92, 97, 93,
                             97, 94, 97, 95, 97, 96, 97, 97, 97, 98, 97, 99, 97, 100, 97, 101, 97, 102,
                             97, 103, 97, 104, 97, 105, 97, 106, 97, 107, 97, 108, 97, 109, 97, 110, 97,
                             111, 97, 112, 97, 113, 97, 114, 97, 115, 97, 116, 97, 117, 97, 118, 97, 119,
                             97, 120, 97, 121, 97, 122, 97, 123, 97, 124, 97, 125, 97, 126, 97, 127, 97,
                             128, 97, 129, 97, 130, 97, 131, 97, 132, 97, 133, 97, 134, 97, 135, 97, 136,
                             97, 137, 97, 138, 97, 139, 97, 140, 97, 141, 97, 142, 97, 143, 97, 144, 97,
                             145, 97, 146, 97, 147, 97, 148, 97, 149, 97, 150, 97, 151, 97, 152, 97, 153,
                             97, 154, 97, 155, 97, 156, 97, 157, 97, 158, 97, 159, 97, 160, 97, 161, 97,
                             162, 97, 163, 97, 164, 97, 165, 97, 166, 97, 167, 97, 168, 97, 169, 97, 170,
                             97, 171, 97, 172, 97, 173, 97, 174, 97, 175, 97, 176, 97, 177, 97, 178, 97,
                             179, 97, 180, 97, 181, 97, 182, 97, 183, 97, 184, 97, 185, 97, 186, 97, 187,
                             97, 188, 97, 189, 97, 190, 97, 191, 97, 192, 97, 193, 97, 194, 97, 195, 97,
                             196, 97, 197, 97, 198, 97, 199, 97, 200, 97, 201, 97, 202, 97, 203, 97, 204,
                             97, 205, 97, 206, 97, 207, 97, 208, 97, 209, 97, 210, 97, 211, 97, 212, 97,
                             213, 97, 214, 97, 215, 97, 216, 97, 217, 97, 218, 97, 219, 97, 220, 97, 221,
                             97, 222, 97, 223, 97, 224, 97, 225, 97, 226, 97, 227, 97, 228, 97, 229, 97,
                             230, 97, 231, 97, 232, 97, 233, 97, 234, 97, 235, 97, 236, 97, 237, 97, 238,
                             97, 239, 97, 240, 97, 241, 97, 242, 97, 243, 97, 244, 97, 245, 97, 246, 97,
                             247, 97, 248, 97, 249, 97, 250, 97, 251, 97, 252, 97, 253, 97, 254, 97, 255,
                             98, 0, 0, 1, 0]
        @test serialize(large_tuple) == large_result
        @test deserialize(serialize(large_tuple)) == large_tuple
    end
    @testset "Nothing" begin
        @test serialize(nothing) == UInt8[131, 106]
    end
end