text:
######################################################################
# montgomerypoints.jl: xz-points for Montgomery curves
######################################################################

export XZPoint, xadd, xdouble, xinfinity, fixedtorsion, isfixedtorsion, XZzero

######################################################################
# Type definitions
######################################################################

"""
Concrete type for xz-points on Montgomery curves.
"""
mutable struct XZPoint{T<:Nemo.RingElem} <: ProjectivePoint{T}
    X::T
    Z::T
    curve::Montgomery{T}
end

base_ring(P::XZPoint) = Nemo.parent(P.X)

base_curve(P::XZPoint) = P.curve

"""
Get a description of an xz-point.
"""
show(io::IO, P::XZPoint) = print(io, "($(P.X):$(P.Z))")

"""
Create an xz-point from a projective point with 3 coordinates.
"""
function XZPoint(P::EllipticPoint{T}) where T<:Nemo.RingElem
    return P.Z == 0 ? xinfinity(P.curve) : XZPoint(P.X, P.Z, P.curve)
end

"""
Get a normalized xz-point from any xz-point.
"""
function normalized(P::XZPoint{T}) where T
    K = Nemo.parent(P.X)
    if isinfinity(P)
        return XZPoint(Nemo.one(K),
                               Nemo.zero(K),
                               P.curve)
    else
        return XZPoint(P.X // P.Z,
                               Nemo.one(K),
                               P.curve)
    end
end

"""
Normalizes a xz-point to a xz-point with Z-coordinate 1, in place.
"""
function normalize!(P::XZPoint{T}) where T<:Nemo.FieldElem
    K = Nemo.parent(P.X)
    if isinfinity(P)
        P.X = Nemo.one(K)
    else
        P.X = P.X // P.Z
        P.Z = Nemo.one(K)
    end
    return
end

"""
Decide whether two xz-points are given by the exact same coordinates.
"""
function samefields(P::XZPoint, Q::XZPoint)
	return (P.curve == Q.curve) & (P.X == Q.X) & (P.Z == Q.Z)
end

"""
Decide whether two xz-points are equal.
"""
==(P::XZPoint, Q::XZPoint) = P.X * Q.Z - P.Z * Q.X == 0

"""
Decide whether an xz-point is the point at infinity.
"""
isinfinity(P::XZPoint) = iszero(P.Z)

"""
Decide whether an xz-point is the fixed 2-torsion point of the Montgomery model (at the origin).
"""
isfixedtorsion(P::XZPoint) = iszero(P.X)

"""
Get the xz-point at infinity on a Montgomery curve.
"""
function xinfinity(E::Montgomery)
    K = Nemo.parent(E.A)
    return XZPoint(Nemo.one(K), Nemo.zero(K), E)
end

function XZzero(E::Montgomery)
	return xinfinity(E)
end

"""
Get the fixed xz- 2-torsion point on a Montgomery curve.
"""
function fixedtorsion(E::Montgomery)
    K = Nemo.parent(E.A)
    return XZPoint(Nemo.zero(K), Nemo.one(K), E)
end

function isvalid(P::XZPoint)
    return !(P.X == P.Z == 0) && issquare( P.curve.B * ( P.X^2 + P.curve.A*P.X*P.Z + P.Z^2 ) * P.X * P.Z )[1]
end

######################################################################
# Montgomery arithmetic
######################################################################

