text:
module MatrixFactorizations
using Base, LinearAlgebra, ArrayLayouts
import Base: axes, axes1, getproperty, iterate, tail
import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!,
                        copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym
import LinearAlgebra.BLAS: libblas
import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen,
            qr, axpy!, ldiv!, rdiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular, inv,
            chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, det, logabsdet,
            AbstractQ, _zeros, _cut_B, _ret_size, require_one_based_indexing, checksquare,
            checknonsingular, ipiv2perm, copytri!, issuccess

import Base: getindex, setindex!, *, +, -, ==, <, <=, >,
                >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds

import Base: convert, size, view, unsafe_indices,
                first, last, size, length, unsafe_length, step,
                to_indices, to_index, show, fill!, promote_op,
                MultiplicativeInverses, OneTo, ReshapedArray,
                Array, Matrix, Vector, AbstractArray, AbstractMatrix, AbstractVector,
                               similar, copy, convert, promote_rule, rand,
                            IndexStyle, real, imag, Slice, pointer, unsafe_convert, copyto!

import ArrayLayouts: reflector!, reflectorApply!, materialize!, @_layoutlmul, @_layoutrmul,
                     MemoryLayout, adjointlayout, AbstractQLayout, QRPackedQLayout,
                     QRCompactWYQLayout, AdjQRCompactWYQLayout, QRPackedLayout, AdjQRPackedQLayout

export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, choleskyinv!, choleskyinv

abstract type LayoutQ{T} <: AbstractQ{T} end
@_layoutlmul LayoutQ
@_layoutlmul Adjoint{<:Any,<:LayoutQ}
@_layoutrmul LayoutQ
@_layoutrmul Adjoint{<:Any,<:LayoutQ}

axes(Q::LayoutQ, dim::Integer) = axes(getfield(Q, :factors), dim == 2 ? 1 : dim)
axes(Q::LayoutQ) = axes(Q, 1), axes(Q, 2)

include("ul.jl")
include("qr.jl")
include("ql.jl")
include("rq.jl")
include("choleskyinv.jl")
include("polar.jl")

end #module


text:
module effective_interaction_rotated

# Optimized functions for calculation of Omega and Gamma