"""
Double any xz-point using the least possible field operations.
"""
function xdouble(P::XZPoint{T}) where T
    E = P.curve
    R = base_ring(E)

    v1 = P.X + P.Z
    v1 = v1^2
    v2 = P.X - P.Z
    v2 = v2^2

    X2 = v1 * v2

    v1 = v1 - v2
    v3 = ((E.A + 2) // 4) * v1
    v3 = v3 + v2

    Z2 = v1 * v3

    return XZPoint(X2, Z2, E)
end

"""
Differential addition on xz-points using the least possible field operations.

This function assumes the difference is not (0:0) or (0:1).
"""
function xadd(P::XZPoint{T}, Q::XZPoint{T}, Minus::XZPoint{T}) where T
    v0 = P.X + P.Z
    v1 = Q.X - Q.Z
    v1 = v1 * v0

    v0 = P.X - P.Z
    v2 = Q.X + Q.Z
    v2 = v2 * v0

    v3 = v1 + v2
    v3 = v3^2
    v4 = v1 - v2
    v4 = v4^2

    Xplus = Minus.Z * v3
    Zplus = Minus.X * v4

    return XZPoint(Xplus, Zplus, P.curve)
end

"""
Montgomery ladder to compute scalar multiplications of generic xz-points, using the least possible field operations.
"""
function xladder(k::Nemo.Integer, P::XZPoint)
    x0, x1 = P, xdouble(P)
    for b in Iterators.drop(Iterators.reverse(digits(Int8, k, base=2)), 1)
        if (b == 0)
        	x0, x1 = xdouble(x0), xadd(x0, x1, P)
        else
        	x0, x1 = xadd(x0, x1, P), xdouble(x1)
        end
    end
    return x0
end

"""
    Top-level function for scalar multiplications with xz-points on Montgomery curves
    """
function *(k::Nemo.Integer, P::XZPoint)
    E = P.curve
    if k == 0
	return xinfinity(E)
    elseif k<0
	return (-k) * P
    else
	if isinfinity(P)
	    return xinfinity(E)
	elseif isfixedtorsion(P)
	    if k % 2 == 0
		return xinfinity(E)
	    else
		return fixedtorsion(E)
	    end
	else
	    return xladder(k, P)
	end
    end
end

*(k::Nemo.fmpz, P::XZPoint) = BigInt(k) * P


text:
const term_c = (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE, Tokens.END, Tokens.ELSE, Tokens.ELSEIF, Tokens.CATCH, Tokens.FINALLY, Tokens.ENDMARKER)

function parse_block(ps::ParseState, ret::Vector{Any} = Any[], closers = (Tokens.END,), docable = false)
    while ps.nt.kind ∉ closers && !ps.errored
        if ps.nt.kind ∈ term_c
            if ps.nt.kind == Tokens.ENDMARKER
                break
            elseif ps.nt.kind == Tokens.RPAREN
                if length(ps.closer.cc) > 1 && :paren == ps.closer.cc[end-1]
                    break
                else
                    push!(ps.errors, Error((ps.nt.startbyte:ps.nt.endbyte) .+ 1 , "Unexpected )."))
                    push!(ret, ErrorToken(INSTANCE(next(ps))))
                end
            elseif ps.nt.kind == Tokens.RBRACE
                if length(ps.closer.cc) > 1 && :brace == ps.closer.cc[end-1]
                    break
                else
                    push!(ps.errors, Error((ps.nt.startbyte:ps.nt.endbyte) .+ 1 , "Unexpected }."))
                    push!(ret, ErrorToken(INSTANCE(next(ps))))
                end
            elseif ps.nt.kind == Tokens.RSQUARE
                if length(ps.closer.cc) > 1 && :square == ps.closer.cc[end-1]
                    break
                else
                    push!(ps.errors, Error((ps.nt.startbyte:ps.nt.endbyte) .+ 1 , "Unexpected ]."))
                    push!(ret, ErrorToken(INSTANCE(next(ps))))
                end
            else
                push!(ps.errors, Error((ps.nt.startbyte:ps.nt.endbyte) .+ 1 , "Unexpected $(ps.nt.kind)."))
                push!(ret, ErrorToken(INSTANCE(next(ps))))
            end
        else
            if docable
                a = parse_doc(ps)
            else
                a = parse_expression(ps)
            end
            push!(ret, a)
        end
    end
    return ret
end

function parse_iter(ps::ParseState)
    startbyte = ps.nt.startbyte
    if ps.nt.kind == Tokens.OUTER && ps.nws.kind != EmptyWS && !Tokens.isoperator(ps.nnt.kind)
        outer = INSTANCE(next(ps))
        arg = @closer ps range @closer ps ws parse_expression(ps)
        if is_range(arg)
            arg.arg1 = EXPR{Outer}([outer, arg.arg1])
            arg.fullspan += outer.fullspan
            arg.span = outer.fullspan + arg.span
        else
            arg = EXPR{ErrorToken}([outer, arg])
        end
    else
        arg = @closer ps range @closer ps ws parse_expression(ps)
    end
    return arg
end

function parse_ranges(ps::ParseState)
    startbyte = ps.nt.startbyte
    arg = parse_iter(ps)

    if (arg isa EXPR{Outer} && !is_range(arg.args[2])) || !is_range(arg)
        push!(ps.errors, Error(ps.nt.startbyte - arg.fullspan:ps.nt.startbyte , "Incorrect iteration specification."))
        arg = ErrorToken(arg)
    elseif ps.nt.kind == Tokens.COMMA
        arg = EXPR{Block}(Any[arg])
        while ps.nt.kind == Tokens.COMMA
            accept_comma(ps, arg)
            nextarg = parse_iter(ps)
            if (nextarg isa EXPR{Outer} && !is_range(nextarg.args[2])) || !is_range(nextarg)
                push!(ps.errors, Error(ps.nt.startbyte - nextarg.fullspan:ps.nt.startbyte , "Incorrect iteration specification."))
                arg = ErrorToken(arg)
            end
            push!(arg, nextarg)
        end
    end
    return arg
end

function is_range(x) false end
function is_range(x::BinarySyntaxOpCall) is_eq(x.op) end
function is_range(x::BinaryOpCall) is_in(x.op) || is_elof(x.op) end

function parse_end(ps::ParseState)
    if ps.closer.square
        ret = IDENTIFIER(ps)
    else
        push!(ps.errors, Error((ps.t.startbyte:ps.t.endbyte) .+ 1 , "Unexpected end."))
        ret = ErrorToken(IDENTIFIER(ps))
    end

    return ret
end

"""
    parse_call(ps, ret)

Parses a function call. Expects to start before the opening parentheses and is passed the expression declaring the function name, `ret`.
"""
function parse_call(ps::ParseState, ret, ismacro = false)
    sb = ps.nt.startbyte - ret.fullspan
    if is_minus(ret) || is_not(ret)
        arg = @closer ps unary @closer ps inwhere @precedence ps 13 parse_expression(ps)
        if arg isa EXPR{TupleH}
            pushfirst!(arg.args, ret)
            fullspan = ps.nt.startbyte - sb
            ret = EXPR{Call}(arg.args, fullspan, fullspan - (last(arg.args).fullspan - last(arg.args).span))
        elseif arg isa WhereOpCall && arg.arg1 isa EXPR{TupleH}
            ret = WhereOpCall(EXPR{Call}(Any[ret; arg.arg1.args]), arg.op, arg.args)
        else
            ret = UnaryOpCall(ret, arg)
        end
    elseif is_and(ret) || is_decl(ret) || is_exor(ret)
        arg = @precedence ps 20 parse_expression(ps)
        if is_exor(ret) && arg isa EXPR{TupleH} && length(arg.args) == 3 && arg.args[2] isa UnarySyntaxOpCall && is_dddot(arg.args[2].arg2)
            arg = EXPR{InvisBrackets}(arg.args)
        end
        ret = UnarySyntaxOpCall(ret, arg)
    elseif is_issubt(ret) || is_issupt(ret)
        arg = @precedence ps PowerOp parse_expression(ps)
        ret = EXPR{Call}(Any[ret; arg.args])
    else
        !ismacro && ret isa EXPR{MacroName} && (ismacro = true)
        args = Any[ret, PUNCTUATION(next(ps))]
        @closeparen ps @default ps parse_comma_sep(ps, args, !ismacro)
        accept_rparen(ps, args)
        fullspan = ps.nt.startbyte - sb
        ret = EXPR{ismacro ? MacroCall : Call}(args, fullspan, fullspan - last(args).fullspan + last(args).span)
    end
    return ret
end

function parse_comma_sep(ps::ParseState, args::Vector{Any}, kw = true, block = false, istuple = false)
    @nocloser ps inwhere @nocloser ps newline @closer ps comma while !closer(ps)
        a = parse_expression(ps)

        if kw && !ps.closer.brace && a isa BinarySyntaxOpCall && is_eq(a.op)
            a = EXPR{Kw}(Any[a.arg1, a.op, a.arg2], a.fullspan, a.span)
        end
        push!(args, a)
        if ps.nt.kind == Tokens.COMMA
            accept_comma(ps, args)
        end
        if ps.ws.kind == SemiColonWS
            break
        end
    end

    if ps.ws.kind == SemiColonWS
        if block && !(istuple && length(args) > 2) && !(length(args) == 1 && args[1] isa PUNCTUATION) && !(last(args) isa UnarySyntaxOpCall && is_dddot(last(args).arg2))
            args1 = Any[pop!(args)]
            @nocloser ps newline @closer ps comma while @nocloser ps semicolon !closer(ps)
                a = parse_expression(ps)
                push!(args1, a)
            end
            body = EXPR{Block}(args1)
            push!(args, body)
            args = body
        else
            parse_parameters(ps, args)
        end
    end
    return args
end

function parse_parameters(ps, args::Vector{Any})
    sb = ps.nt.startbyte
    args1 = Any[]
    @nocloser ps inwhere @nocloser ps newline  @closer ps comma while @nocloser ps semicolon !closer(ps)
        a = parse_expression(ps)
        if !ps.closer.brace && a isa BinarySyntaxOpCall && is_eq(a.op)
            a = EXPR{Kw}(Any[a.arg1, a.op, a.arg2], a.fullspan, a.span)
        end
        push!(args1, a)
        if ps.nt.kind == Tokens.COMMA
            accept_comma(ps, args1)
        end
        if ps.ws.kind == SemiColonWS
            parse_parameters(ps, args1)
        end
    end
    if !isempty(args1)
        fullspan = ps.nt.startbyte - sb
        paras = EXPR{Parameters}(args1, fullspan, fullspan - last(args1).fullspan + last(args1).span)
        push!(args, paras)
    end
    return
end

"""
    parse_macrocall(ps)

Parses a macro call. Expects to start on the `@`.
"""
function parse_macrocall(ps::ParseState)
    sb = ps.t.startbyte
    at = PUNCTUATION(ps)
    if !isemptyws(ps.ws)
        push!(ps.errors, Error((ps.ws.startbyte + 1:ps.ws.endbyte) .+ 1 , "Unexpected whitespace in macrocall."))
        mname = ErrorToken(INSTANCE(next(ps)))
    else
        mname = EXPR{MacroName}(Any[at, IDENTIFIER(next(ps))])
    end

    # Handle cases with @ at start of dotted expressions
    if ps.nt.kind == Tokens.DOT && isemptyws(ps.ws)
        while ps.nt.kind == Tokens.DOT
            op = OPERATOR(next(ps))
            nextarg = IDENTIFIER(next(ps))
            mname = BinarySyntaxOpCall(mname, op, Quotenode(nextarg))
        end
    end

    if ps.nt.kind == Tokens.COMMA
        return EXPR{MacroCall}(Any[mname], mname.fullspan, mname.span)
    elseif isemptyws(ps.ws) && ps.nt.kind == Tokens.LPAREN
        return parse_call(ps, mname, true)
    else
        args = Any[mname]
        insquare = ps.closer.insquare
        @default ps while !closer(ps)
            a = @closer ps inmacro @closer ps ws @closer ps wsop parse_expression(ps)
            push!(args, a)
            if insquare && ps.nt.kind == Tokens.FOR
                break
            end
        end
        # return EXPR{MacroCall}(args)
        fullspan = ps.nt.startbyte - sb
        return EXPR{MacroCall}(args, fullspan, fullspan - last(args).fullspan + last(args).span)
    end
end

"""
parse_generator(ps)

Having hit `for` not at the beginning of an expression return a generator.
Comprehensions are parsed as SQUAREs containing a generator.
"""
function parse_generator(ps::ParseState, @nospecialize ret)
    kw = KEYWORD(next(ps))
    ret = EXPR{Generator}(Any[ret, kw])
    ranges = @closesquare ps parse_ranges(ps)

    if ps.nt.kind == Tokens.IF
        if ranges isa EXPR{Block}
            ranges = EXPR{Filter}(ranges.args)
        else
            ranges = EXPR{Filter}(Any[ranges])
        end
        pushfirst!(ranges, KEYWORD(next(ps)))
        cond = @closer ps range parse_expression(ps)
        pushfirst!(ranges, cond)
        push!(ret, ranges)
    elseif ranges isa EXPR{Block}
        append!(ret, ranges)
    else
        push!(ret, ranges)
    end

    if ret.args[1] isa EXPR{Generator} || ret.args[1] isa EXPR{Flatten}
        ret = EXPR{Flatten}(Any[ret])
    end

    return ret
end

function parse_dot_mod(ps::ParseState, is_colon = false)
    args = Any[]

    while ps.nt.kind == Tokens.DOT || ps.nt.kind == Tokens.DDOT || ps.nt.kind == Tokens.DDDOT
        d = OPERATOR(next(ps))
        if is_dot(d)
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
        elseif is_ddot(d)
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
        elseif is_dddot(d)
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
            push!(args, OPERATOR(1, 1, Tokens.DOT, false))
        end
    end

    # import/export ..
    if ps.nt.kind == Tokens.COMMA || ps.ws.kind == NewLineWS || ps.nt.kind == Tokens.ENDMARKER
        if length(args) == 2
            return Any[INSTANCE(ps)]
        end
    end

    while true
        if ps.nt.kind == Tokens.AT_SIGN
            at = PUNCTUATION(next(ps))
            a = INSTANCE(next(ps))
            push!(args, EXPR{MacroName}(Any[at, a]))
        elseif ps.nt.kind == Tokens.LPAREN
            a = EXPR{InvisBrackets}(Any[PUNCTUATION(next(ps))])
            push!(a, @closeparen ps parse_expression(ps))
            accept_rparen(ps, a)
            push!(args, a)
        elseif ps.nt.kind == Tokens.EX_OR
            a = @closer ps comma parse_expression(ps)
            push!(args, a)
        elseif !is_colon && isoperator(ps.nt)
            next(ps)
            push!(args, OPERATOR(ps.nt.startbyte - ps.t.startbyte,  1 + ps.t.endbyte - ps.t.startbyte, ps.t.kind, false))
        else
            push!(args, INSTANCE(next(ps)))
        end

        if ps.nt.kind == Tokens.DOT
            push!(args, PUNCTUATION(next(ps)))
        elseif isoperator(ps.nt) && (ps.nt.dotop || ps.nt.kind == Tokens.DOT)
            push!(args, PUNCTUATION(Tokens.DOT, 1, 1))
            ps.nt = RawToken(ps.nt.kind, ps.nt.startpos, ps.nt.endpos, ps.nt.startbyte + 1, ps.nt.endbyte, ps.nt.token_error, false)
        else
            break
        end
    end
    args
end

ext: jl


text:
@testset "Updating bounds" begin
    m = MCPModel()
    @variable(m, x >= 0)
    @mapping(m, F, x+2)
    @complementarity(m, F, x)

    status = solveMCP(m)
    @test result_value(x) ≈ 0.

    JuMP.fix(x, 1., force=true)
    status = solveMCP(m)
    @test result_value(x) ≈ 1.

    JuMP.unfix(x)
    JuMP.set_lower_bound(x, 0)
    status = solveMCP(m)
    @test result_value(x) ≈ 0.
end


text:
add_assert_scope(:ModLattice)

@attributes mutable struct ModAlgAss{S, T, U}
  base_ring::S
  dim::Int
  isirreducible::Int # 0 not know
                     # 1 true
                     # 2 false
  degree_splitting_field::Int
                     # same as above
  algebra::U
  action_of_gens::Vector{T}
  action_of_basis::Vector{T}

  function ModAlgAss{T, U}(algebra::U; action_of_basis::Vector{T} = T[], action_of_gens::Vector{T} = T[]) where {T, U}
    S = typeof(base_ring(algebra))
    z = new{S, T, U}()
    if length(action_of_basis) > 0
      z.action_of_basis = action_of_basis
      z.dim = nrows(action_of_basis[1])
    end
    if length(action_of_gens) > 0
      z.action_of_gens = action_of_gens
      if !isdefined(z.dim)
        z.dim = nrows(action_of_gens[1])
      end
    end
    z.algebra = algebra
    z.base_ring = base_ring(algebra)
    if z.dim == 1
      z.isirreducible = 1
      z.degree_splitting_field = 1
    else
      z.isirreducible = 0
      z.degree_splitting_field = 0
    end
    return z
  end

  function ModAlgAss{T}(;action_of_gens::Vector{T} = T[]) where {T}
    K = base_ring(action_of_gens[1])
    U = matrix_algebra_type(K)
    S = typeof(K)
    z = new{S, T, U}()
    z.action_of_gens = action_of_gens
    z.base_ring = base_ring(action_of_gens[1])
    z.dim = nrows(action_of_gens[1])
    if z.dim == 1
      z.isirreducible = 1
      z.degree_splitting_field = 1
    else
      z.isirreducible = 0
      z.degree_splitting_field = 0
    end

    return z
  end
end

function Base.show(io::IO, V::ModAlgAss)
  print(io, "Module over field of dimension ", V.dim)
  if hasalgebra(V)
    print(io, " (with algebra defined))")
  end
end

################################################################################
#
#  Field access
#
################################################################################

function algebra(V::ModAlgAss)
  if !isdefined(V, :algebra)
    error("Algebra of module not defined")
  else
    return V.algebra
  end
end

Markdown.doc"""
    coefficient_ring(V::ModAlgAss) -> Field

Returns the coefficient ring of the module.
"""
coefficient_ring(V::ModAlgAss) = V.base_ring

function dim(V::ModAlgAss)
  return V.dim
end

@doc Markdown.doc"""
    Module(A::AbsAlgAss, M::Vector{<:MatElem})

Given an algebra $A$ over a field $K$ and a list of $\dim(A)$ of square
matrices over $K$, construct the $A$-module with `basis(A)[i]` acting
via `M[i]` from the right.
"""
function Module(A::AbsAlgAss, M::Vector{<:MatElem})
  @assert length(M) == length(basis(A))
  return ModAlgAss{eltype(M), typeof(A)}(A, action_of_basis = M)
end

@doc Markdown.doc"""
    Module(M::Vector{<:MatElem})

Given a list `M` of square matrices over a field $K$, construct the module
for the free algebra with the generators acting via `M[i]` from the right.

Note that the free algebra will not be constructed and the resulting
object has no assocated algebra.
"""
function Module(M::Vector{<:MatElem})
  return ModAlgAss{eltype(M)}(action_of_gens = M)
end

################################################################################
#
#  Predicates
#
################################################################################

@doc Markdown.doc"""
    hasalgebra(V::ModAlgAss) -> Bool

Returns whether the module was defined explicitely using an algebra.
"""
hasalgebra(V::ModAlgAss) = isdefined(V, :algebra)

@doc Markdown.doc"""
    hasmatrix_action(V::ModAlgAss) -> Bool

Returns whether the action on the module is given by matrices.
"""
function hasmatrix_action(V::ModAlgAss{S, T, U}) where {S, T, U}
  if T <: MatElem
    return true
  else
    return false
  end
end

################################################################################
#
#  Action
#
################################################################################

function action_of_basis(V::ModAlgAss)
  if isdefined(V, :action_of_basis)
    return V.action_of_basis
  else
    throw(NotImplemented())
  end
end

function action_of_basis(V::ModAlgAss, i::Int)
  if isdefined(V, :action_of_basis)
    return V.action_of_basis[i]
  else
    throw(NotImplemented())
  end
end

@doc Markdown.doc"""
    action(V::ModAlgAss{_, MatElem, _}, x::AbsAlgAssElem)

Given a $A$-module $V$ for an algebra $A$ with action given by matrices, and an
element $x$ of $A$, returns the action of $x$.
"""
function action(V::ModAlgAss{<:Any, <: MatElem, <:Any}, x::AbsAlgAssElem)
  @req parent(x) == algebra(V) "Algebra of module must be parent of element"
  cx = coefficients(x)
  M = cx[1] * action_of_basis(V, 1)
  for i in 2:length(cx)
    N = cx[i] * action_of_basis(V, i)
    M = add!(M, M, N)
  end
  return M
end

@doc Markdown.doc"""
    action(V::ModAlgAss) -> Vector

Given a module $V$, returns the action on $V$. If no algebra is defined,
these will be generators, otherwise these will be the action of `basis(A)`.
"""
function action(V::ModAlgAss)
  if !hasalgebra(V)
    return V.action_of_Gens
  else
    return action_of_basis(V)
  end
end

function action_of_order_basis(V::ModAlgAss{S, T, U}, O::AlgAssAbsOrd) where {S, T, U}
  s = get_attribute(V, :order_action)
  if s === nothing
    t = Dict{typeof(O), Vector{T}}()
    set_attribute!(V, :order_action => t)
  else
    t = s::Dict{typeof(O), Vector{T}}
    if haskey(t, O)
      return t[O]::Vector{T}
    end
  end
  z = T[action(V, elem_in_algebra(y)) for y in basis(O, copy = false)]
  t[O] = z
  return z
end

# Internal function to give the action in a consistent way
# If there is no algebra, just returns the action
# If there is an algebra, the minimal number of generators
# is returned.
function consistent_action(V::T, W::T) where {T <: ModAlgAss}
  if !hasalgebra(V)
    @assert !hasalgebra(W)
  else
    @assert hasalgebra(W)
    @assert algebra(V) === algebra(W)
  end

  if !hasalgebra(V)
    @assert length(V.action_of_gens) == length(W.action_of_gens)
    return (V.action_of_gens, W.action_of_gens)
  end

  if !isdefined(V.action_of_gens) || !isdefined(W.action_of_gens)
    return (V.action_of_basis, W.action_of_basis)
  else
    return (V.action_of_gens, W.action_of_gens)
  end
end

################################################################################
#
#  Irreducibility
#
################################################################################

function isirreducible_known(M::ModAlgAss)
  return M.isirreducible != 0
end

function isirreducible(M::ModAlgAss)
  if M.isirreducible != 0
    return M.isirreducible == 1
  else
    if !(coefficient_ring(M) isa FinField)
      error("Coefficient ring must be a finite field")
    end
    fl, N = meataxe(M)
    if fl
      M.isirreducible = 1
    else
      M.isirreducible = 2
    end
    return fl
  end
end

#function composition_factors(V::ModAlgAss{<:FinField})
#  act = V.action_of_gens
#
#  cf = stub_composition_factors(act)
#
#  return [Module(c) for c in cf]
#end

function basis_of_hom(V, W)
  x, y = consistent_action(V, W)
  return stub_basis_hom_space(x, y)
end

################################################################################
#
#  Function stubs, filled in by Oscar
#
################################################################################

stub_composition_factors(a) = error("Load Oscar (or GAP) and try again")

stub_basis_hom_space(a, b) = error("Load Oscar (or GAP) and try again")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using Base.Test
using TimerOutputs
const to = TimerOutput()

test_files = String[]
push!(test_files, "test_elements.jl")
push!(test_files, "test_problems.jl")
push!(test_files, "test_temperature.jl")

@testset "FEMTruss.jl" begin
    for fn in test_files
        timeit(to, fn) do
            include(fn)
        end
    end
end
println()
println("Test statistics:")
println(to)

  # Need some boundary conditions and forces
  # simple truss along x axis
  #b1 = Element(Poi1, [1])
  #update!(b1, "displacement 1", 0.0)
  #f1 = Element(Poi1, [2]) # Point load
  #update!(f1, "displacement traction force 1", 10.0)
  #solver = Solver(Linear, p1, f1, b1)
  #solver()
  # u = 10*10/(288*0.1)


text:
using LinearAlgebra, SphericalHarmonics, MultivariatePolynomials, Test

include("sphericalHarmonics.jl")
include("sphericalHarmonicsExpansion.jl")
include("sphericalQuadrature.jl")
include("fastfunc.jl")
include("translation.jl")
include("errorPropagation.jl")


text:
export timeseries

"""
    timeseries(notes::Notes, property, f [, grid]) -> tvec, ts
Produce a timeseries of the `property` of the given notes, optionally
first quantizing on the given `grid`.

After quantization, it is often
the case that many notes are in the same bin of the grid. The **function** `f`
denotes which value of the vector of the `property` of the notes to keep.
Typical values are `minimum, maximum, mean`, etc (*for type stability the returned
timeseries are always `Float64`*).

The `property` can be `:velocity`, `:pitch`, or `:duration`. Grid bins without any
notes are given the value `0`. This **can be problematic** if you request for
`:pitch` and your `notes` also include notes which actually have pitch `0`,
i.e. `C0`.
"""
function timeseries(notes, property, f, grid = 0:1//notes.tpq:1)
    isgrid(grid)
    tpq = notes.tpq
    qnotes = quantize(notes, grid)
    pos = positions(qnotes)
    # Real grid and grid bins
    rgrid = tpq .* grid
    bins = round.(Int, rgrid)[1:end-1]
    # tvec limits
    firstp = Int(qnotes[1].position)
    firstbin = findfirst(x -> x == mod(firstp, tpq), bins)

    # Create tvec
    tvec = (firstp÷tpq)*tpq .+ bins[firstbin:end]
    c = firstp÷tpq + 1
    while tvec[end] < notes[end].position
        append!(tvec, bins .+ c*tpq)
        c += 1
    end

    ts = zeros(Float64, size(tvec))

    i = previdx = 1; L = length(pos)
    while i ≤ L
        # find entries of same grid bin
        j = 1
        while j ≤ L - i && pos[i+j] == pos[i]
            j+=1
        end
        # get timeseries value for this note
        if j > 1
            d = [getfield(qnotes[k], property) for k in i:i+j-1]
            val = Float64(f(d))
        else
            val = Float64(getfield(qnotes[i], property))
        end

        # find pos[i] in tvec
        idx = findfirst(x -> x == pos[i], tvec)
        ts[idx] = val

        # TODO: For pitch, change zeros of intermediate entries
        # to pitch of previous notes

        i += j
    end
    return tvec, ts
end


text:
using Colors, Luxor

function compare_colors(color_a, color_b)
    # compare two colors, looking just at their LUV luminance values
    luv1 = convert(Luv, RGB(color_a[1]/255, color_a[2]/255, color_a[3]/255))
    luv2 = convert(Luv, RGB(color_b[1]/255, color_b[2]/255, color_b[3]/255))
    luv1.l > luv2.l
end

# sort colors into categories
# these are very bikesheddable
function classifycolornames(colornames)
    colorcategories = Dict{String, Array{String, 1}}(
        "whites" => [],
        "pinks" => [],
        "reds" => [],
        "oranges" => [],
        "yellows" => [],
        "greens" => ["lime"],
        "cyans" => [],
        "blues" => ["navy", "indigo", "teal"],
        "purples" => ["fuchsia"],
        "browns" => [],
        "grays" => ["black", "dimgrey", "lightgrey", "darkgrey"])
    for colorname in colornames
        if occursin(r"(^grey)|(^gray)", colorname)
            push!(colorcategories["grays"], colorname)
        elseif occursin(r"turquoise|azure|gains|blue|slate", colorname)
            push!(colorcategories["blues"], colorname)
        elseif occursin(r"lemon|khaki|wheat", colorname)
            push!(colorcategories["yellows"], colorname)
        elseif occursin(r"chartreuse|olive|aqua|green", colorname)
            push!(colorcategories["greens"], colorname)
        elseif occursin(r"mistyrose|plum|thistle|lavender|maroon|violet|orchid|magenta|maroon", colorname)
            push!(colorcategories["purples"], colorname)
        elseif occursin(r".*gold.*", colorname)
            push!(colorcategories["yellows"], colorname)
        elseif occursin(r"chocolate|sienna|wood|moccasin|bisque|peru|peach|papaya|almond|tan", colorname)
            push!(colorcategories["browns"], colorname)
        elseif occursin(r"crimson|firebrick|salmon|coral|tomato", colorname)
            push!(colorcategories["reds"], colorname)
        elseif occursin(r"honey|mint|beige|linen|corn|sea|silver|lace|ivory|snow|white", colorname)
            push!(colorcategories["whites"], colorname)
        else
            # some colors obviously belong to categories
            # ie "deeppink" is in "pinks"
            for k in keys(colorcategories)
                if occursin(replace(k, "s" => ""), colorname)
                    push!(colorcategories[k], colorname)
                end
            end
        end
    end
    # sort the RGB values in Luv.l format
    # TODO perhaps do this conversion just once instead of 1000s of times :)
    for k in keys(colorcategories)
        vals = colorcategories[k]
        sort!(vals, lt = (color1, color2) ->
            compare_colors(Colors.color_names[color1] ./ 255,
                Colors.color_names[color2] ./ 255))
        colorcategories[k] = unique(vals)
    end
    return colorcategories
end

"""
find a readable color for text placed on top of another color
"""
function inversecolor(foregroundcolor, backgroundcolor;
        tolerance = 20)
    foreground = foregroundcolor
    if colordiff(parse(Colorant, foreground), parse(Colorant, backgroundcolor)) < tolerance
        tempcolor = parse(Colorant, backgroundcolor)
        clab = convert(Lab, parse(Colorant, tempcolor))
        if clab.l > 50
            labelbrightness = 0
        else
            labelbrightness = 100
        end
        foreground = convert(RGB, Lab(labelbrightness, clab.b, clab.a))
    end
    return foreground
end

function main()
    colornames = collect(keys(Colors.color_names))
    colordictionary = classifycolornames(colornames)
    categoryorder = ["whites", "yellows", "greens", "cyans", "blues", "purples", "pinks", "reds", "oranges",  "browns", "grays"]
    checkclassified = String[]
    for category in categoryorder
        for c in colordictionary[category]
            push!(checkclassified, c)
        end
    end
    stilltodo = setdiff(colornames, checkclassified)
    length(stilltodo) > 0 && (@show stilltodo; @warn "gotta catch them all!")
    for colorcategory in categoryorder
        numbercells = length(keys(colordictionary[colorcategory]))
        numbercols = 10
        colwidth = 150
        rowheight = 100
        pagewidth = 800
        margin = 15
        numberrows = convert(Int, ceil(numbercells/numbercols) + 1)
        Drawing(pagewidth - margin, margin + (numberrows * rowheight),
            "/tmp/namedcolorchart-$(colorcategory).svg")
        origin()
        background("white")
        setcolor("black")
        fontsize(24)
        fontface("Helvetica-Bold")
        text(titlecase(colorcategory), BoundingBox()[1] + (margin, 30), halign=:left)
        setline(0.5)
        line(Point(-pagewidth/2 + margin, 40 + -currentdrawing.height/2), Point(pagewidth/2 - margin, 40 + -currentdrawing.height/2), :stroke)
        table = Table(numberrows, numbercols,
            (pagewidth - 2margin)/numbercols, # cell width
            rowheight, # row height
            O + (0, 3margin) # initial center position
            )
        pts = first.(collect(table))
        n = 1
        for col in colordictionary[colorcategory]
            backgroundcolor = Colors.RGB(Colors.color_names[col] ./ 255 ...)
            setcolor(backgroundcolor)
            # draw swatch
            box(pts[n], table.colwidths[1], table.rowheights[1]/2, :fill)
            # if overlapping text on swatch, choose color wisely
            # sethue(inversecolor("white", backgroundcolor, tolerance=60))
            sethue("black")
            fontface("DINNextLTPro-BoldCondensed")
            fontsize(12)
            text(col, pts[n] + (0, -table.rowheights[1]/3), halign=:center)
            fontsize(9)
            fontface("InputMonoCompressed-Regular")
            text(string(round(backgroundcolor.r, digits=2), " ",
                        round(backgroundcolor.g, digits=2), " ",
                        round(backgroundcolor.b, digits=2)),
                        pts[n] + (0, table.rowheights[1]/3), halign=:center)
            n += 1
        end
        finish()
    end
    return colordictionary
end

main()


text:
export reverse_geocode

# Remote file: cities1000.zip
const _cities1000r = @RemoteFile(
    "https://download.geonames.org/export/dump/cities1000.zip",
    file = "cities1000.zip",
    dir = "$path_rg_root",
    updates = :yearly
    )

# Remote file: countryInfo.txt
const _countryInfor = @RemoteFile(
    "https://download.geonames.org/export/dump/countryInfo.txt",
    file = "countryInfo.txt",
    dir = "$path_rg_root",
    updates = :yearly
    )

#=
a = """
geonameid         : integer id of record in geonames database
name              : name of geographical point (utf8) varchar(200)
asciiname         : name of geographical point in plain ascii characters, varchar(200)
alternatenames    : alternatenames, comma separated, ascii names automatically transliterated, convenience attribute from alternatename table, varchar(10000)
latitude          : latitude in decimal degrees (wgs84)
longitude         : longitude in decimal degrees (wgs84)
feature class     : see http://www.geonames.org/export/codes.html, char(1)
feature code      : see http://www.geonames.org/export/codes.html, varchar(10)
country code      : ISO-3166 2-letter country code, 2 characters
cc2               : alternate country codes, comma separated, ISO-3166 2-letter country code, 200 characters
admin1 code       : fipscode (subject to change to iso code), see exceptions below, see file admin1Codes.txt for display names of this code; varchar(20)
admin2 code       : code for the second administrative division, a county in the US, see file admin2Codes.txt; varchar(80)
admin3 code       : code for third level administrative division, varchar(20)
admin4 code       : code for fourth level administrative division, varchar(20)
population        : bigint (8 byte int)
elevation         : in meters, integer
dem               : digital elevation model, srtm3 or gtopo30, average elevation of 3''x3'' (ca 90mx90m) or 30''x30'' (ca 900mx900m) area in meters, integer. srtm processed by cgiar/ciat.
timezone          : the iana timezone id (see file timeZone.txt) varchar(40)
modification date : date of last modification in yyyy-MM-dd format
"""
lines = split(a,'\n')
for iline in lines
    println('"',strip(iline[1:18],),'"',',','"',strip(iline[20:end]),'"')
end =#

function read_cities(file::String)

    coordinates1 = Array{Float64,1}()
    locations = Array{Tuple{String,String,},1}()
    # locations1 = Array{CSV.Row,1}()
    suggest_size = 160_000

    sizehint!(coordinates1, suggest_size * 3)
    sizehint!(locations, suggest_size)
    # sizehint!(locations1, suggest_size)
    if isfile(file)
        for row in CSV.File(file, delim='\t',header=[:gid,:name,:ascname,:altname,
                                            :lat,:lon,:fcl,:fco,:cc,
                                            :cc2,:a1c,:a2c,:a3c,:a4c,
                                            :pop,:ele,:dem,:tz,:md])
            e = ECEFfromLLA(wgs84)(LLA(row.lat, row.lon))
            push!(coordinates1, e.x)
            push!(coordinates1, e.y)
            push!(coordinates1, e.z)
            push!(locations, (row.cc, row.name))
            # push!(locations1, row)
            # @show row
            # break
        end
    end
    reshape(coordinates1, 3, :), locations
end

#=
hs = "ISO	ISO3	ISO-Numeric	fips	Country	Capital	Area(in sq km)	Population	Continent	tld	CurrencyCode	CurrencyName	Phone	Postal Code Format	Postal Code Regex	Languages	geonameid	neighbours	EquivalentFipsCode"

a =#
function read_countries(file::String)
    ci = Dict{String,Tuple{String,Int64,String,String,Float64,Float64}}()
    if isfile(file)
        for row in CSV.File(file, delim='\t',  comment="#",
                                header=[:iso,:iso3,:ison,:fips,:name,:capital,:area,:pup,:cont,
                                            :c10,:c11, :c12, :c13, :c14,
                                            :c15, :c16, :c17, :c18, :c19,], )
            # @show file
            # @show row
            ci[row.iso] = row.iso3, row.ison,
                            row.name, ismissing(row.capital) ? "" : row.capital,
                            ismissing(row.area) ? NaN : row.area, row.pup
            # push!(locations1, row)

            # break
        end
    end
    ci
end

function load_geocode(cities1000r::RemoteFile=_cities1000r, countryInfor::RemoteFile=_countryInfor, )

    for ifile in [countryInfor,cities1000r] download(ifile) end

    cities1000 = path(cities1000r)
    if endswith(cities1000, ".zip")
        c = string(cities1000[1:end - 4], ".txt")
        bc = basename(c)
        if isfile(c) && mtime(cities1000) - mtime(c) < 86400
            cities1000 = c
        else
            r = ZipFile.Reader(cities1000);
            for f in r.files
                if f.name == bc
                    println("Filename: $(f.name)")
                    open(c, "w") do fo
                        write(fo, read(f, String));
                    end
                end
            end
            close(r)
            cities1000 = c
        end
    end
    # @show cities1000
    a, b = read_cities(cities1000);
    # @show path(countryInfor)
    c = read_countries(path(countryInfor))

    a, b, c
end

const coordinates, locations, countries =  load_geocode(_cities1000r, _countryInfor)
const kdtree = KDTree(coordinates)

function reverse_geocode(e1::ECEF, kdtree::KDTree=kdtree,
                    locations::Array{Tuple{String,String},1}=locations,
                    countries::Dict{String,Tuple{String,Int64,String,String,Float64,Float64}}=countries )
    idxs, dists = knn(kdtree, e1, 1, true)
    c, city = locations[idxs[1]]
    if haskey(countries, c)
        cc = countries[c]
    else
        cc = ("???", 0, "?", "?", NaN, NaN)
    end
    c, city, dists[1], cc
end


text:
function relaxed_initial_state(mdp::NoCrashProblem, steps=200,
                             rng=Base.GLOBAL_RNG;
                             solver=BehaviorSolver(NORMAL, true, rng))

    mdp = deepcopy(mdp)
    mdp.dmodel.max_dist = Inf
    mdp.dmodel.brake_terminate_thresh = Inf
    mdp.dmodel.lane_terminate = false
    pp = mdp.dmodel.phys_param
    is = MLState(0.0, 0.0, CarState[CarState(pp.lane_length/2, 1.0, pp.v_med, 0.0, NORMAL, 1)])
    sim = HistoryRecorder(max_steps=steps, rng=rng)
    policy = solve(solver, mdp)
    hist = POMDPs.simulate(sim, mdp, policy, is)
    s = last(state_hist(hist))
    s.t = 0.0
    s.x = 0.0
    @assert s.cars[1].y == 1.0
    @assert nothing === (s.terminal)
    return s
end


text:
# From https://chengchingwen.github.io/Transformers.jl/v0.1/tutorial/

# Copy task, point is to learn to copy one string to another from test data
# This is modified so that the sample data is random size between 1-20
using Flux
using Transformers
using Transformers.Basic #for loading the positional embedding

labels = collect(1:10)
startsym = 11
endsym = 12
unksym = 0
labels = [unksym, startsym, endsym, labels...]
vocab = Vocabulary(labels, unksym)

#function for generate training datas
sample_data() = (d = rand(1:10, rand(1:20)); (d,d))
#function for adding start & end symbol
preprocess(x) = [startsym, x..., endsym]   # add start and end sym to each data

@show sample = preprocess.(sample_data())   # apply the preprocess to all samples
@show encoded_sample = vocab(sample[1]) #use Vocabulary to encode the training data

sample = preprocess.(sample_data()) # = ([11, 5, 4, 2, 5, 2, 5, 5, 5, 7, 8, 12], [11, 5, 4, 2, 5, 2, 5, 5, 5, 7, 8, 12])
encoded_sample = vocab(sample[1]) # = [2, 8, 7, 5, 8, 5, 8, 8, 8, 10, 11, 3]

#define a Word embedding layer which turn word index to word vector
embed = Embed(512, length(vocab)) # |> gpu
#define a position embedding layer metioned above
pe = PositionEmbedding(512) # |> gpu

#wrapper for get embedding
function embedding(x)
  we = embed(x, inv(sqrt(512)))
  e = we .+ pe(we)
	return e
end

#define 2 layer of transformer
encode_t1 = Transformer(512, 8, 64, 2048) #|> gpu
encode_t2 = Transformer(512, 8, 64, 2048) #|> gpu

#define 2 layer of transformer decoder
decode_t1 = TransformerDecoder(512, 8, 64, 2048) #|> gpu
decode_t2 = TransformerDecoder(512, 8, 64, 2048) #|> gpu

#define the layer to get the final output probabilities
linear = Positionwise(Dense(512, length(vocab)), logsoftmax) #|> gpu

function encoder_forward(x)
  e = embedding(x)
  t1 = encode_t1(e)
  t2 = encode_t2(t1)
  return t2
end@sa

function decoder_forward(x, m)
  e = embedding(x)
  t1 = decode_t1(e, m)
  t2 = decode_t2(t1, m)
  p = linear(t2)
	return p
end

enc = encoder_forward(encoded_sample)
probs = decoder_forward(encoded_sample, enc)

function smooth(et)
    sm = fill!(similar(et, Float32), 1e-6/size(embed, 2))
    p = sm .* (1 .+ -et)
    label = p .+ et .* (1 - convert(Float32, 1e-6))
    label
end

#define loss function
function loss(x, y)
  label = onehot(vocab, y) #turn the index to one-hot encoding
  label = smooth(label) #perform label smoothing
  enc = encoder_forward(x)
	probs = decoder_forward(y, enc)
  l = logkldivergence(label[:, 2:end, :], probs[:, 1:end-1, :])
  return l
end

#collect all the parameters
ps = params(embed, pe, encode_t1, encode_t2, decode_t1, decode_t2, linear)
opt = ADAM(1e-4)

#function for created batched data
using Transformers.Datasets: batched

#flux function for update parameters
using Flux: gradient
using Flux.Optimise: update!

#define training loop
function train!()
  @info "start training"
  for i = 1:2000
    data = batched([sample_data() for i = 1:32]) #create 32 random sample and batched
		x, y = preprocess.(data[1]), preprocess.(data[2])
    x, y = vocab(x), vocab(y)#encode the data
    x, y = todevice(x, y) #move to gpu
    l = loss(x, y)
    grad = gradient(()->l, ps)
    if i % 8 == 0
    	println("loss = $l")
    end
    update!(opt, ps, grad)
  end
end

train!()

## Save parameters, here pieces are saved independently, not a good idea in general
using BSON: @save

@save "encode_t1.bson" encode_t1
@save "encode_t2.bson" encode_t2
@save "decode_t1.bson" decode_t1
@save "decode_t2.bson" decode_t2
@save "linear.bson" linear
######  RUN

using Flux: onecold

function translate(x)
    ix = todevice(vocab(preprocess(x)))
    seq = [startsym]

    enc = encoder_forward(ix)

    len = length(ix)
    for i = 1:2len
        trg = todevice(vocab(seq))
        dec = decoder_forward(trg, enc)
        #move back to gpu due to argmax wrong result on CuArrays
        ntok = onecold(collect(dec), labels)
        push!(seq, ntok[end])
        ntok[end] == endsym && break
    end
  seq[2:end-1]
end

translate([5,5,6,6,1,2,3,4,7, 10])

translate([2,3,7,8,8,9])


text:
function convolution1d!(tid, iframe, x, nx, rx, gridsize_x, grid_x, bandwidth, f_x, w)
    dx = x[iframe] - grid_x[1]
    ind_min = floor(Int, (dx - rx)/gridsize_x) + 1
    ind_min = ind_min >= 1 ? ind_min : 1
    ind_max = floor(Int, (dx + rx)/gridsize_x) + 2
    ind_max = ind_max <= nx ? ind_max : nx
    for ind in ind_min:ind_max
        dx = (grid_x[ind] - x[iframe])/bandwidth
        f_x[ind, tid] += w*exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth)
    end
end

function convolution1d_box!(tid, iframe, boxsize, x, nx, rx, gridsize_x, grid_x, bandwidth, f_x, w)
    for ind in 1:nx
        dx = x[iframe] - grid_x[ind]
        dx = dx - round(dx/boxsize)*boxsize;
        if abs(dx) < rx
          dx = dx/bandwidth
          f_x[ind, tid] += w*exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth)
        end
    end
end

function convolution2d!(tid, iframe,
    x, nx, rx, gridsize_x, grid_x,
    y, ny, ry, gridsize_y, grid_y,
    bandwidth, f_x, f_y, f_private, w)

    f_x[:, tid] .= 0.0
    f_y[:, tid] .= 0.0

    dx = x[iframe] - grid_x[1]
    ix_min = floor(Int, (dx - rx)/gridsize_x) + 1
    ix_min = ix_min >= 1 ? ix_min : 1
    ix_max = floor(Int, (dx + rx)/gridsize_x) + 2
    ix_max = ix_max <= nx ? ix_max : nx
    for ix in ix_min:ix_max
        dx = (grid_x[ix] - x[iframe])/bandwidth[1]
        f_x[ix, tid] += exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth[1])
    end

    dy = y[iframe] - grid_y[1]
    iy_min = floor(Int, (dy - ry)/gridsize_y) + 1
    iy_min = iy_min >= 1 ? iy_min : 1
    iy_max = floor(Int, (dy + ry)/gridsize_y) + 2
    iy_max = iy_max <= ny ? iy_max : ny
    for iy in iy_min:iy_max
        dy = (grid_y[iy] - y[iframe])/bandwidth[2]
        f_y[iy, tid] += exp(-0.5*dy*dy)/(sqrt(2.0*pi)*bandwidth[2])
    end

    for ix in ix_min:ix_max
        for iy in iy_min:iy_max
            f_private[ix, iy, tid] += w*f_x[ix, tid]*f_y[iy, tid]
        end
    end