function Omega(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    cosξdivξ = cos(ξ)/ξ
    return 3.0/4*((1.0-3.0*cosθpow2) * (sin(ξ)+cosξdivξ)/ξ^2 - (1.0-cosθpow2)*cosξdivξ)
end

function Gamma(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    if ξ<1e-3
        return 3/2*(2/3 + 1/15*ξ^2 *(-2 + cosθpow2) + ξ^4*(1/140-cosθpow2/210))
    else
        sinξ = sin(ξ)/ξ
        return 3/2*((1-cosθpow2)*sinξ + (1-3*cosθpow2) * (cos(ξ)-sinξ)/ξ^2)
    end
end

function OmegaGamma(a, cosθpow2)
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((1.0-3.0*cosθpow2) * (sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - (1.0-cosθpow2)*cosξ/(2*pi*a))
    gamma = 1.5*((1-cosθpow2)*sinξ/(2*pi*a) + (1-3*cosθpow2) * (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function OmegaGamma_orthogonal(a)
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - cosξ/(2*pi*a))
    gamma = 1.5*(sinξ/(2*pi*a) + (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function Omega_orthogonal(a)
    ξ = 2*pi*a
    cosξdξ = cos(ξ)/ξ
    return 3.0/4.0*(-cosξdξ + (sin(ξ)+cosξdξ)/ξ^2)
end

function Gamma_orthogonal(a)
    ξ = 2*a
    sincξ = sinc(ξ)
    return 3.0/2.0*(sincξ + (cospi(ξ)-sincξ)/(pi*ξ)^2)
end

function map_effectiveinteraction(a::Vector{T}, f::Function) where T<:Real
    omega_list = float(T)[]
    gamma_list = float(T)[]
    for ai=a
        omega, gamma = f(ai, e_dipole)
        push!(omega_list, omega)
        push!(gamma_list, gamma)
    end
    return omega_list, gamma_list
end

# Finite symmetric systems

# function triangle_orthogonal(a::Float64)
#     omega_eff = 2*Omega_orthogonal(a)
#     gamma_eff = 2*Gamma_orthogonal(a)
#     return omega_eff, gamma_eff
# end

"""
    effective_interaction_rotated.square_orthogonal(a, Nδ)

Effective Omega and Gamma for a square.

The polarization axis is orthogonal to the square plane.

# Arguments
* `a`: Edge length.
* `Nδ`: Phase shift (Number of atoms in 2π).
"""
function square_orthogonal(a::Real, Nδ::Int)
    @assert 0 <= Nδ < 4
    δ = 2.0*pi*Nδ/4
    Ωcos = 2*Omega_orthogonal(a)*cos(δ) + Omega_orthogonal(sqrt(2.)*a)*cos(2*δ)
    Γcos = 2*Gamma_orthogonal(a)*cos(δ) + Gamma_orthogonal(sqrt(2.)*a)*cos(2*δ)
    return Ωcos, Γcos
end

"""
    effective_interaction_rotated.cube_orthogonal(a, dϕ)

Effective Omega and Gamma for a cube.

The polarization axis is orthogonal to the xy faces.

# Arguments
* `a`: edge length.
* `dϕ`: Phase shift.
"""
function cube_orthogonal(a::Real, dϕ)
    omega_eff, gamma_eff = square_orthogonal(a, 0)
    sqrt2 = sqrt(2.)
    sqrt3 = sqrt(3.)
    Θdiag = atan(sqrt2, 1.)
    omega_eff += (Omega(a, 0.) + 2*Omega(sqrt(2.)*a, pi/4.) + Omega(sqrt3*a, Θdiag))*cos(dϕ)
    gamma_eff += (Gamma(a, 0.) + 2*Gamma(sqrt(2.)*a, pi/4.) + Gamma(sqrt3*a, Θdiag))*cos(dϕ)
    return omega_eff, gamma_eff
end

"""
    effective_interaction_rotated.chain_orthogonal(a, N, dϕ)

Effective Omega and Gamma for an infinite chain.

The polarization axis is orthogonal to the chain and the calculation is
done by adding N spins left and N spins right of a central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
* `dϕ`: Phase shift between neighboring spins.
"""
function chain_orthogonal(a::T, N::Int, dϕ::Real) where T<:Real
    omega_eff::float(T) = 0.
    gamma_eff::float(T) = 0.
    d::float(T) = 0.
    ϕ::float(T) = 0.
    for j=1:N
        d += a
        ϕ += dϕ
        omega_eff += Omega_orthogonal(d)*cos(ϕ)
        gamma_eff += Gamma_orthogonal(d)*cos(ϕ)
    end
    return 2*omega_eff, 2*gamma_eff
end

end # module

ext: jl


text:
module Lexers

include("utilities.jl")

import ..Tokens
import ..Tokens: Token, Kind, UNICODE_OPS, EMPTY_TOKEN, isliteral

import ..Tokens: FUNCTION, ABSTRACT, IDENTIFIER, BAREMODULE, BEGIN, BREAK, CATCH, CONST, CONTINUE,
                 DO, ELSE, ELSEIF, END, EXPORT, FALSE, FINALLY, FOR, FUNCTION, GLOBAL, LET, LOCAL, IF,
                 IMPORT, MACRO, MODULE, OUTER, QUOTE, RETURN, TRUE, TRY, TYPE, USING, WHILE, ISA, IN,
                 MUTABLE, PRIMITIVE, STRUCT, WHERE

export tokenize

@inline ishex(c::Char) = isdigit(c) || ('a' <= c <= 'f') || ('A' <= c <= 'F')
@inline isbinary(c::Char) = c == '0' || c == '1'
@inline isoctal(c::Char) =  '0' ≤ c ≤ '7'
@inline iswhitespace(c::Char) = Base.isspace(c) || c === '\ufeff'

struct StringState
    triplestr::Bool
    raw::Bool
    delim::Char
    paren_depth::Int
end

"""
`Lexer` reads from an input stream and emits a single token each time
`next_token` is called.

Ideally a lexer is stateless but some state is needed here for:
* Disambiguating cases like x' (adjoint) vs 'x' (character literal)
* Tokenizing code within string interpolations
"""
mutable struct Lexer{IO_t <: IO}
    io::IO_t
    io_startpos::Int

    token_start_row::Int
    token_start_col::Int
    token_startpos::Int

    current_row::Int
    current_col::Int
    current_pos::Int

    last_token::Tokens.Kind
    string_states::Vector{StringState}
    charstore::IOBuffer
    chars::Tuple{Char,Char,Char,Char}
    charspos::Tuple{Int,Int,Int,Int}
    doread::Bool
    dotop::Bool
    errored::Bool
end

function Lexer(io::IO)
    c1 = ' '
    p1 = position(io)
    if eof(io)
        c2, p2 = EOF_CHAR, p1
        c3, p3 = EOF_CHAR, p1
        c4, p4 = EOF_CHAR, p1
    else
        c2 = read(io, Char)
        p2 = position(io)
        if eof(io)
            c3, p3 = EOF_CHAR, p1
            c4, p4 = EOF_CHAR, p1
        else
            c3 = read(io, Char)
            p3 = position(io)
            if eof(io)
                c4, p4 = EOF_CHAR, p1
            else
                c4 = read(io, Char)
                p4 = position(io)
            end
        end
    end
    Lexer(io, position(io), 1, 1, position(io), 1, 1, position(io),
                  Tokens.ERROR, Vector{StringState}(), IOBuffer(),
                  (c1,c2,c3,c4), (p1,p2,p3,p4), false, false, false)
end
Lexer(str::AbstractString) = Lexer(IOBuffer(str))

"""
    tokenize(x, T = Token)

Returns an `Iterable` containing the tokenized input. Can be reverted by e.g.
`join(untokenize.(tokenize(x)))`. Setting `T` chooses the type of token
produced by the lexer (`Token` or `Token`).
"""
tokenize(x) = Lexer(x)

# Iterator interface
Base.IteratorSize(::Type{<:Lexer}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:Lexer}) = Base.HasEltype()
Base.eltype(::Type{<:Lexer}) = Token

function Base.iterate(l::Lexer)
    seekstart(l)
    l.token_startpos = position(l)
    l.token_start_row = 1
    l.token_start_col = 1

    l.current_row = 1
    l.current_col = 1
    l.current_pos = l.io_startpos
    t = next_token(l)
    return t, t.kind == Tokens.ENDMARKER
end

function Base.iterate(l::Lexer, isdone::Any)
    isdone && return nothing
    t = next_token(l)
    return t, t.kind == Tokens.ENDMARKER
end

function Base.show(io::IO, l::Lexer)
    print(io, typeof(l), " at position: ", position(l))
end

"""
    startpos(l::Lexer)

Return the latest `Token`'s starting position.
"""
startpos(l::Lexer) = l.token_startpos

"""
    startpos!(l::Lexer, i::Integer)

Set a new starting position.
"""
startpos!(l::Lexer, i::Integer) = l.token_startpos = i

Base.seekstart(l::Lexer) = seek(l.io, l.io_startpos)

"""
    seek2startpos!(l::Lexer)

Sets the lexer's current position to the beginning of the latest `Token`.
"""
seek2startpos!(l::Lexer) = seek(l, startpos(l))

"""
    peekchar(l::Lexer)

Returns the next character without changing the lexer's state.
"""
peekchar(l::Lexer) = l.chars[2]

"""
dpeekchar(l::Lexer)

Returns the next two characters without changing the lexer's state.
"""
dpeekchar(l::Lexer) = l.chars[2], l.chars[3]

"""
peekchar3(l::Lexer)

Returns the next three characters without changing the lexer's state.
"""
peekchar3(l::Lexer) = l.chars[2], l.chars[3], l.chars[4]

"""
    position(l::Lexer)

Returns the current position.
"""
Base.position(l::Lexer) = l.charspos[1]

"""
    eof(l::Lexer)

Determine whether the end of the lexer's underlying buffer has been reached.
"""# Base.position(l::Lexer) = Base.position(l.io)
eof(l::Lexer) = eof(l.io)

Base.seek(l::Lexer, pos) = seek(l.io, pos)

"""
    start_token!(l::Lexer)

Updates the lexer's state such that the next  `Token` will start at the current
position.
"""
function start_token!(l::Lexer)
    l.token_startpos = l.charspos[1]
    l.token_start_row = l.current_row
    l.token_start_col = l.current_col
end

"""
    readchar(l::Lexer)

Returns the next character and increments the current position.
"""
function readchar end

function readchar(l::Lexer)
    c = readchar(l.io)
    l.chars = (l.chars[2], l.chars[3], l.chars[4], c)
    l.charspos = (l.charspos[2], l.charspos[3], l.charspos[4], position(l.io))
    return l.chars[1]
end

"""
    accept(l::Lexer, f::Union{Function, Char, Vector{Char}, String})

Consumes the next character `c` if either `f::Function(c)` returns true, `c == f`
for `c::Char` or `c in f` otherwise. Returns `true` if a character has been
consumed and `false` otherwise.
"""
@inline function accept(l::Lexer, f::Union{Function, Char, Vector{Char}, String})
    c = peekchar(l)
    if isa(f, Function)
        ok = f(c)
    elseif isa(f, Char)
        ok = c == f
    else
        ok = c in f
    end
    ok && readchar(l)
    return ok
end

"""
    accept_batch(l::Lexer, f)

Consumes all following characters until `accept(l, f)` is `false`.
"""
@inline function accept_batch(l::Lexer, f)
    ok = false
    while accept(l, f)
        ok = true
    end
    return ok
end

"""
    emit(l::Lexer, kind::Kind)

Returns a `Token` of kind `kind` with contents `str` and starts a new `Token`.
"""
function emit(l::Lexer, kind::Kind)
    suffix = false
    if optakessuffix(kind)
        while isopsuffix(peekchar(l))
            readchar(l)
            suffix = true
        end
    end

    tok = Token(kind, startpos(l), position(l) - 1, l.dotop, suffix)

    l.dotop = false
    l.last_token = kind
    return tok
end

"""
    emit_error(l::Lexer, err::Kind=Tokens.ERROR)

Returns an `ERROR` token with error `err` and starts a new `Token`.
"""
function emit_error(l::Lexer, err::Kind = Tokens.ERROR)
    l.errored = true
    @assert Tokens.iserror(err)
    return emit(l, err)
end

"""
    next_token(l::Lexer)

Returns the next `Token`.
"""
function next_token(l::Lexer, start = true)
    start && start_token!(l)
    if !isempty(l.string_states)
        lex_string_chunk(l)
    else
        _next_token(l, readchar(l))
    end
end

function _next_token(l::Lexer, c)
    if eof(c)
        return emit(l, Tokens.ENDMARKER)
    elseif iswhitespace(c)
        return lex_whitespace(l, c)
    elseif c == '['
        return emit(l, Tokens.LSQUARE)
    elseif c == ']'
        return emit(l, Tokens.RSQUARE)
    elseif c == '{'
        return emit(l, Tokens.LBRACE)
    elseif c == ';'
        return emit(l, Tokens.SEMICOLON)
    elseif c == '}'
        return emit(l, Tokens.RBRACE)
    elseif c == '('
        return emit(l, Tokens.LPAREN)
    elseif c == ')'
        return emit(l, Tokens.RPAREN)
    elseif c == ','
        return emit(l, Tokens.COMMA)
    elseif c == '*'
        return lex_star(l);
    elseif c == '^'
        return lex_circumflex(l);
    elseif c == '@'
        return emit(l, Tokens.AT_SIGN)
    elseif c == '?'
        return emit(l, Tokens.CONDITIONAL)
    elseif c == '$'
        return lex_dollar(l);
    elseif c == '⊻'
        return lex_xor(l);
    elseif c == '~'
        return emit(l, Tokens.APPROX)
    elseif c == '#'
        return lex_comment(l)
    elseif c == '='
        return lex_equal(l)
    elseif c == '!'
        return lex_exclaim(l)
    elseif c == '>'
        return lex_greater(l)
    elseif c == '<'
        return lex_less(l)
    elseif c == ':'
        return lex_colon(l)
    elseif c == '|'
        return lex_bar(l)
    elseif c == '&'
        return lex_amper(l)
    elseif c == '\''
        return lex_prime(l)
    elseif c == '÷'
        return lex_division(l)
    elseif c == '"'
        return lex_quote(l);
    elseif c == '%'
        return lex_percent(l);
    elseif c == '/'
        return lex_forwardslash(l);
    elseif c == '\\'
        return lex_backslash(l);
    elseif c == '.'
        return lex_dot(l);
    elseif c == '+'
        return lex_plus(l);
    elseif c == '-'
        return lex_minus(l);
    elseif c == '−' # \minus '−' treated as hyphen '-'
        return emit(l, accept(l, '=') ? Tokens.MINUS_EQ : Tokens.MINUS)
    elseif c == '`'
        return lex_backtick(l);
    elseif is_identifier_start_char(c)
        return lex_identifier(l, c)
    elseif isdigit(c)
        return lex_digit(l, Tokens.INTEGER)
    elseif (k = get(UNICODE_OPS, c, Tokens.ERROR)) != Tokens.ERROR
        return emit(l, k)
    else
        emit_error(l)
    end
end

# We're inside a string; possibly reading the string characters, or maybe in
# Julia code within an interpolation.
function lex_string_chunk(l)
    state = last(l.string_states)
    if state.paren_depth > 0
        # Read normal Julia code inside an interpolation but track nesting of
        # parentheses.
        c = readchar(l)
        if c == '('
            l.string_states[end] = StringState(state.triplestr, state.raw, state.delim,
                                               state.paren_depth + 1)
            return emit(l, Tokens.LPAREN)
        elseif c == ')'
            l.string_states[end] = StringState(state.triplestr, state.raw, state.delim,
                                               state.paren_depth - 1)
            return emit(l, Tokens.RPAREN)
        else
            return _next_token(l, c)
        end
    end
    pc = peekchar(l)
    if l.last_token == Tokens.EX_OR
        pc = peekchar(l)
        # Interpolated symbol or expression
        if pc == '('
            readchar(l)
            l.string_states[end] = StringState(state.triplestr, state.raw, state.delim,
                                               state.paren_depth + 1)
            return emit(l, Tokens.LPAREN)
        elseif is_identifier_start_char(pc)
            return lex_identifier(l, readchar(l))
        else
            # Getting here is a syntax error - fall through to reading string
            # characters and let the parser deal with it.
        end
    elseif l.last_token == Tokens.IDENTIFIER &&
            !(eof(pc) || is_operator_start_char(pc) || is_never_id_char(pc))
        # Only allow certain characters after interpolated vars
        # https://github.com/JuliaLang/julia/pull/25234
        return emit_error(l, Tokens.INVALID_INTERPOLATION_TERMINATOR)
    end
    if eof(pc)
        return emit(l, Tokens.ENDMARKER)
    elseif !state.raw && pc == '$'
        # Start interpolation
        readchar(l)
        return emit(l, Tokens.EX_OR)
    elseif !state.raw && pc == '\\' && (pc2 = dpeekchar(l)[2];
                                        pc2 == '\r' || pc2 == '\n')
        # Process escaped newline as whitespace
        readchar(l)
        readchar(l)
        if pc2 == '\r' && peekchar(l) == '\n'
            readchar(l)
        end
        while (pc = peekchar(l); pc == ' ' || pc == '\t')
            readchar(l)
        end
        return emit(l, Tokens.WHITESPACE)
    elseif pc == state.delim && string_terminates(l, state.delim, state.triplestr)
        # Terminate string
        pop!(l.string_states)
        readchar(l)
        if state.triplestr
            readchar(l); readchar(l)
            return emit(l, state.delim == '"' ?
                        Tokens.TRIPLE_DQUOTE : Tokens.TRIPLE_BACKTICK)
        else
            return emit(l, state.delim == '"' ? Tokens.DQUOTE : Tokens.BACKTICK)
        end
    end
    # Read a chunk of string characters
    if state.raw
        # Raw strings treat all characters as literals with the exception that
        # the closing quotes can be escaped with an odd number of \ characters.
        while true
            pc = peekchar(l)
            if string_terminates(l, state.delim, state.triplestr) || eof(pc)
                break
            elseif state.triplestr && (pc == '\n' || pc == '\r')
                # triple quoted newline splitting
                readchar(l)
                if pc == '\r' && peekchar(l) == '\n'
                    readchar(l)
                end
                break
            end
            c = readchar(l)
            if c == '\\'
                n = 1
                while peekchar(l) == '\\'
                    readchar(l)
                    n += 1
                end
                if peekchar(l) == state.delim && !iseven(n)
                    readchar(l)
                end
            end
        end
    else
        while true
            pc = peekchar(l)
            if pc == '$' || eof(pc)
                break
            elseif state.triplestr && (pc == '\n' || pc == '\r')
                # triple quoted newline splitting
                readchar(l)
                if pc == '\r' && peekchar(l) == '\n'
                    readchar(l)
                end
                break
            elseif pc == state.delim && string_terminates(l, state.delim, state.triplestr)
                break
            elseif pc == '\\'
                # Escaped newline
                pc2 = dpeekchar(l)[2]
                if pc2 == '\r' || pc2 == '\n'
                    break
                end
            end
            c = readchar(l)
            if c == '\\'
                c = readchar(l)
                eof(c) && break
                continue
            end
        end
    end
    return emit(l, state.delim == '"' ?  Tokens.STRING : Tokens.CMD)
end

# Lex whitespace, a whitespace char `c` has been consumed
function lex_whitespace(l::Lexer, c)
    k = Tokens.WHITESPACE
    while true
        if c == '\n'
            k = Tokens.NEWLINE_WS
        end
        pc = peekchar(l)
        # stop on non whitespace and limit to a single newline in a token
        if !iswhitespace(pc) || (k == Tokens.NEWLINE_WS && pc == '\n')
            break
        end
        c = readchar(l)
    end
    return emit(l, k)
end

function lex_comment(l::Lexer, doemit=true)
    if peekchar(l) != '='
        while true
            pc = peekchar(l)
            if pc == '\n' || eof(pc)
                return doemit ? emit(l, Tokens.COMMENT) : EMPTY_TOKEN
            end
            readchar(l)
        end
    else
        pc = '#'
        c = readchar(l) # consume the '='
        n_start, n_end = 1, 0
        while true
            if eof(c)
                return doemit ? emit_error(l, Tokens.EOF_MULTICOMMENT) : EMPTY_TOKEN
            end
            nc = readchar(l)
            if c == '#' && nc == '='
                n_start += 1
            elseif c == '=' && nc == '#' && pc != '#'
                n_end += 1
            end
            if n_start == n_end
                return doemit ? emit(l, Tokens.COMMENT) : EMPTY_TOKEN
            end
            pc = c
            c = nc
        end
    end
end

# Lex a greater char, a '>' has been consumed
function lex_greater(l::Lexer)
    if accept(l, '>') # >>
        if accept(l, '>') # >>>
            if accept(l, '=') # >>>=
                return emit(l, Tokens.UNSIGNED_BITSHIFT_EQ)
            else # >>>?, ? not a =
                return emit(l, Tokens.UNSIGNED_BITSHIFT)
            end
        elseif accept(l, '=') # >>=
            return emit(l, Tokens.RBITSHIFT_EQ)
        else # '>>'
            return emit(l, Tokens.RBITSHIFT)
        end
    elseif accept(l, '=') # >=
        return emit(l, Tokens.GREATER_EQ)
    elseif accept(l, ':') # >:
        return emit(l, Tokens.ISSUPERTYPE)
    else  # '>'
        return emit(l, Tokens.GREATER)
    end
end

# Lex a less char, a '<' has been consumed
function lex_less(l::Lexer)
    if accept(l, '<') # <<
        if accept(l, '=') # <<=
            return emit(l, Tokens.LBITSHIFT_EQ)
        else # '<<?', ? not =, ' '
            return emit(l, Tokens.LBITSHIFT)
        end
    elseif accept(l, '=') # <=
        return emit(l, Tokens.LESS_EQ)
    elseif accept(l, ':')
        return emit(l, Tokens.ISSUBTYPE)
    elseif accept(l, '|') # <|
        return emit(l, Tokens.LPIPE)
    elseif dpeekchar(l) == ('-', '-') # <-- or <-->
        readchar(l); readchar(l)
        if accept(l, '>')
            return emit(l, Tokens.DOUBLE_ARROW)
        else
            return emit(l, Tokens.LEFT_ARROW)
        end
    else
        return emit(l, Tokens.LESS) # '<'
    end
end

# Lex all tokens that start with an = character.
# An '=' char has been consumed
function lex_equal(l::Lexer)
    if accept(l, '=') # ==
        if accept(l, '=') # ===
            emit(l, Tokens.EQEQEQ)
        else
            emit(l, Tokens.EQEQ)
        end
    elseif accept(l, '>') # =>
        emit(l, Tokens.PAIR_ARROW)
    else
        emit(l, Tokens.EQ)
    end
end

# Lex a colon, a ':' has been consumed
function lex_colon(l::Lexer)
    if accept(l, ':') # '::'
        return emit(l, Tokens.DECLARATION)
    elseif accept(l, '=') # ':='
        return emit(l, Tokens.COLON_EQ)
    else
        return emit(l, Tokens.COLON)
    end
end

function lex_exclaim(l::Lexer)
    if accept(l, '=') # !=
        if accept(l, '=') # !==
            return emit(l, Tokens.NOT_IS)
        else # !=
            return emit(l, Tokens.NOT_EQ)
        end
    else
        return emit(l, Tokens.NOT)
    end
end

function lex_percent(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.REM_EQ)
    else
        return emit(l, Tokens.REM)
    end
end

function lex_bar(l::Lexer)
    if accept(l, '=') # |=
        return emit(l, Tokens.OR_EQ)
    elseif accept(l, '>') # |>
        return emit(l, Tokens.RPIPE)
    elseif accept(l, '|') # ||
        return emit(l, Tokens.LAZY_OR)
    else
        emit(l, Tokens.OR) # '|'
    end
end

function lex_plus(l::Lexer)
    if accept(l, '+')
        return emit(l, Tokens.PLUSPLUS)
    elseif accept(l, '=')
        return emit(l, Tokens.PLUS_EQ)
    end
    return emit(l, Tokens.PLUS)
end

function lex_minus(l::Lexer)
    if accept(l, '-')
        if accept(l, '>')
            return emit(l, Tokens.RIGHT_ARROW)
        else
            return emit_error(l, Tokens.INVALID_OPERATOR) # "--" is an invalid operator
        end
    elseif accept(l, '>')
        return emit(l, Tokens.ANON_FUNC)
    elseif accept(l, '=')
        return emit(l, Tokens.MINUS_EQ)
    end
    return emit(l, Tokens.MINUS)
end

function lex_star(l::Lexer)
    if accept(l, '*')
        return emit_error(l, Tokens.INVALID_OPERATOR) # "**" is an invalid operator use ^
    elseif accept(l, '=')
        return emit(l, Tokens.STAR_EQ)
    end
    return emit(l, Tokens.STAR)
end

function lex_circumflex(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.CIRCUMFLEX_EQ)
    end
    return emit(l, Tokens.CIRCUMFLEX_ACCENT)
end

function lex_division(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.DIVISION_EQ)
    end
    return emit(l, Tokens.DIVISION_SIGN)
end

function lex_dollar(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.EX_OR_EQ)
    end
    return emit(l, Tokens.EX_OR)
end

function lex_xor(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.XOR_EQ)
    end
    return emit(l, Tokens.XOR)
end

function accept_number(l::Lexer, f::F) where {F}
    while true
        pc, ppc = dpeekchar(l)
        if pc == '_' && !f(ppc)
            return
        elseif f(pc) || pc == '_'
            readchar(l)
        else
            return
        end
    end
end

# A digit has been consumed
function lex_digit(l::Lexer, kind)
    accept_number(l, isdigit)
    pc,ppc = dpeekchar(l)
    if pc == '.'
        if kind === Tokens.FLOAT
            # If we enter the function with kind == FLOAT then a '.' has been parsed.
            readchar(l)
            return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
        elseif ppc == '.'
            return emit(l, kind)
        elseif is_operator_start_char(ppc) && ppc !== ':'
            readchar(l)
            return emit_error(l)
        elseif (!(isdigit(ppc) ||
            iswhitespace(ppc) ||
            is_identifier_start_char(ppc)
            || ppc == '('
            || ppc == ')'
            || ppc == '['
            || ppc == ']'
            || ppc == '{'
            || ppc == '}'
            || ppc == ','
            || ppc == ';'
            || ppc == '@'
            || ppc == '`'
            || ppc == '"'
            || ppc == ':'
            || ppc == '?'
            || eof(ppc)))
            kind = Tokens.INTEGER

            return emit(l, kind)
        end
        readchar(l)

        kind = Tokens.FLOAT
        accept_number(l, isdigit)
        pc, ppc = dpeekchar(l)
        if (pc == 'e' || pc == 'E' || pc == 'f') && (isdigit(ppc) || ppc == '+' || ppc == '-' || ppc == '−')
            kind = Tokens.FLOAT
            readchar(l)
            accept(l, "+-−")
            if accept_batch(l, isdigit)
                pc,ppc = dpeekchar(l)
                if pc === '.' && !dotop2(ppc, ' ')
                    accept(l, '.')
                    return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
                end
            else
                return emit_error(l)
            end
        elseif pc == '.' && (is_identifier_start_char(ppc) || eof(ppc))
            readchar(l)
            return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
        end

    elseif (pc == 'e' || pc == 'E' || pc == 'f') && (isdigit(ppc) || ppc == '+' || ppc == '-' || ppc == '−')
        kind = Tokens.FLOAT
        readchar(l)
        accept(l, "+-−")
        if accept_batch(l, isdigit)
            pc,ppc = dpeekchar(l)
            if pc === '.' && !dotop2(ppc, ' ')
                accept(l, '.')
                return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
            end
        else
            return emit_error(l)
        end
    elseif position(l) - startpos(l) == 1 && l.chars[1] == '0'
        kind == Tokens.INTEGER
        if pc == 'x'
            kind = Tokens.HEX_INT
            isfloat = false
            readchar(l)
            !(ishex(ppc) || ppc == '.') && return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
            accept_number(l, ishex)
            if accept(l, '.')
                accept_number(l, ishex)
                isfloat = true
            end
            if accept(l, "pP")
                kind = Tokens.FLOAT
                accept(l, "+-−")
                accept_number(l, isdigit)
            elseif isfloat
                return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
            end
        elseif pc == 'b'
            !isbinary(ppc) && return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
            readchar(l)
            accept_number(l, isbinary)
            kind = Tokens.BIN_INT
        elseif pc == 'o'
            !isoctal(ppc) && return emit_error(l, Tokens.INVALID_NUMERIC_CONSTANT)
            readchar(l)
            accept_number(l, isoctal)
            kind = Tokens.OCT_INT
        end
    end
    return emit(l, kind)
end

function lex_prime(l, doemit = true)
    if l.last_token == Tokens.IDENTIFIER ||
        Tokens.iscontextualkeyword(l.last_token) ||
        Tokens.iswordoperator(l.last_token) ||
        l.last_token == Tokens.DOT ||
        l.last_token ==  Tokens.RPAREN ||
        l.last_token ==  Tokens.RSQUARE ||
        l.last_token ==  Tokens.RBRACE ||
        l.last_token == Tokens.PRIME ||
        l.last_token == Tokens.END || isliteral(l.last_token)
        return emit(l, Tokens.PRIME)
    else
        if accept(l, '\'')
            if accept(l, '\'')
                return doemit ? emit(l, Tokens.CHAR) : EMPTY_TOKEN
            else
                # Empty char literal
                # Arguably this should be an error here, but we generally
                # look at the contents of the char literal in the parser,
                # so we defer erroring until there.
                return doemit ? emit(l, Tokens.CHAR) : EMPTY_TOKEN
            end
        end
        while true
            c = readchar(l)
            if eof(c)
                return doemit ? emit_error(l, Tokens.EOF_CHAR) : EMPTY_TOKEN
            elseif c == '\\'
                if eof(readchar(l))
                    return doemit ? emit_error(l, Tokens.EOF_CHAR) : EMPTY_TOKEN
                end
            elseif c == '\''
                return doemit ? emit(l, Tokens.CHAR) : EMPTY_TOKEN
            end
        end
    end
end

function lex_amper(l::Lexer)
    if accept(l, '&')
        return emit(l, Tokens.LAZY_AND)
    elseif accept(l, '=')
        return emit(l, Tokens.AND_EQ)
    else
        return emit(l, Tokens.AND)
    end
end

# Parse a token starting with a quote.
# A '"' has been consumed
function lex_quote(l::Lexer)
    raw = l.last_token == Tokens.IDENTIFIER ||
          Tokens.iscontextualkeyword(l.last_token) ||
          Tokens.iswordoperator(l.last_token)
    pc, dpc = dpeekchar(l)
    triplestr = pc == '"' && dpc == '"'
    push!(l.string_states, StringState(triplestr, raw, '"', 0))
    if triplestr
        readchar(l)
        readchar(l)
        emit(l, Tokens.TRIPLE_DQUOTE)
    else
        emit(l, Tokens.DQUOTE)
    end
end

function string_terminates(l, delim::Char, triplestr::Bool)
    if triplestr
        c1, c2, c3 = peekchar3(l)
        c1 === delim && c2 === delim && c3 === delim
    else
        peekchar(l) === delim
    end
end

# Parse a token starting with a forward slash.
# A '/' has been consumed
function lex_forwardslash(l::Lexer)
    if accept(l, '/') # //
        if accept(l, '=') # //=
            return emit(l, Tokens.FWDFWD_SLASH_EQ)
        else
            return emit(l, Tokens.FWDFWD_SLASH)
        end
    elseif accept(l, '=') # /=
        return emit(l, Tokens.FWD_SLASH_EQ)
    else
        return emit(l, Tokens.FWD_SLASH)
    end
end

function lex_backslash(l::Lexer)
    if accept(l, '=')
        return emit(l, Tokens.BACKSLASH_EQ)
    end
    return emit(l, Tokens.BACKSLASH)
end

# TODO .op
function lex_dot(l::Lexer)
    if accept(l, '.')
        if accept(l, '.')
            return emit(l, Tokens.DDDOT)
        else
            return emit(l, Tokens.DDOT)
        end
    elseif Base.isdigit(peekchar(l))
        return lex_digit(l, Tokens.FLOAT)
    else
        pc, dpc = dpeekchar(l)
        if dotop1(pc)
            l.dotop = true
            return _next_token(l, readchar(l))
        elseif pc =='+'
            l.dotop = true
            readchar(l)
            return lex_plus(l)
        elseif pc =='-'
            l.dotop = true
            readchar(l)
            return lex_minus(l)
        elseif pc == '−'
            l.dotop = true
            readchar(l)
            return emit(l, accept(l, '=') ? Tokens.MINUS_EQ : Tokens.MINUS)
        elseif pc =='*'
            l.dotop = true
            readchar(l)
            return lex_star(l)
        elseif pc =='/'
            l.dotop = true
            readchar(l)
            return lex_forwardslash(l)
        elseif pc =='\\'
            l.dotop = true
            readchar(l)
            return lex_backslash(l)
        elseif pc =='^'
            l.dotop = true
            readchar(l)
            return lex_circumflex(l)
        elseif pc =='<'
            l.dotop = true
            readchar(l)
            return lex_less(l)
        elseif pc =='>'
            l.dotop = true
            readchar(l)
            return lex_greater(l)
        elseif pc =='&'
            l.dotop = true
            readchar(l)
            if accept(l, '=')
                return emit(l, Tokens.AND_EQ)
            else
                if accept(l, '&')
                    return emit(l, Tokens.LAZY_AND)
                end
                return emit(l, Tokens.AND)
            end
        elseif pc =='%'
            l.dotop = true
            readchar(l)
            return lex_percent(l)
        elseif pc == '=' && dpc != '>'
            l.dotop = true
            readchar(l)
            return lex_equal(l)
        elseif pc == '|'
            l.dotop = true
            readchar(l)
            if accept(l, '|')
                return emit(l, Tokens.LAZY_OR)
            end
            return lex_bar(l)
        elseif pc == '!' && dpc == '='
            l.dotop = true
            readchar(l)
            return lex_exclaim(l)
        elseif pc == '⊻'
            l.dotop = true
            readchar(l)
            return lex_xor(l)
        elseif pc == '÷'
            l.dotop = true
            readchar(l)
            return lex_division(l)
        elseif pc == '=' && dpc == '>'
            l.dotop = true
            readchar(l)
            return lex_equal(l)
        end
        return emit(l, Tokens.DOT)
    end
end

# A ` has been consumed
function lex_backtick(l::Lexer)
    pc, dpc = dpeekchar(l)
    triplestr = pc == '`' && dpc == '`'
    # Backticks always contain raw strings only. See discussion on bug
    # https://github.com/JuliaLang/julia/issues/3150
    raw = true
    push!(l.string_states, StringState(triplestr, raw, '`', 0))
    if triplestr
        readchar(l)
        readchar(l)
        emit(l, Tokens.TRIPLE_BACKTICK)
    else
        emit(l, Tokens.BACKTICK)
    end
end

const MAX_KW_LENGTH = 10
function lex_identifier(l::Lexer, c)
    h = simple_hash(c, UInt64(0))
    n = 1
    while true
        pc, ppc = dpeekchar(l)
        if (pc == '!' && ppc == '=') || !is_identifier_char(pc)
            break
        end
        c = readchar(l)
        h = simple_hash(c, h)
        n += 1
    end

    if n > MAX_KW_LENGTH
        emit(l, IDENTIFIER)
    else
        emit(l, get(kw_hash, h, IDENTIFIER))
    end
end

# This creates a hash for chars in [a-z] using 5 bit per char.
# Requires an additional input-length check somewhere, because
# this only works up to ~12 chars.
@inline function simple_hash(c::Char, h::UInt64)
    bytehash = (clamp(c - 'a' + 1, -1, 30) % UInt8) & 0x1f
    h << 5 + bytehash
end

function simple_hash(str)
    ind = 1
    h = UInt64(0)
    while ind <= length(str)
        h = simple_hash(str[ind], h)
        ind = nextind(str, ind)
    end
    h
end

kws = [
Tokens.BAREMODULE,
Tokens.BEGIN,
Tokens.BREAK,
Tokens.CATCH,
Tokens.CONST,
Tokens.CONTINUE,
Tokens.DO,
Tokens.ELSE,
Tokens.ELSEIF,
Tokens.END,
Tokens.EXPORT,
Tokens.FINALLY,
Tokens.FOR,
Tokens.FUNCTION,
Tokens.GLOBAL,
Tokens.IF,
Tokens.IMPORT,
Tokens.LET,
Tokens.LOCAL,
Tokens.MACRO,
Tokens.MODULE,
Tokens.QUOTE,
Tokens.RETURN,
Tokens.STRUCT,
Tokens.TRY,
Tokens.USING,
Tokens.WHILE,
Tokens.IN,
Tokens.ISA,
Tokens.WHERE,
Tokens.TRUE,
Tokens.FALSE,

Tokens.ABSTRACT,
Tokens.AS,
Tokens.DOC,
Tokens.MUTABLE,
Tokens.OUTER,
Tokens.PRIMITIVE,
Tokens.TYPE,
Tokens.VAR,
]

const kw_hash = Dict(simple_hash(lowercase(string(kw))) => kw for kw in kws)

end # module


text:
# Процессы для которых требуется указать только число точек

poisson(n; dim=2) = rand(dim, n)

function circle_uniform(n)
    ws = range(0., length=n, step=2π / n)
    [f(w) for f in (cos, sin), w in ws]
end

function circle_poisson(n)
    ws = rand(n) * 2π
    [f(w) for f in (cos, sin), w in ws]
end

normal(n; dim=2) = rand(Normal(), dim, n)


text:
################################################################################

"""
Keep track of scores & indexes

# Fields
    score::Float64: score
    index::Int64: index

See also: [`HMM`](@ref), [`HMMParams`](@ref)
"""
struct ScorePair
  score::Float64
  index::Int64
end

"""
Hidden Markov model object

# Fields
    dataM::Array: hold data
    tbM::Array: collect traceback

See also: [`HMMParams`](@ref), [`setup!`](@ref),
"""
mutable struct HMM
  dataM::Array
  tbM::Array
end

"""
Hidden Markov model parameters

# Fields
    pen::Float64: penalty for hidden Markov model
    distance::Function: function used to calculate distance

See also: [`HMM`](@ref), [`setup!`](@ref),
"""
@with_kw struct HMMParams
  pen::Float64                 = 200.      # penalty
  distance::Function           = distance  # distance function
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
using Test, StructuralDynamicsODESolvers

include("problems.jl")
include("api.jl")
include("solvers.jl")


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

function read_geant4_ascii end
export read_geant4_ascii

struct Geant4CSVInput <: AbstractLegendInput
    stream::BufferedInputStream
end

export Geant4CSVInput

Base.open(filename::AbstractString, ::Type{Geant4CSVInput}) =
    Geant4CSVInput(BufferedInputStream(open(filename)))

Base.close(input::Geant4CSVInput) = close(input.stream)

function Base.getindex(input::Geant4CSVInput, ::Colon)
    read(input)
end

function Base.read(input::Geant4CSVInput)
    input_stream = input.stream
    title_expr = r"""^#title (.+)$"""
    fsep_expr = r"""^#separator ([0-9]+)$"""
    vsep_expr = r"""^#vector_separator ([0-9]+)$"""
    colspec_expr = r"""^#column ([^ ]+) ([^ ]+)$"""

    readline(input_stream) == "#class tools::wcsv::ntuple" || throw(ErrorException("Input doesn't seem to be Geant4 ntuple CSV data"))

    title = _match_or_throw(title_expr, readline(input_stream))[1]
    separator = Char(parse(Int, _match_or_throw(fsep_expr, readline(input_stream))[1]))
    vector_separator = Char(parse(Int, _match_or_throw(vsep_expr, readline(input_stream))[1]))

    @debug "File header information: " title separator vector_separator

    cols = Dict{String,Int}()
    ncols::Int = 0
    reading_header::Bool = true
    while reading_header
        c = Char(BufferedStreams.peek(input_stream))
        if c == '#'
            colname = _match_or_throw(colspec_expr, readline(input_stream))[2]
            ncols += 1
            cols[colname] = ncols
        else
            reading_header = false
        end
    end
    @assert ncols == length(cols)

    @debug "Columns in file: " cols length(cols)

    colno_evtno = cols["Event#"]
    colno_detno = get(cols, "Detector_ID", 0)
    colno_thit = get(cols, "Time", 0)
    colno_edep = cols["EnergyDeposit"]
    colno_x = cols["X"]
    colno_y = cols["Y"]
    colno_z = cols["Z"]

    @debug "Colums numbers:" colno_evtno colno_detno colno_thit colno_edep colno_x colno_y colno_z

    colno_xyz = colno_x:colno_x+2
    colno_y == colno_xyz[2] && colno_z == colno_xyz[3] || throw(ErrorException("Expected columns Y and Z directly after X"))

    @debug readline(input_stream)

    data = readdlm(input_stream, separator, Float64)::Matrix{Float64}
    data_ncols = size(data, 2)
    data_ncols == ncols || throw(ErrorException("Expected $ncols columns, but data has $data_ncols columns"))

    evtno = Int32.(data[:, colno_evtno])
    edep = Float32.(data[:, colno_edep]) .* u"keV"
    pos_mat = Float32.(data[:, colno_xyz]) .* u"mm"
    pos = nestedview(Array(pos_mat'), SVector{3})

    detno = colno_detno != 0 ? Int32.(data[:, colno_detno]) : fill!(similar(evtno), 1)
    thit = (colno_thit != 0 ? Float32.(data[:, colno_thit]) : fill!(similar(evtno, Float32), NaN)) .* u"s"

    hits = DetectorHits((
        evtno = evtno,
        detno = detno,
        thit = thit,
        edep = edep,
        pos = pos,
    ))

    events = group_by_evtno(hits)
end


text:
using PlutoProfile
using Test

@testset "PlutoProfile.jl" begin
    # Write your tests here.
end


text:
function get_a(obj::OGMFile)
    return jcall(obj, "getA", jdouble, ())
end

function get_anomaly(obj::OGMFile)
    return jcall(obj, "getAnomaly", jdouble, ())
end

function get_anomaly_type(obj::OGMFile)
    return jcall(obj, "getAnomalyType", PositionAngle, ())
end

function get_cov_ref_frame(obj::OGMFile)
    return jcall(obj, "getCovRefFrame", Frame, ())
end

function get_cov_ref_lof_type(obj::OGMFile)
    return jcall(obj, "getCovRefLofType", LOFType, ())
end

function get_covariance_comment(obj::OGMFile)
    return jcall(obj, "getCovarianceComment", List, ())
end

function get_covariance_matrix(obj::OGMFile)
    return jcall(obj, "getCovarianceMatrix", RealMatrix, ())
end

function get_drag_area(obj::OGMFile)
    return jcall(obj, "getDragArea", jdouble, ())
end

function get_drag_coeff(obj::OGMFile)
    return jcall(obj, "getDragCoeff", jdouble, ())
end

function get_e(obj::OGMFile)
    return jcall(obj, "getE", jdouble, ())
end

function get_epoch(obj::OGMFile)
    return jcall(obj, "getEpoch", AbsoluteDate, ())
end

function get_epoch_comment(obj::OGMFile)
    return jcall(obj, "getEpochComment", List, ())
end

function get_i(obj::OGMFile)
    return jcall(obj, "getI", jdouble, ())
end

function get_keplerian_elements_comment(obj::OGMFile)
    return jcall(obj, "getKeplerianElementsComment", List, ())
end

function get_mass(obj::OGMFile)
    return jcall(obj, "getMass", jdouble, ())
end

function get_meta_data(obj::OGMFile)
    return jcall(obj, "getMetaData", ODMMetaData, ())
end

function get_pa(obj::OGMFile)
    return jcall(obj, "getPa", jdouble, ())
end

function get_raan(obj::OGMFile)
    return jcall(obj, "getRaan", jdouble, ())
end

function get_solar_rad_area(obj::OGMFile)
    return jcall(obj, "getSolarRadArea", jdouble, ())
end

function get_solar_rad_coeff(obj::OGMFile)
    return jcall(obj, "getSolarRadCoeff", jdouble, ())
end

function get_spacecraft_comment(obj::OGMFile)
    return jcall(obj, "getSpacecraftComment", List, ())
end

function get_user_defined_parameters(obj::OGMFile)
    return jcall(obj, "getUserDefinedParameters", Map, ())
end

function has_covariance_matrix(obj::OGMFile)
    return jcall(obj, "hasCovarianceMatrix", jboolean, ())
end

function has_keplerian_elements(obj::OGMFile)
    return jcall(obj, "hasKeplerianElements", jboolean, ())
end


text:
using MeasureTheory
export RadioLikelihood, logdensity
using LinearAlgebra

struct RadioLikelihood{T} <: MeasureBase.AbstractMeasure
    lklhds::T
    function RadioLikelihood(data...)
        ls = map(makelikelihood, data)
        new{typeof(ls)}(ls)
    end
end

function Base.show(io::IO, d::RadioLikelihood{T}) where {T}
    println(io, "RadioLikelihood{$T}")
    println(io, "\tNumber of data products: ", length(d.lklhds))
end

@inline MeasureBase.logdensity(d::RadioLikelihood, m) = sum(x->logdensity(x, m), d.lklhds)

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTVisibilityAmplitudeDatum})
    u,v = getdata(data, :u), getdata(data, :v)
    σ = Diagonal(getdata(data, :error))
    f(m) = amplitude.(Ref(m), u, v)
    k = kernel(MvNormal{(:μ, :σ)},
                    σ = x->σ,
                    μ = x->f(x)
              )
    amp = getdata(data, :amp)
    return Likelihood(k, amp)
end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTLogClosureAmplitudeDatum})
    u1, v1 = getdata(data, :u1), getdata(data, :v1)
    u2, v2 = getdata(data, :u2), getdata(data, :v2)
    u3, v3 = getdata(data, :u3), getdata(data, :v3)
    u4, v4 = getdata(data, :u4), getdata(data, :v4)
    τ = Diagonal(getdata(data, :error))
    f(m) = logclosure_amplitudes(m, u1, v1, u2, v2, u3, v3, u4, v4)
    k = kernel(MvNormal{(:μ, :σ)},
                σ = x->Diagonal(τ),
                μ = x->f(x)
              )

    amp = getdata(data, :amp)
    return Likelihood(k, amp)
end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTClosurePhaseDatum})
    u1, v1 = getdata(data, :u1), getdata(data, :v1)
    u2, v2 = getdata(data, :u2), getdata(data, :v2)
    u3, v3 = getdata(data, :u3), getdata(data, :v3)
    τ = inv.(getdata(data, :error)).^2
    f(m) = closure_phases(m, u1, v1, u2, v2, u3, v3)
    d = CPVonMises{(:μ, :κ)}
    k = kernel(d,
                κ = x->τ,
                μ = x->f(x)
              )

    amp = getdata(data, :phase)
    return Likelihood(k, amp)
end


text:
# ---
# title: 220. Contains Duplicate III
# id: problem220
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort, Ordered Map
# link: <https://leetcode.com/problems/contains-duplicate-iii/description/>
# hidden: true
# ---
#
# Given an array of integers, find out whether there are two distinct indices
# _i_ and _j_ in the array such that the **absolute** difference between
# **nums[i]** and **nums[j]** is at most _t_ and the **absolute** difference
# between _i_ and _j_ is at most _k_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,1], k = 3, t = 0
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,0,1,1], k = 1, t = 2
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,5,9,1,5,9], k = 2, t = 3
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 2 * 104`
#   * `-231 <= nums[i] <= 231 - 1`
#   * `0 <= k <= 104`
#   * `0 <= t <= 231 - 1`
#
#
## @lc code=start
using LeetCode

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


text:
#=Packages needed
LsqFit
Plots
DelimitedFiles
DSP
Statistics
=#

include("load_convolve_data.jl")
include("DE-infinite-TR.jl")
using LsqFit
using Plots

IRF_filename = "IRF_10mm_SDS_750nm.asc"
DTOF_filename = "IL5_10mm_SDS_750nm_20ml.asc"
IRF = loadData(IRF_filename)
DTOF = loadData1(DTOF_filename)

function fitDTOF(DTOF,IRF)

    maxval,maxin = findmax(DTOF.counts)

    function conv_DT_IRF(t,β)
        #function we are fitting raw data to on line 30
        ρ=1
        Rt = (time = t, counts = refl_DT1(t,β,ρ))
        convDT = convolveTR_IRF(Rt,IRF)
        ~,maxinConv = findmax(convDT.counts)
        timefit = convDT.time .- convDT.time[maxinConv]
        log.(convDT.counts[maxinConv-maxin:maxinConv+(length(DTOF.counts)-maxin-1)])
    end

    lb = [0.001, 1]
    ub = [0.5, 80]
    w  = Vector(1:1:length(DTOF.counts)).^1
    β0 = [0.3, 40]
    fit = curve_fit(conv_DT_IRF, IRF.time[1:length(DTOF.time)], log.(DTOF.counts./maximum(DTOF.counts)),w,β0,lower = lb,upper = ub)
end

maxval,maxin = findmax(DTOF.counts)
fit = fitDTOF(DTOF,IRF)
xfit = IRF.time[1:length(DTOF.time)]
yfit = conv_DT_IRF(xfit,fit.param)

scatter(xfit,DTOF.counts./maximum(DTOF.counts), color="black",yaxis =:log,label = "Expt",markersize=3,alpha=0.8)
plot!(xfit,yfit, color="red", linewidth=3, label = "DT-fit",alpha = 0.8)
xlabel!("time (ns)")
ylabel!("Reflectance [counts]")
title!("\\mu_{a} = $(round(fit.param[1],digits=4)) cm^{-1}, \\mu_{s}' = $(round(fit.param[2],digits=3)) cm^{-1}")
plot!(ytickfontsize=10,xtickfontsize=10,margin=5Plots.mm)
#savefig("fit_plot.png")
#savefig("fit_plot.pdf")

max_stars_repo_name: heltonmc/load-asc-data

max_issues_repo_name: heltonmc/load-asc-data

max_forks_repo_name: heltonmc/load-asc-data


text:
"""
Description of the Rastrigin problem, and the setup used to solve it.

Could accept optional parameters (tournament size, pop size, etc).
"""
function rastrigin()
  def = algorithm.genetic() do alg
    alg.population = population.simple() do pop
      pop.size = 100
      pop.species = species.simple() do sp
        sp.fitness = fitness.scalar()
        sp.representation = representation.float_vector(100, -5.12, 5.12)
      end
      pop.breeder = breeder.simple() do br
        br.selection = selection.tournament(2)
        br.mutation = mutation.gaussian() do mut
          mut.rate = 0.01
          mut.mean = 0.0
          mut.std  = 1.0
        end
        br.crossover = crossover.one_point(1.0)
      end
    end
    alg.evaluator = evaluator.simple(; threads = 8) do fs, phenome
      f = reduce(phenome, 1000.0) do f, x
        f += x*x - 10.0 * cos(2.0 * pi * x)
      end
      assign(fs, f)
    end
    alg.termination = [criterion.generations(1000)]
  end
  algorithm.compose!(def)
end


text:
using Documenter, GeoMakie

deploy = length(ARGS) == 1 && ARGS[1] == "deploy"
if deploy && !haskey(ENV, "GITHUB_TOKEN")
    @warn("Not deploying, no GITHUB_TOKEN not found in ENV")
    deploy = false
end

makedocs(;
    modules=[GeoMakie],
    doctest=false,
    format=Documenter.HTML(; prettyurls=deploy, collapselevel=3),
    pages=["GeoMakie.jl" => "index.md"],
    sitename="GeoMakie.jl",
    authors="Makie.jl contributors")

deploy && deploydocs(; repo="github.com/JuliaPlots/GeoMakie.jl", target="build", push_preview=true)


text:
type Sound
    ptr::Ptr{Void}
    buffer

    function Sound(ptr::Ptr{Void})
        s = new(ptr, nothing)
        finalizer(s, destroy)
        s
    end
end

function Sound()
    Sound(ccall((:sfSound_create, libcsfml_audio), Ptr{Void}, ()))
end

function Sound(buffer::SoundBuffer)
    s = Sound()
    set_buffer(s, buffer)
    s
end

function copy(sound::Sound)
    return Sound(ccall((:sfSound_copy, libcsfml_audio), Ptr{Void}, (Ptr{Void},), sound.ptr))
end

function destroy(sound::Sound)
    ccall((:sfSound_destroy, libcsfml_audio), Void, (Ptr{Void},), sound.ptr)
end

function play(sound::Sound)
    ccall((:sfSound_play, libcsfml_audio), Void, (Ptr{Void},), sound.ptr)
end

function pause(sound::Sound)
    ccall((:sfSound_pause, libcsfml_audio), Void, (Ptr{Void},), sound.ptr)
end

function stop(sound::Sound)
    ccall((:sfSound_stop, libcsfml_audio), Void, (Ptr{Void},), sound.ptr)
end

function get_status(sound::Sound)
    ccall((:sfSound_getStatus, libcsfml_audio), Int32, (Ptr{Void},), sound.ptr)
end

function set_buffer(sound::Sound, sound_buffer::SoundBuffer)
    ccall((:sfSound_setBuffer, libcsfml_audio), Void, (Ptr{Void}, Ptr{Void},), sound.ptr, sound_buffer.ptr)
    sound.buffer = sound_buffer
end

function get_buffer(sound::Sound)
    return SoundBuffer(ccall((:sfSound_getBuffer, libcsfml_audio), Void, (Ptr{Void},), sound.ptr))
end

function set_loop(sound::Sound, loop::Bool)
    ccall((:sfSound_setLoop, libcsfml_audio), Void, (Ptr{Void}, Int32,), sound.ptr, loop)
end

function get_loop(sound::Sound)
    return Bool(ccall((:sfSound_getLoop, libcsfml_audio), Int32, (Ptr{Void},), sound.ptr))
end

function set_pitch(sound::Sound, pitch::Real)
    ccall((:sfSound_setPitch, libcsfml_audio), Void, (Ptr{Void}, Cfloat,), sound.ptr, pitch)
end

function get_pitch(sound::Sound)
    return Real(ccall((:sfSound_getPitch, libcsfml_audio), Cfloat, (Ptr{Void},), sound.ptr))
end

function set_volume(sound::Sound, volume::Real)
    ccall((:sfSound_setVolume, libcsfml_audio), Void, (Ptr{Void}, Cfloat,), sound.ptr, volume)
end

function get_volume(sound::Sound)
    return Real(ccall((:sfSound_getVolume, libcsfml_audio), Cfloat, (Ptr{Void},), sound.ptr))
end


text:
import PySCF: pyimport

using TensorOperations: @tensor

#########################################################################
# Additional Concepts
# also see: "Hartree-Fock methods of Quantum Chemistry" by Janos Angyan
#########################################################################
# One-electron properties
function dipole_moment(D, mol)
   # eq. (365) Janos
   # It does the same as PySCF does.
   #ao_dip = @pywith mol.with_common_orig([0,0,0]) begin
   #   mol.intor_symmetric("int1e_r", comp=3)
   #end
   common_orig_bak = mol._env[1:3]
   mol = mol.set_common_orig([0,0,0])
   ao_dip = mol.intor_symmetric("int1e_r", comp=3)
   mol = mol.set_common_orig(common_orig_bak)

   @tensor el_dip[x] := ao_dip[x,i,j] * D[i,j]
   el_dip = real.(el_dip)
   charges  = mol.atom_charges()
   coords   = mol.atom_coords()
   nucl_dip = vec(sum(charges .* coords, dims=1))
   mol_dip  = nucl_dip - el_dip
   _au2debye = 2.5417464157449032
   mol_dip *= _au2debye
end

# Population Analysis/Atomic Charges
@doc raw"""
   `population_analysis(D,S,α,mol)`

Generalized Mulliken, Lowdin, ... population analysis
with `D` the AO density matrix, `S` the overlap matrix,
`α` defines the analysis type, and mol is a PySCF molecule.

Qa = Za - ∑μ∈a [S^α P S^(1-α)]μμ  with 0 < α < 1.

In particular:
* α = 1 => Mulliken doi:10.1063/1.1740588
* α = ½ => Löwdin
"""
function population_analysis(D,S,α,mol;verbose=true)
   s_α, s_mα = S^α, S^(1-α) # Matrix powers (via squaring)
   @assert isapprox(mol.nelectron, tr(s_α*D*s_mα), atol=1e-8) # eq. (363) Janos
   PopAO = diag(s_α * D * s_mα)
   Q = zeros(mol.natm)
   for (μ,record) in enumerate(mol.ao_labels(fmt=nothing))
      a = record[1] + 1 |> n -> convert(Int,n)
      Q[a] -= PopAO[μ] # eqs. (362,364) Janos
   end
   Q += mol.atom_charges()
   if verbose
      for a in eachindex(Q)
         println(" Atom $a has charge $(Q[a])")
      end
   end
   @assert isapprox(sum(Q), mol.charge, atol=1e-8)
   Q
end

function _covalent_bond_order(P,S,mol)
   natm = mol.natm
   nelec = mol.nelectron
   @show nelec
   @show Npairs = sum( (2P * S) * (2P * S) ) # RHF
   ao_atom = first.(mol.ao_labels(fmt=nothing))
   B_ab = zeros(natm,natm)
   for a in 1:natm, b in a:natm
      inds = findall( x -> isequal(x+1,a) || isequal(x+1,b), ao_atom)
      #@show inds
      #Pab = reshape(P[inds,inds],length(inds),length(inds))
      #Sab = reshape(S[inds,inds],length(inds),length(inds))
      Pab = P[inds,inds]
      Sab = S[inds,inds]
      # eq. (375) Janos
      #Bab = sum( (Pab * Sab)*(Pab * Sab)' + (Pab * Sab)*(Pab * Sab)' ) # RHF
      Bab = 2dot( Pab * Sab, Pab * Sab) # RHF
      println("$a  $b  $Bab")
   end
end

# Hartree-Fock Reduced Density Matrices
@doc raw"""
   fock_dirac_denmat(C,S,D)

1RDM in MO basis. For Hartree-Fock it is idempotent, diagonal with 1 or 0 elements.
Therefore it is idempotent. And it is called Fock-Dirac density matrix.

`C` is the matrix of cannonical coefficients, `S` is the overlap matrix
and `D` is the density matrix in the basis of AO.
 """
function fock_dirac_denmat(C,S,D)
   C' * S * D * S * C
end

"""
Nj^2. McWeeny convention (wo. 1/2)
# Arguments
- D in AO basis
- S is ovlp in AO basis
"""
function number_pairs_j(D,S)
   #@tensor scalar(D[i,j] * D[k,l] * S[i,j] * S[k,l])
   #@tensor scalar( (D[i,j] * S[i,j]) * (D[k,l] * S[k,l]))
   @tensor scalar(D[i,j] * S[i,j] * D[k,l] * S[k,l])
end

"""
Nxc. McWeeny convention (wo. 1/2)

# Arguments
- D in AO basis
- S is ovlp in AO basis
"""
function number_pairs_xc(D,S)
   #@tensor 0.5scalar(D[i,j] * D[k,l] * S[i,l] * S[k,j])
   #@tensor 0.5scalar(D[i,j] * D[k,l] * S[l,i] * S[k,j]) # D is symmetric
   #@tensor 0.5scalar(D[i,j] * S[l,i] * D[k,l] * S[k,j]) # D is symmetric
   @tensor 0.5scalar(D[i,j] * S[i,l] * D[k,l] * S[k,j]) # j ↔ l
end

function number_pairs(D,S)
   number_pairs_j(D,S) - number_pairs_xc(D,S)
end

"""
2RD
"""
function _second_reduced_density(rdm1)
   #rdm2D = zeros(nao,nao)
   #for i in 1:nao, j in 1:nao
   rdm2D = similar(rdm1)
   for (i,j) in Tuple.(CartesianIndices(rdm1))
      rdm2D[i,j] = rdm1[i,i]*rdm1[j,j] - rdm1[i,j]*rdm1[j,i]
   end
   rdm2D
end

"""
2RDM
"""
function _second_denmat(D)
   # Choose convention: McWeeny or Löwdin (0.5 factor)
   # Clearly explained in sections 11.15, 13.3 of Janos Angyan notes
   #rdm2 = zeros(nsmo,nsmo,nsmo,nsmo)
   #for i in 1:nao, j in 1:nao, k in 1:nao, l in 1:nao
   #for i in 1:nao,j in 1:nao,k in 1:nao,l in 1:nao
      # 1st RHS is the Coulomb component ∫ = N^2
      # 2nd RHS is the exchange component ∫ = N
      #rdm2[l,k,i,j] =     D[l,i] * D[k,j] - D[l,j] * D[k,i] # eq. (338) Janos
      #lk, ij = index(l,k), index(i,j)
      #rdm2[l,k,i,j] = D[l,i]*D[j,k] - D[l,j]*D[k,i]
      #rdm2[l,j,i,k] = D[l,i]*D[k,j] - D[l,k]*D[j,i]
      #rdm2[i,j,k,l] = 0.5(D[i,j] * D[k,l] - D[i,l] * D[k,j]) # eq. (6) J. of Math. Chem. 28, Nos. 1–3, 2000
   #end
   #rdm2 = zeros(size(D)...,size(D)...)
   #for (i,j,k,l) in Tuple.(CartesianIndices(rdm2))
   #   rdm2[i,j,k,l] = D[i,j] * D[k,l] - D[i,l] * D[k,j]
   #end
   #rdm2
   @tensor rdm2[i,j,k,l] := D[i,j] * D[k,l] - D[i,l] * D[k,j]
end

"""
2RDM Coulomb part (J to distinguish from c for correlation)
"""
function _second_denmat_j(D)
   @tensor rdm2_j[i,j,k,l] := D[i,j] * D[k,l]
end