end

function convolution2d_box!(tid, iframe, boxsize,
    x, nx, rx, gridsize_x, grid_x,
    y, ny, ry, gridsize_y, grid_y,
    bandwidth, f_x, f_y, f_private, w,
    ix_array, iy_array)

    f_x[:, tid] .= 0.0
    f_y[:, tid] .= 0.0

    ix_count = 0
    for ix in 1:nx
        dx = x[iframe] - grid_x[ix]
        dx = dx - round(dx/boxsize[1])*boxsize[1];
        if abs(dx) < rx
          dx = dx/bandwidth[1]
          f_x[ix, tid] += exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth[1])
          ix_count += 1
          ix_array[ix_count, tid] = ix
        end
    end

    iy_count = 0
    for iy in 1:ny
        dy = y[iframe] - grid_y[iy]
        dy = dy - round(dy/boxsize[2])*boxsize[2];
        if abs(dy) < ry
          dy = dy/bandwidth[2]
          f_y[iy, tid] += exp(-0.5*dy*dy)/(sqrt(2.0*pi)*bandwidth[2])
          iy_count += 1
          iy_array[iy_count, tid] = iy
        end
    end

    for i in 1:ix_count
        for j in 1:iy_count
            ix = ix_array[i, tid]
            iy = iy_array[j, tid]
            f_private[ix, iy, tid] += w*f_x[ix, tid]*f_y[iy, tid]
        end
    end
end

function convolution3d!(tid, iframe,
    x, nx, rx, gridsize_x, grid_x,
    y, ny, ry, gridsize_y, grid_y,
    z, nz, rz, gridsize_z, grid_z,
    bandwidth, f_x, f_y, f_z, f_private, w)

    f_x[:, tid] .= 0.0
    f_y[:, tid] .= 0.0
    f_z[:, tid] .= 0.0

    dx = x[iframe] - grid_x[1]
    ix_min = floor(Int, (dx - rx)/gridsize_x) + 1
    ix_min = ix_min >= 1 ? ix_min : 1
    ix_max = floor(Int, (dx + rx)/gridsize_x) + 2
    ix_max = ix_max <= nx ? ix_max : nx
    for ix in ix_min:ix_max
        dx = (grid_x[ix] - x[iframe])/bandwidth[1]
        f_x[ix, tid] += exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth[1])
    end

    dy = y[iframe] - grid_y[1]
    iy_min = floor(Int, (dy - ry)/gridsize_y) + 1
    iy_min = iy_min >= 1 ? iy_min : 1
    iy_max = floor(Int, (dy + ry)/gridsize_y) + 2
    iy_max = iy_max <= ny ? iy_max : ny
    for iy in iy_min:iy_max
        dy = (grid_y[iy] - y[iframe])/bandwidth[2]
        f_y[iy, tid] += exp(-0.5*dy*dy)/(sqrt(2.0*pi)*bandwidth[2])
    end

    dz = z[iframe] - grid_z[1]
    iz_min = floor(Int, (dz - rz)/gridsize_z) + 1
    iz_min = iz_min >= 1 ? iz_min : 1
    iz_max = floor(Int, (dz + rz)/gridsize_z) + 2
    iz_max = iz_max <= nz ? iz_max : nz
    for iz in iz_min:iz_max
        dz = (grid_z[iz] - z[iframe])/bandwidth[3]
        f_z[iz, tid] += exp(-0.5*dz*dz)/(sqrt(2.0*pi)*bandwidth[3])
    end

    for ix in ix_min:ix_max
        for iy in iy_min:iy_max
            for iz in iz_min:iz_max
                f_private[ix, iy, iz, tid] += w*f_x[ix, tid]*f_y[iy, tid]*f_z[iz, tid]
            end
        end
    end
end

function convolution3d_box!(tid, iframe, boxsize,
    x, nx, rx, gridsize_x, grid_x,
    y, ny, ry, gridsize_y, grid_y,
    z, nz, rz, gridsize_z, grid_z,
    bandwidth, f_x, f_y, f_z, f_private, w,
    ix_array, iy_array, iz_array)

    f_x[:, tid] .= 0.0
    f_y[:, tid] .= 0.0
    f_z[:, tid] .= 0.0

    ix_count = 0
    for ix in 1:nx
        dx = x[iframe] - grid_x[ix]
        dx = dx - round(dx/boxsize[1])*boxsize[1];
        if abs(dx) < rx
          dx = dx/bandwidth[1]
          f_x[ix, tid] += exp(-0.5*dx*dx)/(sqrt(2.0*pi)*bandwidth[1])
          ix_count += 1
          ix_array[ix_count, tid] = ix
        end
    end

    iy_count = 0
    for iy in 1:ny
        dy = y[iframe] - grid_y[iy]
        dy = dy - round(dy/boxsize[2])*boxsize[2];
        if abs(dy) < ry
          dy = dy/bandwidth[2]
          f_y[iy, tid] += exp(-0.5*dy*dy)/(sqrt(2.0*pi)*bandwidth[2])
          iy_count += 1
          iy_array[iy_count, tid] = iy
        end
    end

    iz_count = 0
    for iz in 1:nz
        dz = z[iframe] - grid_z[iz]
        dz = dz - round(dz/boxsize[3])*boxsize[3];
        if abs(dz) < rz
          dz = dz/bandwidth[3]
          f_z[iz, tid] += exp(-0.5*dz*dz)/(sqrt(2.0*pi)*bandwidth[3])
          iz_count += 1
          iz_array[iz_count, tid] = iz
        end
    end

    for i in 1:ix_count
        for j in 1:iy_count
            for k in 1:iz_count
                ix = ix_array[i, tid]
                iy = iy_array[j, tid]
                iz = iz_array[k, tid]
                f_private[ix, iy, iz, tid] += w*f_x[ix, tid]*f_y[iy, tid]*f_z[iz, tid]
            end
        end
    end
end

"""
ksdensity

kernel density estimator with gaussian kernel
"""
function ksdensity(x::Vector{Float64};
    grid_x,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Union{Float64,Missing}=missing,
    bandwidth::Union{Float64,Missing}=missing)
    nframe = length(x)
    if isempty(grid_x)
        grid_x2 = range(minimum(x), stop=maximum(x), length=100)
    else typeof(grid_x) != Vector{Float64}
        grid_x2 = map(Float64, grid_x)
    end
    nx = length(grid_x2)
    gridsize_x = abs(grid_x2[2]-grid_x2[1])
    if isempty(weight)
        weight2 = ones(Float64, nframe)/nframe
    else
        weight2 = weight
    end
    if typeof(bandwidth) == Missing
        bandwidth = gridsize_x
    end
    rx = bandwidth*5.0
    nth = Threads.nthreads()
    f_x = zeros(nx, nth)
    d = div(nframe, nth)
    if typeof(boxsize) == Missing
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution1d!(tid, iframe, x, nx, rx, gridsize_x, grid_x2, bandwidth, f_x, weight2[iframe])
            end
        end
    else
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution1d_box!(tid, iframe, boxsize, x, nx, rx, gridsize_x, grid_x2, bandwidth, f_x, weight2[iframe])
            end
        end
    end
    f = zeros(nx)
    for ind in 1:nx
        f[ind] = sum(f_x[ind, :])
    end
    f, grid_x2
end

function ksdensity_serial(x::Vector{Float64};
    grid_x,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Union{Float64,Missing}=missing,
    bandwidth::Union{Float64,Missing}=missing)
    nframe = length(x)
    if isempty(grid_x)
        grid_x2 = range(minimum(x), stop=maximum(x), length=100)
    else typeof(grid_x) != Vector{Float64}
        grid_x2 = map(Float64, grid_x)
    end
    nx = length(grid_x2)
    gridsize_x = abs(grid_x2[2]-grid_x2[1])
    if isempty(weight)
        weight2 = ones(Float64, nframe)/nframe
    else
        weight2 = weight
    end
    if typeof(bandwidth) == Missing
        bandwidth = gridsize_x
    end
    rx = bandwidth*5.0
    nth = 1
    f_x = zeros(nx, nth)
    if typeof(boxsize) == Missing
        for iframe in 1:nframe
            convolution1d!(1, iframe, x, nx, rx, gridsize_x, grid_x2, bandwidth, f_x, weight2[iframe])
        end
    else
        for iframe in 1:nframe
            convolution1d_box!(1, iframe, boxsize, x, nx, rx, gridsize_x, grid_x2, bandwidth, f_x, weight2[iframe])
        end
    end
    f = zeros(nx)
    for ind in 1:nx
        f[ind] = sum(f_x[ind, :])
    end
    f, grid_x2
end

function ksdensity(x::Vector{Float64}, y::Vector{Float64};
    grid_x, grid_y,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Vector{Float64}=Vector{Float64}(undef, 0),
    bandwidth::Vector{Float64}=Vector{Float64}(undef, 0))
    nframe = length(x)
    if isempty(grid_x)
        grid_x2 = range(minimum(x), stop=maximum(x), length=100)
    else typeof(grid_x) != Vector{Float64}
        grid_x2 = map(Float64, grid_x)
    end
    if isempty(grid_y)
        grid_y2 = range(minimum(y), stop=maximum(y), length=100)
    else typeof(grid_y) != Vector{Float64}
        grid_y2 = map(Float64, grid_y)
    end
    nx = length(grid_x2)
    ny = length(grid_y2)
    gridsize_x = abs(grid_x2[2]-grid_x2[1])
    gridsize_y = abs(grid_y2[2]-grid_y2[1])
    if isempty(weight)
        weight2 = ones(Float64, nframe)/nframe
    else
        weight2 = weight
    end
    if isempty(bandwidth)
        bandwidth = [gridsize_x, gridsize_y]
    end
    rx = bandwidth[1]*5.0
    ry = bandwidth[2]*5.0
    nth = Threads.nthreads()
    f_x = zeros(nx, nth)
    f_y = zeros(ny, nth)
    f_private = zeros(nx, ny, nth)
    d = div(nframe, nth)
    if isempty(boxsize)
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution2d!(tid, iframe,
                    x, nx, rx, gridsize_x, grid_x2,
                    y, ny, ry, gridsize_y, grid_y2,
                    bandwidth, f_x, f_y, f_private, weight2[iframe])
            end
        end
    else
        ix_array = zeros(Int, nx, nth)
        iy_array = zeros(Int, ny, nth)
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution2d_box!(tid, iframe, boxsize,
                    x, nx, rx, gridsize_x, grid_x2,
                    y, ny, ry, gridsize_y, grid_y2,
                    bandwidth, f_x, f_y, f_private, weight2[iframe],
                    ix_array, iy_array)
            end
        end

    end
    f = zeros(nx, ny)
    for ix in 1:nx
        for iy in 1:ny
            f[ix, iy] = sum(f_private[ix, iy, :])
        end
    end
    transpose(f), grid_x2, grid_y2
end

function ksdensity(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64};
    grid_x, grid_y, grid_z,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Vector{Float64}=Vector{Float64}(undef, 0),
    bandwidth::Vector{Float64}=Vector{Float64}(undef, 0))
    nframe = length(x)
    if isempty(grid_x)
        grid_x2 = range(minimum(x), stop=maximum(x), length=100)
    else typeof(grid_x) != Vector{Float64}
        grid_x2 = map(Float64, grid_x)
    end
    if isempty(grid_y)
        grid_y2 = range(minimum(y), stop=maximum(y), length=100)
    else typeof(grid_y) != Vector{Float64}
        grid_y2 = map(Float64, grid_y)
    end
    if isempty(grid_z)
        grid_z2 = range(minimum(z), stop=maximum(z), length=100)
    else typeof(grid_z) != Vector{Float64}
        grid_z2 = map(Float64, grid_z)
    end
    nx = length(grid_x2)
    ny = length(grid_y2)
    nz = length(grid_z2)
    gridsize_x = abs(grid_x2[2]-grid_x2[1])
    gridsize_y = abs(grid_y2[2]-grid_y2[1])
    gridsize_z = abs(grid_z2[2]-grid_z2[1])
    if isempty(weight)
        weight2 = ones(Float64, nframe)/nframe
    else
        weight2 = weight
    end
    if isempty(bandwidth)
        bandwidth = [gridsize_x, gridsize_y, gridsize_z]
    end
    rx = bandwidth[1]*5.0
    ry = bandwidth[2]*5.0
    rz = bandwidth[3]*5.0
    nth = Threads.nthreads()
    f_x = zeros(nx, nth)
    f_y = zeros(ny, nth)
    f_z = zeros(nz, nth)
    f_private = zeros(nx, ny, nz, nth)
    d = div(nframe, nth)
    if isempty(boxsize)
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution3d!(tid, iframe,
                    x, nx, rx, gridsize_x, grid_x2,
                    y, ny, ry, gridsize_y, grid_y2,
                    z, nz, rz, gridsize_z, grid_z2,
                    bandwidth, f_x, f_y, f_z, f_private, weight2[iframe])
            end
        end
    else
        ix_array = zeros(Int, nx, nth)
        iy_array = zeros(Int, ny, nth)
        iz_array = zeros(Int, nz, nth)
        Threads.@threads for id in 1:nth
            tid = Threads.threadid()
            frames = (tid-1)*d+1:tid*d
            for iframe in frames
                convolution3d_box!(tid, iframe, boxsize,
                    x, nx, rx, gridsize_x, grid_x2,
                    y, ny, ry, gridsize_y, grid_y2,
                    z, nz, rz, gridsize_z, grid_z2,
                    bandwidth, f_x, f_y, f_z, f_private, weight2[iframe],
                    ix_array, iy_array, iz_array)
            end
        end
    end
    f = zeros(nx, ny, nz)
    for ix in 1:nx
        for iy in 1:ny
            for iz in 1:nz
                f[ix, iy, iz] = sum(f_private[ix, iy, iz, :])
            end
        end
    end
    f, grid_x2, grid_y2, grid_z2
end

"""
compute_pmf

Potential of mean force estimator by a kernel density estimator
"""
function compute_pmf(x::Vector{Float64};
    grid_x,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Union{Float64,Missing}=missing,
    bandwidth::Union{Float64,Missing}=missing)

    pdf, grid_x = ksdensity(x, grid_x=grid_x, weight=weight, boxsize=boxsize, bandwidth=bandwidth)
    #pdf[pdf .< eps()] .= missing
    pmf = -log.(pdf)
    pmf = pmf .- minimum(pmf)
    pmf, grid_x
end

function compute_pmf(x::Vector{Float64}, y::Vector{Float64};
    grid_x, grid_y,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Vector{Float64}=Vector{Float64}(undef, 0),
    bandwidth::Vector{Float64}=Vector{Float64}(undef, 0))

    pdf, grid_x, grid_y = ksdensity(x, y, grid_x=grid_x, grid_y=grid_y, weight=weight, boxsize=boxsize, bandwidth=bandwidth)
    #pdf[pdf .< eps()] .= NaN
    pmf = -log.(pdf)
    pmf = pmf .- minimum(pmf)
    pmf, grid_x, grid_y
end

function compute_pmf(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64};
    grid_x, grid_y, grid_z,
    weight::Vector{Float64}=Vector{Float64}(undef, 0),
    boxsize::Vector{Float64}=Vector{Float64}(undef, 0),
    bandwidth::Vector{Float64}=Vector{Float64}(undef, 0))

    pdf, grid_x, grid_y, grid_z = ksdensity(x, y, z, grid_x=grid_x, grid_y=grid_y, grid_z=grid_z, weight=weight, boxsize=boxsize, bandwidth=bandwidth)
    #pdf[pdf .< eps()] .= NaN
    pmf = -log.(pdf)
    pmf = pmf .- minimum(pmf)
    pmf, grid_x, grid_y, grid_z
end


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking
using StatsPlots

ProjDir = @__DIR__

for i in 5:7
  include("$(ProjDir)/m5.$i.jl")
end

# ### snippet 5.29

println()
df = DataFrame!(CSV.File(rel_path("..", "data", "milk.csv"), delim=';'));
df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
df[!, :neocortex_perc] = parse.(Float64, df[:, :neocortex_perc])
df[!, :lmass] = log.(df[:, :mass])
#first(df, 5) |> display

# ### snippet 5.1

scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])
println()

m_5_7 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] K; // Outcome
 vector[N] NC; // Predictor
 vector[N] M; // Predictor
}

parameters {
 real a; // Intercept
 real bM; // Slope (regression coefficients)
 real bN; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);           //Priors
  bN ~ normal(0, 0.5);
  bM ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bM * M + bN * NC;
  K ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5_7s = SampleModel("m5.7", m_5_7);

# Input data for cmdstan

m5_7_data = Dict("N" => size(df, 1), "M" => df[!, :lmass_s],
    "K" => df[!, :kcal_per_g_s], "NC" => df[!, :neocortex_perc_s]);

# Sample using StanSample

rc = stan_sample(m5_7s, data=m5_7_data);

if success(rc)

  # Describe the draws

  dfa7 = read_samples(m5_7s; output_format=:dataframe)
  println("Normal estimate:")
  p = Particles(dfa7) |> display
  println("Quap estimates:")
  quap(dfa7) |> display

  r1 = plotcoef([m5_5s, m5_6s, m5_7s], [:a, :bN, :bM], "$(ProjDir)/Fig-38a.png",
    "bN & bM Normal estimates")
  r1 |> display

  r2 = plotcoef([m5_5s, m5_6s, m5_7s], [:a, :bN, :bM], "$(ProjDir)/Fig-38b.png",
    "bN & bM Quap estimates", quap)
  r2 |> display
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module OEISUtils
using HTTP, Nemo

export ModuleOEISUtils
export oeis_writebfile, oeis_trimdata, oeis_remote, oeis_local, is_oeis_installed
export oeis_notinstalled, oeis_path, oeis_search, oeis_readbfile

"""

A collection of utilities for handling OEIS related tasks.
"""
const ModuleOEISUtils = ""

# Directory of OEIS data.
srcdir = realpath(joinpath(dirname(@__FILE__)))
ROOTDIR = dirname(srcdir)
datadir = joinpath(ROOTDIR, "data")

"""

Returns the path where the oeis data is expected.
"""
oeis_path() = joinpath(datadir, "stripped")

"""

Indicates if the local copy of the OEIS data (the so-called 'stripped' file) is installed (in ../data).
"""
is_oeis_installed() = isfile(oeis_path())

"""

Indicates if the local copy of the OEIS data (the so-called 'stripped' file) is not installed and warns.
"""
function oeis_notinstalled()
    if !is_oeis_installed()
        @warn("OEIS data not installed! Download stripped.gz from oeis.org,")
        @warn("expand it and put it in the directory ../data.")
        return true
    end
    return false
end

"""

Write a so-called b-file for submission to the OEIS. The file is saved in the 'data' directory.
"""
function oeis_writebfile(anum::String, fun::Function, range::OrdinalRange)

    if !occursin(r"^A[0-9]{6}$", anum)
        @warn("Not a valid A-number!")
        return
    end

    filename = joinpath(datadir, "b" * anum[2:end] * ".txt")
    @info("Writing " * anum * " to " * filename)

    f = open(filename, "w")
    for n ∈ range
        println(f, n, " ", fun(n))
    end
    close(f)
end

function oeis_writebfile(anum, list, offset = 1)

    if !occursin(r"^A[0-9]{6}$", anum)
        @warn("Not a valid A-number!")
        return
    end

    filename = joinpath(datadir, "b" * anum[2:end] * ".txt")
    @info("Writing " * anum * " to " * filename)

    n = offset
    f = open(filename, "w")
    for l ∈ list
        println(f, n, " ", list[n])
        n += 1
    end
    close(f)
end

"""

Read a file in the so-called b-file format of the OEIS.
"""
function oeis_readbfile(filename)

    @info("Reading " * filename)

    file = open(filename, "r")
    lines = readlines(file)

    for l ∈ lines
        s = split(l, ' ')
        i = parse(Int, s[1])
        v = parse(Int, s[2])
        println(i, " ", v)
    end
    close(file)
end

"""

Make sure that the length of the data section of an OEIS entry does not exceed 260 characters.
"""
function oeis_trimdata(fun, offset::Int)
    len = n = 0
    S = ""
    while true
        st = string(fun(offset + n))
        len += length(st)
        len > 260 && break
        S *= st * ", "
        len += 2
        n += 1
    end
    println(n, " terms")
    return S
end

"""

Download the sequence with A-number 'anum' from the OEIS to a file in json format. The file is saved in the 'data' directory.
"""
function oeis_remote(anum)
    if !occursin(r"^A[0-9]{6}$", anum)
        @warn("Not a valid A-number!")
        return
    end

    filename = joinpath(datadir, anum * ".json")

    url = HTTP.URI("http://oeis.org/search?q=id:" * anum * "&fmt=json")
    tries = 3
    r = nothing
    for i = 1:tries
        try
            r = HTTP.get(url; readtimeout = 2)
            getfield(r, :status) == 200 && break
            getfield(r, :status) == 302 && break
        catch e
            @warn(e)
        end
        sleep(2)
    end
    if r ≠ nothing && getfield(r, :status) == 200
        open(filename, "w") do f
            write(f, getfield(r, :body))
        end
        @info("Dowloaded " * basename(filename) * " to " * datadir)
    else
        if r === nothing
            @warn("Could not download $url, connection timed out.\n")
        else
            @warn("Could not download $url\nStatus: $(getfield(r, :status))")
        end
    end
end

"""

Get the sequence with A-number 'anum' from a local copy of the expanded 'stripped.gz' file which can be downloaded from the OEIS. 'bound' is an upper bound for the number of terms returned. The 'stripped' file is assumed to be in the '../data' directory.
"""
function oeis_local(anum::String, bound::Int)

    if !occursin(r"^A[0-9]{6}$", anum)
        @warn(anum * " is not a valid A-number!")
        return []
    end

    oeis_notinstalled() && return []

    A = Array{String}
    data = open(oeis_path())
    for ln ∈ eachline(data)
        if startswith(ln, anum)
            A = split(chop(chomp(ln)), ","; limit = bound + 2)
            break
        end
    end
    close(data)

    [convert(fmpz, parse(BigInt, n)) for n ∈ A[2:min(bound + 1, end)]]
end

"""

Search for a sequence in the local OEIS database ('../data/stripped'). Input the sequence as a comma separated string. The search is redone 'restarts' times in the case that no match was found with the first remaining term removed from the search string. Prints the matches.
"""
function oeis_search(seq::String, restarts::Int, verbose=false)

    oeis_notinstalled() && return []

    found = false
    seq = replace(seq, ' ' => "")
    verbose && println("Searching for:")
    verbose && println(seq)

    data = open(oeis_path())
    foundMAX = 4  # bound the search (sequences with many '0' make trouble).
    for line ∈ eachline(data)
        index = findfirst(seq, line)
        index === nothing && continue
        println("Searched for $seq")
        println("$line in range ", index)
        found = true
        foundMAX -= 1
        if foundMAX <= 0
            @warn("Truncating the search! There might be too many others sequences!")
            break
        end
    end
    close(data)

    if !found && restarts > 0
        ind = findfirst(isequal(','), seq)
        if !(ind === nothing) && (length(seq) > ind)
            seq = seq[ind+1:end]
            verbose && println("Restarting with fewer terms.")
            oeis_search(seq, restarts-1, verbose)
        end
    end
end

"""

Search for a sequence in the local OEIS database starting from the second term.
"""
function oeis_search(seq::Array{fmpz,1}, restarts::Int, verbose=false)
    oeis_search("$seq"[2:end-1], restarts, verbose)
end

function oeis_search(seq::Array{Int,1}, restarts::Int, verbose=false)
    oeis_search("$seq"[2:end-1], restarts, verbose)
end

function oeis_search(seq::Array{BigInt,1}, restarts::Int, verbose=false)
    oeis_search("$seq"[2:end-1], restarts, verbose)
end

#START-TEST-########################################################

function test()
    if is_oeis_installed()
        @info("OEIS data is installed as: " * oeis_path())
    end
end

function demo()
    oeis_writebfile("A000290", n -> n * n, 0:10)
    bfilepath = joinpath(datadir, "b000290.txt")
    oeis_readbfile(bfilepath)

    println(oeis_trimdata(n -> (-1)^n * n^5, 10))

    oeis_remote("A123456")

    println(oeis_local("A123456", 12))
    println(oeis_local("A015108", 11))

    seq = "0, 1, 2, 3, 6, 9, 14, 22, 32, 46, 66, 93, 128, 176, 238, 319"
    oeis_search(seq, 2)
    oeis_search([0, 1, 2, 3, 6, 9, 14, 22, 32, 46, 66, 93, 128, 176], 2)
end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
import Base:
  # Types
  Array,
  CartesianIndex,
  IndexStyle,
  Tuple,
  # Symbols
  +,
  -,
  *,
  # Methods
  checkbounds,
  complex,
  convert,
  conj,
  copy,
  copyto!,
  eachindex,
  eltype,
  fill,
  fill!,
  getindex,
  hash,
  isempty,
  isless,
  iterate,
  length,
  ndims,
  permutedims,
  permutedims!,
  promote_rule,
  randn,
  reshape,
  setindex,
  setindex!,
  show,
  size,
  stride,
  strides,
  summary,
  to_indices,
  unsafe_convert,
  view,
  zeros

import Base.Broadcast: Broadcasted, BroadcastStyle

import LinearAlgebra: exp

import TupleTools: isperm


text:
using MonteCarlo, MonteCarloObservable
using Test
using Random

@testset "All Tests" begin
    @testset "Lattices" begin
        include("lattices.jl")
    end

    @testset "Models" begin
        include("modeltests_IsingModel.jl")
        include("modeltests_HubbardModelAttractive.jl")
    end

    @testset "Flavors" begin
        # include("flavortests_MC.jl")
        include("flavortests_DQMC.jl")
    end

    @testset "Measurements" begin
        include("measurements.jl")
    end

    @testset "Intergration tests" begin
        include("integration_tests.jl")
    end

    @testset "Exact Diagonalization" begin
        include("ED/ED_tests.jl")
    end
end


text:
###
#   Utilities and show
###

left(x :: Interval) = x.lo
right(x :: Interval) = x.hi

left(x :: IntervalUnion) = left.(x.v)
right(x :: IntervalUnion) = right.(x.v)

isequal(x :: IntervalUnion, y ::IntervalUnion) = all(isequal.(sort(x.v),sort(y.v)))

function in( x:: IntervalUnion, y :: Vector{IntervalUnion})
    for i =1:length(y)
        if x == y[i]
            return true
        end
    end
    return false
end

function Base.show(io::IO, this::IntervalUnion)
    if length(this.v) == 1;
        print(io, "$(this.v[1])ᵤ");
    else
        v = sort(this.v)
        print(io, join(v, " ∪ "));
    end
end


text:
module EAUtils

export code_escapes, @code_escapes, __clear_cache!

import ..EscapeAnalysis
const EA = EscapeAnalysis
const CC = Core.Compiler

# entries
# -------

import InteractiveUtils: gen_call_with_extracted_types_and_kwargs

"""
    @code_escapes [options...] f(args...)

Evaluates the arguments to the function call, determines its types, and then calls
[`code_escapes`](@ref) on the resulting expression.
As with `@code_typed` and its family, any of `code_escapes` keyword arguments can be given
as the optional arguments like `@code_escapes debuginfo=:source myfunc(myargs...)`.
"""
macro code_escapes(ex0...)
    return gen_call_with_extracted_types_and_kwargs(__module__, :code_escapes, ex0)
end

"""
    code_escapes(f, argtypes=Tuple{}; [world], [interp], [debuginfo]) -> result::EscapeResult
    code_escapes(tt::Type{<:Tuple}; [world], [interp], [debuginfo]) -> result::EscapeResult

Runs the escape analysis on optimized IR of a generic function call with the given type signature.
Note that the escape analysis runs after inlining, but before any other optimizations.

```julia
julia> mutable struct SafeRef{T}
           x::T
       end

julia> Base.getindex(x::SafeRef) = x.x;

julia> Base.isassigned(x::SafeRef) = true;

julia> get′(x) = isassigned(x) ? x[] : throw(x);

julia> result = code_escapes((String,String,String,String)) do s1, s2, s3, s4
           r1 = Ref(s1)
           r2 = Ref(s2)
           r3 = SafeRef(s3)
           try
               s1 = get′(r1)
               ret = sizeof(s1)
           catch err
               global g = err # will definitely escape `r1`
           end
           s2 = get′(r2)      # still `r2` doesn't escape fully
           s3 = get′(r3)      # still `r3` doesn't escape fully
           s4 = sizeof(s4)    # the argument `s4` doesn't escape here
           return s2, s3, s4
       end
#1(X _2::String, ↑ _3::String, ↑ _4::String, ✓ _5::String) in Main at REPL[6]:2
X  1 ── %1  = %new(Base.RefValue{String}, _2)::Base.RefValue{String}
*′ │    %2  = %new(Base.RefValue{String}, _3)::Base.RefValue{String}
✓′ └─── %3  = %new(SafeRef{String}, _4)::SafeRef{String}
◌  2 ── %4  = \$(Expr(:enter, #8))
✓′ │    %5  = ϒ (%3)::SafeRef{String}
*′ │    %6  = ϒ (%2)::Base.RefValue{String}
✓  └─── %7  = ϒ (_5)::String
◌  3 ── %8  = Base.isdefined(%1, :x)::Bool
◌  └───       goto #5 if not %8
X  4 ──       Base.getfield(%1, :x)::String
◌  └───       goto #6
◌  5 ──       Main.throw(%1)::Union{}
◌  └───       unreachable
◌  6 ──       nothing::typeof(Core.sizeof)
◌  │          nothing::Int64
◌  └───       \$(Expr(:leave, 1))
◌  7 ──       goto #10
✓′ 8 ── %18 = φᶜ (%5)::SafeRef{String}
*′ │    %19 = φᶜ (%6)::Base.RefValue{String}
✓  │    %20 = φᶜ (%7)::String
◌  └───       \$(Expr(:leave, 1))
X  9 ── %22 = \$(Expr(:the_exception))::Any
◌  │          (Main.g = %22)::Any
◌  └───       \$(Expr(:pop_exception, :(%4)))::Any
✓′ 10 ┄ %25 = φ (#7 => %3, #9 => %18)::SafeRef{String}
*′ │    %26 = φ (#7 => %2, #9 => %19)::Base.RefValue{String}
✓  │    %27 = φ (#7 => _5, #9 => %20)::String
◌  │    %28 = Base.isdefined(%26, :x)::Bool
◌  └───       goto #12 if not %28
↑  11 ─ %30 = Base.getfield(%26, :x)::String
◌  └───       goto #13
◌  12 ─       Main.throw(%26)::Union{}
◌  └───       unreachable
↑  13 ─ %34 = Base.getfield(%25, :x)::String
◌  │    %35 = Core.sizeof::typeof(Core.sizeof)
◌  │    %36 = (%35)(%27)::Int64
↑  │    %37 = Core.tuple(%30, %34, %36)::Tuple{String, String, Int64}
◌  └───       return %37
```

The symbols in the side of each call argument and SSA statements represents the following meaning:
- `◌` (plain): this value is not analyzed because escape information of it won't be used anyway (when the object is `isbitstype` for example)
- `✓` (green or cyan): this value never escapes (`has_no_escape(result.state[x])` holds), colored blue if it has arg escape also (`has_arg_escape(result.state[x])` holds)
- `↑` (blue or yellow): this value can escape to the caller via return (`has_return_escape(result.state[x])` holds), colored yellow if it has unhandled thrown escape also (`has_thrown_escape(result.state[x])` holds)
- `X` (red): this value can escape to somewhere the escape analysis can't reason about like escapes to a global memory (`has_all_escape(result.state[x])` holds)
- `*` (bold): this value's escape state is between the `ReturnEscape` and `AllEscape` in the partial order of [`EscapeInfo`](@ref), colored yellow if it has unhandled thrown escape also (`has_thrown_escape(result.state[x])` holds)
- `′`: this value has additional object field / array element information in its `AliasInfo` property

For testing, escape information of each call argument and SSA value can be inspected programmatically as like:
```julia
julia> result.state[Core.Argument(3)]
ReturnEscape

julia> result.state[Core.SSAValue(3)]
NoEscape′
```
"""
function code_escapes(@nospecialize(args...);
                      world::UInt = get_world_counter(),
                      interp::Core.Compiler.AbstractInterpreter = Core.Compiler.NativeInterpreter(world),
                      debuginfo::Symbol = :none,
                      optimize::Bool = true)
    interp = EscapeAnalyzer(interp, optimize)
    empty!(LOCAL_ESCAPE_CACHE)
    results = code_typed(args...; optimize=true, world, interp)
    isone(length(results)) || throw(ArgumentError("`code_escapes` only supports single analysis result"))
    return EscapeResult(interp.ir, interp.state, interp.linfo, debuginfo===:source)
end

__clear_cache!() = (empty!(GLOBAL_ESCAPE_CACHE); empty!(GLOBAL_CODE_CACHE))

# AbstractInterpreter
# -------------------

# imports
import .CC:
    AbstractInterpreter,
    NativeInterpreter,
    WorldView,
    WorldRange,
    InferenceParams,
    OptimizationParams,
    get_world_counter,
    get_inference_cache,
    lock_mi_inference,
    unlock_mi_inference,
    add_remark!,
    may_optimize,
    may_compress,
    may_discard_trees,
    verbose_stmt_info,
    code_cache,
    @timeit,
    get_inference_cache,
    convert_to_ircode,
    slot2reg,
    compact!,
    ssa_inlining_pass!,
    sroa_pass!,
    adce_pass!,
    type_lift_pass!,
    JLOptions,
    verify_ir,
    verify_linetable
# usings
import Core:
    CodeInstance, MethodInstance, CodeInfo
import .CC:
    InferenceResult, OptimizationState, IRCode, copy as cccopy
import .EA:
    analyze_escapes, cache_escapes!, LOCAL_ESCAPE_CACHE, GLOBAL_ESCAPE_CACHE

mutable struct EscapeAnalyzer{State} <: AbstractInterpreter
    native::NativeInterpreter
    optimize::Bool
    ir::IRCode
    state::State
    linfo::MethodInstance
    EscapeAnalyzer(native::NativeInterpreter, optimize::Bool) =
        new{EscapeState}(native, optimize)
end

CC.InferenceParams(interp::EscapeAnalyzer)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::EscapeAnalyzer) = OptimizationParams(interp.native)
CC.get_world_counter(interp::EscapeAnalyzer)  = get_world_counter(interp.native)

CC.lock_mi_inference(::EscapeAnalyzer,   ::MethodInstance) = nothing
CC.unlock_mi_inference(::EscapeAnalyzer, ::MethodInstance) = nothing

CC.add_remark!(interp::EscapeAnalyzer, sv, s) = add_remark!(interp.native, sv, s)

CC.may_optimize(interp::EscapeAnalyzer)      = may_optimize(interp.native)
CC.may_compress(interp::EscapeAnalyzer)      = may_compress(interp.native)
CC.may_discard_trees(interp::EscapeAnalyzer) = may_discard_trees(interp.native)
CC.verbose_stmt_info(interp::EscapeAnalyzer) = verbose_stmt_info(interp.native)

CC.get_inference_cache(interp::EscapeAnalyzer) = get_inference_cache(interp.native)

const GLOBAL_CODE_CACHE = IdDict{MethodInstance,CodeInstance}()

function CC.code_cache(interp::EscapeAnalyzer)
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(GlobalCache(), worlds)
end

struct GlobalCache end

CC.haskey(wvc::WorldView{GlobalCache}, mi::MethodInstance) = haskey(GLOBAL_CODE_CACHE, mi)

CC.get(wvc::WorldView{GlobalCache}, mi::MethodInstance, default) = get(GLOBAL_CODE_CACHE, mi, default)

CC.getindex(wvc::WorldView{GlobalCache}, mi::MethodInstance) = getindex(GLOBAL_CODE_CACHE, mi)

function CC.setindex!(wvc::WorldView{GlobalCache}, ci::CodeInstance, mi::MethodInstance)
    GLOBAL_CODE_CACHE[mi] = ci
    add_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_cache!)
        end
    end
    return nothing
end

function invalidate_cache!(replaced, max_world, depth = 0)
    delete!(GLOBAL_CODE_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !haskey(GLOBAL_CODE_CACHE, mi)
                continue # otherwise fall into infinite loop
            end
            invalidate_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

function CC.optimize(interp::EscapeAnalyzer,
    opt::OptimizationState, params::OptimizationParams, caller::InferenceResult)
    ir = run_passes_with_ea(interp, opt.src, opt, caller)
    return CC.finish(interp, opt, params, ir, caller)
end

function CC.cache_result!(interp::EscapeAnalyzer, caller::InferenceResult)
    GLOBAL_ESCAPE_CACHE[caller.linfo] = LOCAL_ESCAPE_CACHE[caller]
    return Base.@invoke CC.cache_result!(interp::AbstractInterpreter, caller::InferenceResult)
end

function run_passes_with_ea(interp::EscapeAnalyzer, ci::CodeInfo, sv::OptimizationState,
    caller::InferenceResult)
    @timeit "convert"   ir = convert_to_ircode(ci, sv)
    @timeit "slot2reg"  ir = slot2reg(ir, ci, sv)
    # TODO: Domsorting can produce an updated domtree - no need to recompute here
    @timeit "compact 1" ir = compact!(ir)
    nargs = let def = sv.linfo.def; isa(def, Method) ? Int(def.nargs) : 0; end
    local state
    try
        @timeit "[IPO EA]" state, callinfo = analyze_escapes(ir, nargs, true)
        Core.eval(Main, :(callinfo = $callinfo; ir = $(cccopy(ir))))
    catch err
        @error "error happened within [IPO EA], insepct `Main.ir` and `Main.nargs`"
        @eval Main (ir = $ir; nargs = $nargs)
        rethrow(err)
    end
    # cache this result
    cache_escapes!(caller, state, cccopy(ir))
    if !interp.optimize
        # return back the result
        interp.ir = cccopy(ir)
        interp.state = state
        interp.linfo = sv.linfo
    end
    @timeit "Inlining"  ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
    # @timeit "verify 2" verify_ir(ir)
    @timeit "compact 2" ir = compact!(ir)
    try
        @timeit "[Local EA]" state, = analyze_escapes(ir, nargs, false)
    catch err
        @error "error happened within [Local EA], insepct `Main.ir` and `Main.nargs`"
        @eval Main (ir = $ir; nargs = $nargs)
        rethrow(err)
    end
    if interp.optimize
        # return back the result
        interp.ir = cccopy(ir)
        interp.state = state
        interp.linfo = sv.linfo
    end
    @timeit "SROA"      ir = sroa_pass!(ir)
    @timeit "ADCE"      ir = adce_pass!(ir)
    @timeit "type lift" ir = type_lift_pass!(ir)
    @timeit "compact 3" ir = compact!(ir)
    if JLOptions().debug_level == 2
        @timeit "verify 3" (verify_ir(ir); verify_linetable(ir.linetable))
    end
    return ir
end

# printing
# --------

import Core: Argument, SSAValue
import .CC: widenconst, singleton_type
import .EA: EscapeInfo, EscapeState, EscapeCache

# in order to run a whole analysis from ground zero (e.g. for benchmarking, etc.)
__clear_caches!() = (__clear_code_cache!(); EA.__clear_escape_cache!())

function get_name_color(x::EscapeInfo, symbol::Bool = false)
    getname(x) = string(nameof(x))
    if x === EA.⊥
        name, color = (getname(EA.NotAnalyzed), "◌"), :plain
    elseif EA.has_no_escape(EA.ignore_argescape(x))
        if EA.has_arg_escape(x)
            name, color = (getname(EA.ArgEscape), "✓"), :cyan
        else
            name, color = (getname(EA.NoEscape), "✓"), :green
        end
    elseif EA.has_all_escape(x)
        name, color = (getname(EA.AllEscape), "X"), :red
    elseif EA.has_return_escape(x)
        name = (getname(EA.ReturnEscape), "↑")
        color = EA.has_thrown_escape(x) ? :yellow : :blue
    else
        name = (nothing, "*")
        color = EA.has_thrown_escape(x) ? :yellow : :bold
    end
    name = symbol ? last(name) : first(name)
    if name !== nothing && !isa(x.AliasInfo, Bool)
        name = string(name, "′")
    end
    return name, color
end

# pcs = sprint(show, collect(x.EscapeSites); context=:limit=>true)
function Base.show(io::IO, x::EscapeInfo)
    name, color = get_name_color(x)
    if isnothing(name)
        Base.@invoke show(io::IO, x::Any)
    else
        printstyled(io, name; color)
    end
end
function Base.show(io::IO, ::MIME"application/prs.juno.inline", x::EscapeInfo)
    name, color = get_name_color(x)
    if isnothing(name)
        return x # use fancy tree-view
    else
        printstyled(io, name; color)
    end
end

struct EscapeResult
    ir::IRCode
    state::EscapeState
    linfo::Union{Nothing,MethodInstance}
    source::Bool
    function EscapeResult(ir::IRCode, state::EscapeState,
        linfo::Union{Nothing,MethodInstance} = nothing,
        source::Bool=false)
        return new(ir, state, linfo, source)
    end
end
Base.show(io::IO, result::EscapeResult) = print_with_info(io, result)
@eval Base.iterate(res::EscapeResult, state=1) =
    return state > $(fieldcount(EscapeResult)) ? nothing : (getfield(res, state), state+1)

Base.show(io::IO, cached::EscapeCache) =
    show(io, EscapeResult(cached.ir, cached.state, nothing))

# adapted from https://github.com/JuliaDebug/LoweredCodeUtils.jl/blob/4612349432447e868cf9285f647108f43bd0a11c/src/codeedges.jl#L881-L897
function print_with_info(io::IO, (; ir, state, linfo, source)::EscapeResult)
    # print escape information on SSA values
    function preprint(io::IO)
        ft = ir.argtypes[1]
        f = singleton_type(ft)
        if f === nothing
            f = widenconst(ft)
        end
        print(io, f, '(')
        for i in 1:state.nargs
            arg = state[Argument(i)]
            i == 1 && continue
            c, color = get_name_color(arg, true)
            printstyled(io, c, ' ', '_', i, "::", ir.argtypes[i]; color)
            i ≠ state.nargs && print(io, ", ")
        end
        print(io, ')')
        if !isnothing(linfo)
            def = linfo.def
            printstyled(io, " in ", (isa(def, Module) ? (def,) : (def.module, " at ", def.file, ':', def.line))...; color=:bold)
        end
        println(io)
    end

    # print escape information on SSA values
    # nd = ndigits(length(ssavalues))
    function preprint(io::IO, idx::Int)
        c, color = get_name_color(state[SSAValue(idx)], true)
        # printstyled(io, lpad(idx, nd), ' ', c, ' '; color)
        printstyled(io, rpad(c, 2), ' '; color)
    end

    print_with_info(preprint, (args...)->nothing, io, ir, source)
end

function print_with_info(preprint, postprint, io::IO, ir::IRCode, source::Bool)
    io = IOContext(io, :displaysize=>displaysize(io))
    used = Base.IRShow.stmts_used(io, ir)
    if source
        line_info_preprinter = function (io::IO, indent::String, idx::Int)
            r = Base.IRShow.inline_linfo_printer(ir)(io, indent, idx)
            idx ≠ 0 && preprint(io, idx)
            return r
        end
    else
        line_info_preprinter = Base.IRShow.lineinfo_disabled
    end
    line_info_postprinter = Base.IRShow.default_expr_type_printer
    preprint(io)
    bb_idx_prev = bb_idx = 1
    for idx = 1:length(ir.stmts)
        preprint(io, idx)
        bb_idx = Base.IRShow.show_ir_stmt(io, ir, idx, line_info_preprinter, line_info_postprinter, used, ir.cfg, bb_idx)
        postprint(io, idx, bb_idx != bb_idx_prev)
        bb_idx_prev = bb_idx
    end
    max_bb_idx_size = ndigits(length(ir.cfg.blocks))
    line_info_preprinter(io, " "^(max_bb_idx_size + 2), 0)
    postprint(io)
    return nothing
end

end # module EAUtils


text:
using DataDeps

include("download_utils.jl")
include("loader.jl")

xsum_dependency = DataDep(
    "XSum-Preprocessed-BERT",
    """
    Preprocessed XSum summary data for BERT-based summarization.

    XSum dataset, preprocessed for summarization by Liu et al.
    Original and preprocessed dataset released by the University of Edinburgh.

    Citing:
    @inproceedings{narayan2018xsum,
        author={Shashi Narayan and Shay B. Cohen and Mirella Lapata},
        title={Don't Give Me the Details, Just the Summary! Topic-Aware Convolutional Neural Networks for Extreme Summarization},
        booktitle={Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing},
        year={2018},
        address={Brussels, Belgium},
    }
    @inproceedings{liu2019presumm,
        author={Yang Liu and Mirella Lapata},
        title={Text Summarization with Pretrained Encoders},
        booktitle={Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing}
        year={2019},
        address={Hong Kong, China},
    }
    """,
    "1BWBN1coTWGBqrWoOfRc5dhojPHhatbYs",
    "26bc6ff7dafe040fc3a14578e4dbb9193f276d34f1e7a310edb44b8bfd77384c",
    fetch_method=fetch_google_drive_id,
    post_fetch_method=unpack,
)

function xsum(type::String)::Channel{SummaryPair}
    return summaries(
        joinpath(
            datadep"XSum-Preprocessed-BERT",
            "bert_data_xsum_new"
        ),
        type
    )
end

cnn_dm_dependency = DataDep(
    "CNN-Daily-Mail-Preprocessed-BERT",
    """
    Preprocessed CNN / Daily Mail summary data for BERT-based summarization.

    CNN / Daily Mail dataset, preprocessed for summarization by Liu et al.
    Original dataset released by Google DeepMind and the University of Oxford.

    Citing:
    @inproceedings{hermann2015cnndm,
        author={Karl Moritz Hermann and Tom\\'a\\v{s} Ko\\v{c}isk\\'y and Edward Grefenstette and Lasse Espeholt and Will Kay and Mustafa Suleyman and Phil Blunsom},
        title={Teaching Machines to Read and Comprehend},
        booktitle={Advances in Neural Information Processing Systems},
        year={2015},
        address={Montreal, Quebec, Canada},
    }
    @inproceedings{liu2019presumm,
        author={Yang Liu and Mirella Lapata},
        title={Text Summarization with Pretrained Encoders},
        booktitle={Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing}
        year={2019},
        address={Hong Kong, China},
    }
    """,
    "1DN7ClZCCXsk2KegmC6t4ClBwtAf5galI",
    "e1921b399012b4ce7e585eabe8a73fea1e9541fc30892b21de2c503ab9a1f9a9",
    fetch_method=fetch_google_drive_id,
    post_fetch_method=unpack,
)

function cnn_dm(type::String)::Channel{SummaryPair}
    return summaries(
        joinpath(
            datadep"CNN-Daily-Mail-Preprocessed-BERT",
            "bert_data_cnndm_final"
        ),
        type
    )
end


text:
import ..GPU

import ..YOLO: loadWeights!

export loadWeights!

"""
    flipkernel(x)

Flips given kernel
"""
flipkernel(x) = x[end:-1:1, end:-1:1, :, :]

"""
    updateconv!(c, gamma, mean, variance)

Applies batch-normalization to given convolutional layer
"""
function updateconv!(c, gamma, mean, variance)
    gamma4 = reshape(gamma, 1, 1, 1, :)
    variance4 = reshape(variance, 1, 1, 1, :)
    gamma3 = reshape(gamma, 1, 1, :, 1)
    mean3 = reshape(mean, 1, 1, :, 1)
    variance3 = reshape(variance, 1, 1, :, 1)
    #update c.w
    a = convert(Float32, 0.001)
    c.w = (c.w .* gamma4) ./ sqrt.(variance4 .+ a)
    #update c.b
    c.b = c.b .- (gamma3 .* mean3 ./ sqrt.(variance3 .+ a))
    return c.w, c.b
end

"""
    loadWeights!(model::YOLO.v2_tiny.v2Chain, weights_filepath::String)

Loads layers' weights from given weights file path.
"""
function loadWeights!(model::Chain_layers, settings::Settings, nr_constants = 4)
    outlength = settings.cell_bboxes * (5 + settings.num_classes)

    io = open(settings.weights_filepath, "r") #do io
        readconstants!(io, nr_constants)
        #First Conv layer
        loadconv!(model.layers[1].layers[1], io, 3, 3, 3, 32)#25
        #Second Conv layer
        loadconv!(model.layers[1].layers[3], io, 3, 3, 32, 64)#37
        #Third Conv layer
        loadconv!(model.layers[1].layers[5], io, 3, 3, 64, 128)#49
        #4th Conv layer
        loadconv!(model.layers[1].layers[6], io, 1, 1, 128, 64)#57
        #5th Conv layer
        loadconv!(model.layers[1].layers[7], io, 3, 3, 64, 128)#65
        #6th Conv layer
        loadconv!(model.layers[1].layers[9], io, 3, 3, 128, 256)#77
	#7th Conv layer
        loadconv!(model.layers[1].layers[10], io, 1, 1, 256, 128)#85
	#8th Conv layer
        loadconv!(model.layers[1].layers[11], io, 3, 3, 128, 256)#93
	#9th Conv layer
        loadconv!(model.layers[1].layers[13], io, 3, 3, 256, 512)#105
	#10th Conv layer
        loadconv!(model.layers[1].layers[14], io, 1, 1, 512, 256)#113
	#11th Conv layer
        loadconv!(model.layers[1].layers[15], io, 3, 3, 256, 512)#121
	#12th Conv layer
        loadconv!(model.layers[1].layers[16], io, 1, 1, 512, 256)#129
	#13th Conv layer
        loadconv!(model.layers[1].layers[17], io, 3, 3, 256, 512)#137
	#14th Conv layer
        loadconv!(model.layers[2].layers[2], io, 3, 3, 512, 1024)#149
	#15th Conv layer
        loadconv!(model.layers[2].layers[3], io, 1, 1, 1024, 512)#157
	#16th Conv layer
        loadconv!(model.layers[2].layers[4], io, 3, 3, 512, 1024)#165
	#17th Conv layer
        loadconv!(model.layers[2].layers[5], io, 1, 1, 1024, 512)#173
	#18th Conv layer
        loadconv!(model.layers[2].layers[6], io, 3, 3, 512, 1024)#181
	#19th Conv layer
        loadconv!(model.layers[2].layers[7], io, 3, 3, 1024, 1024)#192
	#20th Conv layer
        loadconv!(model.layers[2].layers[8], io, 3, 3, 1024, 1024)#200
	#21st Conv layer
        loadconv!(model.layers[3].layers[1], io, 1, 1, 512, 64)#211
	#22nd Conv layer
        loadconv!(model.layers[4].layers[1], io, 3, 3, 1280, 1024)#225
	#23rd Conv layer
        loadconv_last!(model.layers[4].layers[2], io, 1, 1, 1024, 125)#233 no batch norm here
        if GPU >= 0
          end
    # end
end

"""
    loadconv!(c, io, d1, d2, d3, d4)

Loads the io to given convolutional layer and updates it by batch-normalization
"""
function loadconv!(c, io::IOStream, d1, d2, d3, d4)
    read!(io, c.b)
    gamma = Array{Float32}(undef, d4)
    mean = Array{Float32}(undef, d4)
    variance = Array{Float32}(undef, d4)
    read!(io, gamma)
    read!(io, mean)
    read!(io, variance)
    toRead = Array{Float32}(undef, d4 * d3 * d2 * d1)
    read!(io, toRead)
    toRead = reshape(toRead, d1, d2, d3, d4)
    c.w = permutedims(toRead, [2, 1, 3, 4])
    c.w, c.b = updateconv!(c, gamma, mean, variance)
    c.w = flipkernel(c.w)
    if GPU >= 0
        c.w = KnetArray(c.w)
        c.b = KnetArray(c.b)
    end
end

# to be simplified. the weights file does not have batch norm data for the last layer
function loadconv_last!(c, io::IOStream, d1, d2, d3, d4)
    read!(io, c.b)
    toRead = Array{Float32}(undef, d4 * d3 * d2 * d1)
    read!(io, toRead)
    toRead = reshape(toRead, d1, d2, d3, d4)
    c.w = permutedims(toRead, [2, 1, 3, 4])
    c.w = flipkernel(c.w)
    if GPU >= 0
        c.w = KnetArray(c.w)
        c.b = KnetArray(c.b)
    end
end
"""
    readconstants!(io)

Read constant and unnecessary numbers from the io stream
"""
function readconstants!(io::IOStream, nr_constants = 4)
    # major = read(io, Int32)
    # minor = read(io, Int32)
    # revision = read(io, Int32)
    # iseen = read(io, Int32)
	# dummy = read(io, Float32)
	for i = 1:nr_constants
		read(io, Float32)
	end
end

## initial
# io = open(settings.weights_filepath, "r") #do io
# 	readconstants!(io)
# 	#First Conv layer
# 	loadconv!(model.layers[1], io, 3, 3, 3, 32)#25
# 	#Second Conv layer
# 	loadconv!(model.layers[3], io, 3, 3, 32, 64)#37
# 	#Third Conv layer
# 	loadconv!(model.layers[5], io, 3, 3, 64, 128)#49
# 	#4th Conv layer
# 	loadconv!(model.layers[6], io, 1, 1, 128, 64)#57
# 	#5th Conv layer
# 	loadconv!(model.layers[7], io, 3, 3, 64, 128)#65
# 	#6th Conv layer
# 	loadconv!(model.layers[9], io, 3, 3, 128, 256)#77
# 	#7th Conv layer
# 	loadconv!(model.layers[10], io, 1, 1, 256, 128)#85
# 	#8th Conv layer
# 	loadconv!(model.layers[11], io, 3, 3, 128, 256)#93
# 	#9th Conv layer
# 	loadconv!(model.layers[13], io, 3, 3, 256, 512)#105
# 	#10th Conv layer
# 	loadconv!(model.layers[14], io, 1, 1, 512, 256)#113
# 	#11th Conv layer
# 	loadconv!(model.layers[15], io, 3, 3, 256, 512)#121
# 	#12th Conv layer
# 	loadconv!(model.layers[16], io, 1, 1, 512, 256)#129
# 	#13th Conv layer
# 	loadconv!(model.layers[17], io, 3, 3, 256, 512)#137
# 	#14th Conv layer
# 	loadconv!(model.layers[19], io, 3, 3, 512, 1024)#149
# 	#15th Conv layer
# 	loadconv!(model.layers[20], io, 1, 1, 1024, 512)#157
# 	#16th Conv layer
# 	loadconv!(model.layers[21], io, 3, 3, 512, 1024)#165
# 	#17th Conv layer
# 	loadconv!(model.layers[22], io, 1, 1, 1024, 512)#173
# 	#18th Conv layer
# 	loadconv!(model.layers[23], io, 3, 3, 512, 1024)#181
# 	#19th Conv layer
# 	loadconv!(model.layers[24], io, 3, 3, 1024, 1024)#192
# 	#20th Conv layer
# 	loadconv!(model.layers[25], io, 3, 3, 1024, 1024)#200
# 	#21st Conv layer
# 	loadconv!(model.layers[26], io, 1, 1, 512, 64)#211
# 	#22nd Conv layer
# 	loadconv!(model.layers[27], io, 3, 3, 1280, 1024)#225
# 	#23rd Conv layer
# 	loadconv_last!(model.layers[28], io, 1, 1, 1024, 125)#233
# 	if GPU >= 0
# 	  end


text:
include("settings.jl")

makedocs(; sitename=sitename, settings...)

deploydocs(; repo="github.com/GTorlai/PastaQ.jl")


text:
123 == 123

123 == 100 + 23

123 == -123

typeof(true)

typeof(false)

typeof(45 == 45)

23 ≤ -5

4 ≠ 44

1.0 == 1

1.000000000001 == 1

3 < 4 && -15 > -10

3 < 4 || -15 > -10

!(-4 ≥ 4) && (1 > 2 || 1 < 2)

1 < 2 <= 2 < 3 == 3 > 2 >= 1 == 1 < 3 != 5

x = -2.3
-3 ≤ x ≤ 4

# Examples of if-else statements

function check_x(x)
    if 0 ≤ x ≤ 10
        println("x is between 0 and 10")
    elseif !(x < 0)
        println("x is not between 0 and 10, and x is not negative (meaning x > 10)")
    else
        println("x is not between 0 and 10, and x is not not negative (meaning x is negative)")
    end
end

check_x(π)

check_x(10.000001)

check_x(2^63)        # trick question - causes Int64 overflow

x = -3
if x ≥ 0 && sqrt(x) < 10
    # do something
else
    println("Sorry, need a non-negative number with square root strictly less than 10")
end

function fibonacci(n)
    # First special cases
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        # General case n ≥ 2 - initialize and iterate
        Fkold = 0
        Fk = 1
        for k = 2:n
            Fknew = Fk + Fkold
            # Swap previous values
            Fkold = Fk
            Fk = Fknew
        end
        return Fk
    end
end

fibonacci(10)

fibonacci(1000)   # Again be careful with overflow of Int64, but this one is OK

function PE1(n)
    s = 0
    for i = 1:n-1
        if i % 3 == 0 || i % 5 == 0
            s += i
        end
    end
    s
end

PE1(10)

PE1(1000)


text:
using FastTransforms, Compat
using Compat.Test

import FastTransforms: normalizecolumns!, maxcolnorm, sphevaluatepi

@testset "Pointwise evaluation" begin
    n = 256
    A = sphrandn(Float64, n, n);
    normalizecolumns!(A);
    B = sph2fourier(A);

    for θ in (0.123, 0.456)
        S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,1] for ℓ in 1:n)
        SA = sum(sphevaluatepi(θ,ℓ-1,0)*A[ℓ,1] for ℓ in 1:n)
        @test norm(S0-SA) < 1000eps()
        for m in 3:2:n
            S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-2] for ℓ in 1:n)
            SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n)
            @test norm(S0-SA) < 1000eps()
            S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-1] for ℓ in 1:n)
            SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n)
            @test norm(S0-SA) < 1000eps()
        end
    end

    for θ in (0.123, 0.456), m in 2:2:n
        S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-2] for ℓ in 1:n)
        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n)
        @test norm(S1-SA) < 1000eps()
        S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-1] for ℓ in 1:n)
        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n)
        @test norm(S1-SA) < 1000eps()
    end
end