text:
import Base.show

export CustomParser, Quoted

using Compat

@compat abstract type AbstractToken{T} end
fieldtype{T}(::AbstractToken{T}) = T
fieldtype{T}(::Type{AbstractToken{T}}) = T
fieldtype{T<:AbstractToken}(::Type{T}) = fieldtype(supertype(T))

"""
`tryparsenext{T}(tok::AbstractToken{T}, str, i, till, localopts)`

Parses the string `str` starting at position `i` and ending at or before position `till`. `localopts` is a [LocalOpts](@ref) object which contains contextual options for quoting and NA parsing. (see [LocalOpts](@ref) documentation)

`tryparsenext` returns a tuple `(result, nextpos)` where `result` is of type `Nullable{T}`, null if parsing failed, non-null containing the parsed value if it succeeded. If parsing succeeded, `nextpos` is the position the next token, if any, starts at. If parsing failed, `nextpos` is the position at which the parsing failed.
"""
function tryparsenext end

## options passed down for tokens (specifically NAToken, StringToken)
## inside a Quoted token
"""
    LocalOpts

Options local to the token currently being parsed.
- `endchar`: Till where to parse. (e.g. delimiter or quote ending character)
- `spacedelim`: Treat spaces as delimiters
- `quotechar`: the quote character
- `escapechar`: char that escapes the quote
- `includequotes`: whether to include quotes while parsing
- `includenewlines`: whether to include newlines while parsing
"""
immutable LocalOpts
    endchar::Char         # End parsing at this char
    spacedelim::Bool
    quotechar::Char       # Quote char
    escapechar::Char      # Escape char
    includequotes::Bool   # Whether to include quotes in string parsing
    includenewlines::Bool # Whether to include newlines in string parsing
end

const default_opts = LocalOpts(',', false, '"', '\\', false, false)
# helper function for easy testing:
@inline function tryparsenext(tok::AbstractToken, str, opts::LocalOpts=default_opts)
    tryparsenext(tok, str, start(str), endof(str), opts)
end

# fallback for tryparsenext methods which don't care about local opts
@inline function tryparsenext(tok::AbstractToken, str, i, len, locopts)
    tryparsenext(tok, str, i, len)
end

immutable WrapLocalOpts{T, X<:AbstractToken} <: AbstractToken{T}
    opts::LocalOpts
    inner::X
end

WrapLocalOpts(opts, inner) = WrapLocalOpts{fieldtype(inner), typeof(inner)}(opts, inner)

@inline function tryparsenext(tok::WrapLocalOpts, str, i, len, opts::LocalOpts=default_opts)
    tryparsenext(tok.inner, str, i, len, tok.opts)
end

# needed for promoting guessses
immutable Unknown <: AbstractToken{Union{}} end
fromtype(::Type{Union{}}) = Unknown()
const nullableNA = Nullable{DataValue{Union{}}}(DataValue{Union{}}())
function tryparsenext(::Unknown, str, i, len, opts)
    nullableNA, i
end
show(io::IO, ::Unknown) = print(io, "<unknown>")
immutable CustomParser{T, F} <: AbstractToken{T}
    f::Function
end

"""
    CustomParser(f, T)

Provide a custom parsing mechanism.

# Arguments:

- `f`: the parser function
- `T`: The type of the parsed value

The parser function must take the following arguments:
- `str`: the entire string being parsed
- `pos`: the position in the string at which to start parsing
- `len`: the length of the string the maximum position where to parse till
- `opts`: a [LocalOpts](@ref) object with options local to the current field.

The parser function must return a tuple of two values:

- `result`: A `Nullable{T}`. Set to null if parsing must fail, containing the value otherwise.
- `nextpos`: If parsing succeeded this must be the next position after parsing finished, if it failed this must be the position at which parsing failed.
"""
CustomParser(f, T) = CustomParser{T,typeof(f)}(f)

show{T}(io::IO, c::CustomParser{T}) = print(io, "{{custom:$T}}")

@inline function tryparsenext(c::CustomParser, str, i, len, opts)
    c.f(str, i, len, opts)
end

# Numberic parsing
"""
parse numbers of type T
"""
immutable Numeric{T} <: AbstractToken{T}
    decimal::Char
    thousands::Char
end
show{T}(io::IO, c::Numeric{T}) = print(io, "<$T>")

Numeric{T}(::Type{T}, decimal='.', thousands=',') = Numeric{T}(decimal, thousands)
fromtype{N<:Number}(::Type{N}) = Numeric(N)

### Unsigned integers

function tryparsenext{T<:Signed}(::Numeric{T}, str, i, len)
    R = Nullable{T}
    @chk2 sign, i = tryparsenext_sign(str, i, len)
    @chk2 x, i = tryparsenext_base10(T, str, i, len)

    @label done
    return R(sign*x), i

    @label error
    return R(), i
end

@inline function tryparsenext{T<:Unsigned}(::Numeric{T}, str, i, len)
    tryparsenext_base10(T,str, i, len)
end

@inline function tryparsenext{F<:AbstractFloat}(::Numeric{F}, str, i, len)
    R = Nullable{F}
    f = 0.0
    @chk2 sign, i = tryparsenext_sign(str, i, len)
    x=0

    i > len && @goto error
    c, ii = next(str, i)
    if c == '.'
        i=ii
        @goto dec
    end
    @chk2 x, i = tryparsenext_base10(Int, str, i, len)
    i > len && @goto done
    @inbounds c, ii = next(str, i)

    c != '.' && @goto parse_e
    @label dec
    @chk2 y, i = tryparsenext_base10(Int, str, ii, len) done
    f = y / 10.0^(i-ii)

    i > len && @goto done
    c, ii = next(str, i)

    @label parse_e
    if c == 'e' || c == 'E'
        @chk2 exp, i = tryparsenext(Numeric(Int), str, ii, len)
        return R(sign*(x+f) * 10.0^exp), i
    end

    @label done
    return R(sign*(x+f)), i

    @label error
    return R(), i
end

immutable Percentage <: AbstractToken{Float64}
end

const floatparser = Numeric(Float64)
function tryparsenext(::Percentage, str, i, len, opts)
    num, ii = tryparsenext(floatparser, str, i, len, opts)
    if isnull(num)
        return num, ii
    else
        # parse away the % char
        ii = eatwhitespaces(str, ii, len)
        c, k = next(str, ii)
        if c != '%'
            return Nullable{Float64}(), ii # failed to parse %
        else
            return Nullable{Float64}(num.value / 100.0), k # the point after %
        end
    end
end

"""
Parses string to the AbstractString type `T`. If `T` is `StrRange` returns a
`StrRange` with start position (`offset`) and `length` of the substring.
It is used internally by `csvparse` for avoiding allocating strings.
"""
immutable StringToken{T} <: AbstractToken{T}
end

function StringToken{T}(t::Type{T})
    StringToken{T}()
end
show(io::IO, c::StringToken) = print(io, "<string>")

fromtype{S<:AbstractString}(::Type{S}) = StringToken(S)

function tryparsenext{T}(s::StringToken{T}, str, i, len, opts)
    R = Nullable{T}
    p = ' '
    i0 = i
    if opts.includequotes && i <= len
        c, ii = next(str, i)
        if c == opts.quotechar
            i = ii # advance counter so that
                   # the while loop doesn't react to opening quote
        end
    end

    while i <= len
        c, ii = next(str, i)
        if opts.spacedelim && (c == ' ' || c == '\t')
            break
        elseif !opts.spacedelim && c == opts.endchar
            if opts.endchar == opts.quotechar
                # this means we're inside a quoted string
                if opts.quotechar == opts.escapechar
                    # sometimes the quotechar is the escapechar
                    # in that case we need to see the next char
                    if ii > len
                        if opts.includequotes
                            i=ii
                        end
                        break
                    end
                    nxt, j = next(str, ii)
                    if nxt == opts.quotechar
                        # the current character is escaping the
                        # next one
                        i = j # skip next char as well
                        p = nxt
                        continue
                    end
                elseif p == opts.escapechar
                    # previous char escaped this one
                    i = ii
                    p = c
                    continue
                end
            end
            if opts.includequotes
                i = ii
            end
            break
        elseif (!opts.includenewlines && isnewline(c))
            break
        end
        i = ii
        p = c
    end

    return R(_substring(T, str, i0, i-1)), i
end

@inline function _substring(::Type{String}, str, i, j)
    str[i:j]
end

if VERSION <= v"0.6.0-dev"
    # from lib/Str.jl
    @inline function _substring(::Type{Str}, str, i, j)
        Str(pointer(Vector{UInt8}(str))+(i-1), j-i+1)
    end
end

@inline function _substring{T<:SubString}(::Type{T}, str, i, j)
    T(str, i, j)
end

fromtype(::Type{StrRange}) = StringToken(StrRange)

@inline function alloc_string(str, r::StrRange)
    unsafe_string(_pointer(str, 1+r.offset), r.length)
end

@inline function _substring(::Type{StrRange}, str, i, j)
    StrRange(i-1, j-i+1)
end

export Quoted

immutable Quoted{T, S<:AbstractToken} <: AbstractToken{T}
    inner::S
    required::Bool
    stripwhitespaces::Bool
    includequotes::Bool
    includenewlines::Bool
    quotechar::Nullable{Char}
    escapechar::Nullable{Char}
end

function show(io::IO, q::Quoted)
    c = quotechar(q, default_opts)
    print(io, "$c")
    show(io, q.inner)
    print(io, "$c")
end

"""
`Quoted(inner::AbstractToken; <kwargs>...)`

# Arguments:
- `inner`: The token inside quotes to parse
- `required`: are quotes required for parsing to succeed? defaults to `false`
- `includequotes`: include the quotes in the output. Defaults to `false`
- `includenewlines`: include newlines that appear within quotes. Defaults to `true`
- `quotechar`: character to use to quote (default decided by `LocalOpts`)
- `escapechar`: character that escapes the quote char (default set by `LocalOpts`)
"""
function Quoted{S<:AbstractToken}(inner::S;
    required=false,
    stripwhitespaces=fieldtype(S)<:Number,
    includequotes=false,
    includenewlines=true,
    quotechar=Nullable{Char}(),   # This is to allow file-wide config
    escapechar=Nullable{Char}())

    T = fieldtype(S)
    Quoted{T,S}(inner, required, stripwhitespaces, includequotes,
                includenewlines, quotechar, escapechar)
end

@inline quotechar(q::Quoted, opts) = get(q.quotechar, opts.quotechar)
@inline escapechar(q::Quoted, opts) = get(q.escapechar, opts.escapechar)

Quoted(t::Type; kwargs...) = Quoted(fromtype(t); kwargs...)

function tryparsenext{T}(q::Quoted{T}, str, i, len, opts)
    if i > len
        q.required && @goto error
        # check to see if inner thing is ok with an empty field
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts) error
        @goto done
    end
    c, ii = next(str, i)
    quotestarted = false
    if quotechar(q, opts) == c
        quotestarted = true
        if !q.includequotes
            i = ii
        end

        if q.stripwhitespaces
            i = eatwhitespaces(str, i)
        end
    else
        q.required && @goto error
    end

    if quotestarted
        qopts = LocalOpts(quotechar(q, opts), false, quotechar(q, opts), escapechar(q, opts),
                         q.includequotes, q.includenewlines)
        @chk2 x, i = tryparsenext(q.inner, str, i, len, qopts)
    else
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts)
    end

    if i > len
        if quotestarted && !q.includequotes
            @goto error
        end
        @goto done
    end

    if q.stripwhitespaces
        i = eatwhitespaces(str, i)
    end
    c, ii = next(str, i)

    if quotestarted && !q.includequotes
        c != quotechar(q, opts) && @goto error
        i = ii
    end

    @label done
    return Nullable{T}(x), i

    @label error
    return Nullable{T}(), i
end

## Date and Time
immutable DateTimeToken{T,S<:DateFormat} <: AbstractToken{T}
    format::S
end

"""
    DateTimeToken(T, fmt::DateFormat)

Parse a date time string of format `fmt` into type `T` which is
either `Date`, `Time` or `DateTime`.
"""
DateTimeToken{S<:DateFormat}(T::Type, df::S) = DateTimeToken{T, S}(df)
DateTimeToken{S<:DateFormat}(df::S) = DateTimeToken{DateTime, S}(df)
fromtype(df::DateFormat) = DateTimeToken(DateTime, df)
fromtype(::Type{DateTime}) = DateTimeToken(DateTime, ISODateTimeFormat)
fromtype(::Type{Date}) = DateTimeToken(Date, ISODateFormat)

function fromtype(nd::Union{Nullable{DateFormat}, DataValue{DateFormat}})
    if !isnull(nd)
        NAToken(DateTimeToken(DateTime, get(nd)))
    else
        fromtype(DataValue{DateTime})
    end
end

function tryparsenext{T}(dt::DateTimeToken{T}, str, i, len, opts)
    R = Nullable{T}
    nt, i = tryparsenext_internal(T, str, i, len, dt.format, opts.endchar)
    if isnull(nt)
        return R(), i
    else
        return R(T(nt.value...)), i
    end
end

### Nullable

const nastrings_upcase = ["NA", "NULL", "N/A","#N/A", "#N/A N/A", "#NA",
                          "-1.#IND", "-1.#QNAN", "-NaN", "-nan",
                          "1.#IND", "1.#QNAN", "N/A", "NA", "NaN", "nan"]

const NA_STRINGS = sort!(vcat(nastrings_upcase, map(lowercase, nastrings_upcase)))

immutable NAToken{T, S<:AbstractToken} <: AbstractToken{T}
    inner::S
    emptyisna::Bool
    endchar::Nullable{Char}
    nastrings::Vector{String}
end

"""
`NAToken(inner::AbstractToken; options...)`

Parses a Nullable item.

# Arguments
- `inner`: the token to parse if non-null.
- `emptyisna`: should an empty item be considered NA? defaults to true
- `nastrings`: strings that are to be considered NA. Defaults to `$NA_STRINGS`
"""
function NAToken{S}(
    inner::S,
  ; emptyisna=true
  , endchar=Nullable{Char}()
  , nastrings=NA_STRINGS)

    T = fieldtype(inner)
    NAToken{DataValue{T}, S}(inner, emptyisna, endchar, nastrings)
end

function show(io::IO, na::NAToken)
    show(io, na.inner)
    print(io, "?")
end

endchar(na::NAToken, opts) = get(na.endchar, opts.endchar)

function tryparsenext{T}(na::NAToken{T}, str, i, len, opts)
    R = Nullable{T}
    i = eatwhitespaces(str, i)
    if i > len
        if na.emptyisna
            @goto null
        else
            @goto error
        end
    end

    c, ii=next(str,i)
    if (c == endchar(na, opts) || isnewline(c)) && na.emptyisna
       @goto null
    end

    if isa(na.inner, Unknown)
        @goto maybe_null
    end
    @chk2 x,ii = tryparsenext(na.inner, str, i, len, opts) maybe_null

    @label done
    return R(T(x)), ii

    @label maybe_null
    naopts = LocalOpts(endchar(na,opts), opts.spacedelim, opts.quotechar,
                       opts.escapechar, false, opts.includenewlines)
    @chk2 nastr, ii = tryparsenext(StringToken(String), str, i, len, naopts)
    if !isempty(searchsorted(na.nastrings, nastr))
        i=ii
        i = eatwhitespaces(str, i)
        @goto null
    end
    return R(), i

    @label null
    return R(T()), i

    @label error
    return R(), i
end

fromtype{N<:Nullable}(::Type{N}) = NAToken(fromtype(eltype(N)))
fromtype{N<:DataValue}(::Type{N}) = NAToken(fromtype(eltype(N)))

### Field parsing

@compat abstract type AbstractField{T} <: AbstractToken{T} end # A rocord is a collection of abstract fields

immutable Field{T,S<:AbstractToken} <: AbstractField{T}
    inner::S
    ignore_init_whitespace::Bool
    ignore_end_whitespace::Bool
    eoldelim::Bool
end

function Field{S}(inner::S; ignore_init_whitespace=true, ignore_end_whitespace=true, eoldelim=false)
    T = fieldtype(inner)
    Field{T,S}(inner, ignore_init_whitespace, ignore_end_whitespace, eoldelim)
end

function Field(f::Field; inner=f.inner, ignore_init_whitespace=f.ignore_init_whitespace,
                  ignore_end_whitespace=f.ignore_end_whitespace,
                  eoldelim=f.eoldelim)
    T = fieldtype(inner)
    Field{T,typeof(inner)}(inner, ignore_init_whitespace,
                           ignore_end_whitespace, eoldelim)
end

function swapinner(f::Field, inner::AbstractToken;
        ignore_init_whitespace= f.ignore_end_whitespace
      , ignore_end_whitespace=f.ignore_end_whitespace
      , eoldelim=f.eoldelim
  )
    Field(inner;
        ignore_init_whitespace=ignore_end_whitespace
      , ignore_end_whitespace=ignore_end_whitespace
      , eoldelim=eoldelim
     )

end

function tryparsenext{T}(f::Field{T}, str, i, len, opts)
    R = Nullable{T}
    i > len && @goto error
    if f.ignore_init_whitespace
        while i <= len
            @inbounds c, ii = next(str, i)
            !isspace(c) && break
            i = ii
        end
    end
    @chk2 res, i = tryparsenext(f.inner, str, i, len, opts)

    if f.ignore_end_whitespace
        i0 = i
        while i <= len
            @inbounds c, ii = next(str, i)
            !opts.spacedelim && opts.endchar == '\t' && c == '\t' && (i =ii; @goto done)
            !isspace(c) && c != '\t' && break
            i = ii
        end

        opts.spacedelim && i > i0 && @goto done
    end
    # todo don't ignore whitespace AND spacedelim

    if i > len
        if f.eoldelim
            @goto done
        else
            @goto error
        end
    end

    @inbounds c, ii = next(str, i)
    opts.spacedelim && (isspace(c) || c == '\t') && (i=ii; @goto done)
    !opts.spacedelim && opts.endchar == c && (i=ii; @goto done)

    if f.eoldelim
        if c == '\r'
            i=ii
            if i <= len
                @inbounds c, ii = next(str, i)
                if c == '\n'
                    i=ii
                end
            end
            @goto done
        elseif c == '\n'
            i=ii
            if i <= len
                @inbounds c, ii = next(str, i)
                if c == '\r'
                    i=ii
                end
            end
            @goto done
        end
    end

    @label error
    return R(), i

    @label done
    if R <: Nullable{DataValue{Union{}}}
        # optimization to remove allocation
        return nullableNA, i
    end
    return R(res), i
end


text:
# http://web.mit.edu/puzzle/www/2013/coinheist.com/rubik/substance_abuse/index.html
@testset "substance abuse" begin
    f = best_feature(wordlist("""
            KBNSCSICRBNA
            PSICAFCRLIO
            CRLISCNCA
            LIKBBECALI
            NCABPARNSCBNA
            CABNEHBKNSCHESCP
            BCSIOCLNCBNA
            NEHBNATIHEBCASLIO
            BHESIPOSIPNCSICA
            SINASCNCLI
            NNALISCSICLI
            NAARCASIOSIMGSIOCR
            NEHNOSCAL
            NCRCRSICBN
            NASIOHCKHCR
            CAFLI
            PBSCNAARBECALICKLI
            NANPHENBNABC
            SCARFCRSICA
            HELIOSISCSICBC
            PHEBCASINAFBEBC
            LINAHESCNHEF
            SNCBCACABC
        """), 5)
    @test f.description == "can be completely broken down into chemical element symbols"
end

ext: jl


text:
# Autogenerated wrapper script for MPIR_jll for x86_64-w64-mingw32-cxx11
export libmpir

JLLWrappers.@generate_wrapper_header("MPIR")
JLLWrappers.@declare_library_product(libmpir, "libmpir-23.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmpir,
        "bin\\libmpir-23.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    StableNumber(base::Number, exponent)

Stores results from munipulations of `StableArray`.
"""
struct StableNumber{E, A <: Number}
    base::A
    exponent::E
end

base(x::StableNumber) = x.base
exponent(x::StableNumber) = x.exponent

function _plus_ordered(x1::StableNumber, x2::StableNumber)
    xbase = base(x1) * exp(exponent(x1) - exponent(x2)) + base(x2)
    StableNumber(xbase, exponent(x2))
end

function Base.:(+)(x1::StableNumber, x2::StableNumber)
    indx = argmax((x1.exponent, x2.exponent))
    if indx == 2
        _plus_ordered(x1, x2)
    else
        _plus_ordered(x2, x1)
    end
end

function Base.:(+)(x1::StableNumber, x2::Number)
    x1 + StableNumber(x2, 1)
end

function Base.:(+)(x1::Number, x2::StableNumber)
    StableNumber(x1, 1) + x2
end

function Base.:(-)(x::StableNumber)
    StableNumber(-base(x), exponent(x))
end

function Base.:(-)(x1::StableNumber, x2::StableNumber)
    x1 + (-x2)
end

function Base.:(-)(x1::StableNumber, x2)
    x1 + (-x2)
end

function Base.:(-)(x1, x2::StableNumber)
    x1 + (-x2)
end

function Base.:(*)(x1::StableNumber, x2::StableNumber)
    StableNumber(base(x1) * base(x2), exponent(x1) + exponent(x2))
end

function Base.:(/)(x1::StableNumber, x2::StableNumber)
    StableNumber(base(x1) / base(x2), exponent(x1) - exponent(x2))
end

function Base.:(/)(x1::StableNumber, x2::Number)
    StableNumber(base(x1) * sign(x2), exponent(x1) - log(abs(x2)))
end

function Base.:(<)(x1::StableNumber, x2::StableNumber)
    xd = x1 - x2
    base(xd) < 0
end

function Base.:(<)(x1, x2::StableNumber)
    xd = x1 - x2
    base(xd) < 0
end

function Base.:(<)(x1::StableNumber, x2)
    xd = x1 - x2
    base(xd) < 0
end

function Base.log(x::StableNumber)
    log(base(x)) + exponent(x)
end

function Base.abs(x::StableNumber)
    StableNumber(abs(base(x)), exponent(x))
end

function Base.isnan(x::StableNumber)
    isnan(base(x)) || isnan(exponent(x))
end

function Base.show(io::IO, ::MIME"text/plain", x::StableNumber)
    print(io, "Stabled($(base(x))×exp($(exponent(x))))")
end

function Base.sign(x::StableNumber)
    sign(base(x))
end

function logview(x::StableNumber)
    StableNumber(sign(x), exponent(x)+log(abs(base(x))))
end


text:
using SQLiteGraph
using Test
using JSON3

#-----------------------------------------------------------------------------# setup
db = DB()

#-----------------------------------------------------------------------------# Nodes
@testset "Nodes" begin
    @testset "Round Trips" begin
        for n in [
                Node(1),
                Node(2, "lab"),
                Node(3, "lab1", "lab2"),
                Node(4, "lab"; x=1),
                Node(5, "lab1", "lab2"; x=1, y=2)
            ]
            insert!(db, n)
            @test db[n.id] == n
            @test_throws Exception insert!(db, n)
        end
    end
    @testset "replace!" begin
        replace!(db, Node(1, "lab"))
        @test db[1].labels == ["lab"]
    end
    @testset "simple query" begin
        q = db[:]
        for (i,n) in enumerate(q)
            @test n.id == i
        end
        @test length(collect(db[:])) == 5
    end
end

#-----------------------------------------------------------------------------# Edges
@testset "Edges" begin
    @testset "Round Trips" begin
        for e in [
                Edge(1,2,"type"),
                Edge(1,3,"type"; x=1),
                Edge(1,4,"type 2"; x=1,y=2,z=3)
            ]
            insert!(db, e)
            @test db[e.source, e.target, e.type] == e
            @test_throws Exception insert!(db, e)
        end
    end
    @testset "replace!" begin
        replace!(db, Edge(1,2,"type"; x=1))
        @test db[1,2,"type"].props.x == 1
    end
    @testset "simple query" begin
        @test db[1,2,"type"] == Edge(1,2,"type"; x=1)
        @test length(collect(db[1,2,:])) == 1
        @test length(collect(db[1,:,"type"])) == 2
        @test length(collect(db[:,2,"type"])) == 1
        @test length(collect(db[:,:,"type"])) == 2
        @test length(collect(db[:, 4, :])) == 1
        @test length(collect(db[:,:,"type"])) == 2
        @test length(collect(db[:,:,:])) == 3
    end
    @testset "Adj Matrix" begin
        m = SQLiteGraph.adjacency_matrix(db, "type")
        @test sum(m) == 2
        @test m[1, 2]
        @test m[1, 3]

        m2 = SQLiteGraph.adjacency_matrix(db, "type 2")
        @test sum(m2) == 1
        @test m2[1, 4]
    end
end


text:
using StatsBase

import StatsBase.zscore!
zscore!(X,dim::Tuple{Int64}) = begin
    d = dim[1]
    for i in indices(X,d)
        t = slicedimview(X, 2, i)
        t .= (t.-mean(t))./std(t)
    end
    X
end

zscore(X,dim::Tuple{Int64}) = begin
    Y = Float64.(X)
    zscore!(Y,dim)
    return Y
end

function zscore_base(X::Vector{T},inds) where T<:Real
    (X .- mean(X[inds]))/std(X[inds])
end

function zscore_base(X::Array{T,N},inds,dim::Tuple{Int64}) where {T<:Real,N}
    Y = Float64.(X)
    zscore_base!(Y,inds,dim)
    return Y
end

function zscore_base!(X::Array{T,N},inds,dim::Tuple{Int64}) where {T<:Real,N}
    d = dim[1]
    for i in axes(X,d)
        t = slicedimview(X, 2, i)
        t .= (t.-mean(t[inds]))./std(t[inds])
    end
    X
end

function _bootstrap(f::Function, args...; nreps=100)
    μ   = f(args...)
    res = typeof(μ)[]
    lena = length(args[1])
    @assert all(lena .== length.(args))
    for _ in 1:nreps
        samples = sample(1:lena, lena; replace=true)
        args_bootstrap = [a[samples] for a in args]
        push!(res, f(args_bootstrap...))
    end
    return μ, res
end

"""
     bootstrap_confint(f::Function, args...; nreps=100, perc=95)
Bootstrap a function over multiple arguements *with the same length*.
Returns a three-tuple:
- mean (result of function on complete args)
- hi (\$perc'th percentile of results of function applied on re-sampled arguments)
- lo (1-\$perc'th percetile of results of function applied on re-sampled arguments)
"""
function bootstrap_confint(f::Function, args...; nreps=100, perc=95)
    μ,res = _bootstrap(f, args..., nreps=nreps)
    return μ, percentile(res,perc), percentile(res,100-perc)
end


text: """Create sample weather scenario."""


text:
# Based upon the paper for parallel implicit sindy
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

function cart_pole(u, p, t)
    du = similar(u)
    F = -0.2 + 0.5*sin(6*t) # the input
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -(19.62*sin(u[1])+sin(u[1])*cos(u[1])*u[3]^2+F*cos(u[1]))/(2-cos(u[1])^2)
    du[4] = -(sin(u[1])*u[3]^2 + 9.81*sin(u[1])*cos(u[1])+F)/(2-cos(u[1])^2)
    return du
end

u0 = [0.3; 0; 1.0; 0]
tspan = (0.0, 16.0)
dt = 0.001
cart_pole_prob = ODEProblem(cart_pole, u0, tspan)
solution = solve(cart_pole_prob, Tsit5(), saveat = dt)

plot(solution)

# Create the differential data
X = solution[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = cart_pole(xi, [], solution.t[i])
end

plot(DX')

@variables u[1:4] t
polys = Operation[]
# Lots of basis functions -> sindy pi can handle more than ADM()
for i ∈ 0:4
    if i == 0
        push!(polys, u[1]^0)
    else
        if i < 2
            push!(polys, u.^i...)
        else
            push!(polys, u[3:4].^i...)
        end

    end
end
push!(polys, sin.(u[1])...)
push!(polys, cos.(u[1]))
push!(polys, sin.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*u[3:4].^2...)
push!(polys, cos.(u[1]).^2...)
push!(polys, sin.(u[1]).*cos.(u[1])...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4].^2...)
push!(polys, -0.2+0.5*sin(6*t))
basis= Basis(polys, u, iv = t)

# Simply use any optimizer you would use for sindy
λ = exp10.(-4:0.1:-1)
g(x) = norm([1e-3; 10.0] .* x, 2)
Ψ = ISINDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
println(Ψ)
print_equations(Ψ, show_parameter = true)

# Transform into ODE System
sys = ODESystem(Ψ)
dudt = ODEFunction(sys)
ps = parameters(Ψ)

# Simulate
estimator = ODEProblem(dudt, u0, tspan, ps)
sol_ = solve(estimator, Tsit5(), saveat = dt)

# Yeah! We got it right
plot(solution.t[:], solution[:,:]', color = :red, label = nothing)
plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation")

plot(solution.t, abs.(solution-sol_)')
norm(solution[:,:]-sol_[:,:], 2) # ≈ 0.018


text:
######################################################################
# edwardspoints.jl: Projective and extended points for Edwards curves
######################################################################

export EdwardsPoint, ExtendedEdwardsPoint, neutral, isneutral, extended_neutral

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

"""
Concrete type for classical projective points on Edwards curves.
"""
mutable struct EdwardsPoint{T<:Nemo.RingElem} <: ProjectivePoint{T}
    X::T
	Y::T
    Z::T
    curve::Edwards{T}
end

"""
Concrete type for extended projective points on Edwards curves.
"""
mutable struct ExtendedEdwardsPoint{T<:Nemo.RingElem} <: ProjectivePoint{T}
    X::T
	Y::T
    Z::T
	U::T
    curve::Edwards{T}
end

base_ring(P::EdwardsPoint) = parent(P.X)

base_curve(P::EdwardsPoint) = P.curve

base_ring(P::ExtendedEdwardsPoint) = parent(P.X)

base_curve(P::ExtendedEdwardsPoint) = P.curve

"""
Get a description of a classical projective point on Edwards curves.
"""
function show(io::IO, P::EdwardsPoint)
	print(io, "Point ($(P.X):&(P.Y):$(P.Z)) on ")
	show(io, P.curve)
end

"""
Get a description of an extended projective point on Edwards curves.
"""
function show(io::IO, P::ExtendedEdwardsPoint)
	print(io, "Extended point ($(P.X):&(P.Y):$(P.Z):$P.U) on ")
	show(io, P.curve)
end

######################################################################
# Conversions
######################################################################

EdwardsPoint(x, y, E::Edwards) = EdwardsPoint(x, y, one(parent(X)), E)

ExtendedEdwardsPoint(P::EdwardsPoint) = ExtendedEdwardsPoint(P.X * P.Z, P.Y * P.Z, P.Z^2, P.X * P.Y, P.curve)

EdwardsPoint(P::ExtendedEdwardsPoint) = EdwardsPoint(P.X, P.Y, P.Z, P.curve)

######################################################################
# Normalization, equality and validity
######################################################################

normalized(P::EdwardsPoint) = EdwardsPoint(P.X // P.Z, P.Y // P.Z, one(base_ring(P)), P.curve)

function isvalid(P::EdwardsPoint)
	Eq = projective_equation(P.curve)
	return (P.Z != 0) && Eq(P.X, P.Y, P.Z) == 0
end

function samefields(P::EdwardsPoint, Q::EdwardsPoint)
	return (P.X == Q.X) && (P.Y == Q.Y) && (P.Z == Q.Z) && (P.curve == Q.curve)
end

function ==(P::EdwardsPoint, Q::EdwardsPoint)
	c = P.Z // Q.Z
	return (P.X = c * Q.X) && (P.Y = c * Q.Y) && (P.curve == Q.curve)
end

function isvalid(P::ExtendedEdwardsPoint)
	Q = EdwardsPoint(P)
	return isvalid(Q) && (P.U * P.Z == P.X * P.Y)
end

function samefields(P::ExtendedEdwardsPoint, Q::ExtendedEdwardsPoint)
	return (P.X == Q.X) && (P.Y == Q.Y) && (P.Z == Q.Z) && (P.U == Q.U) && (P.curve == Q.curve)
end

function ==(P::ExtendedEdwardsPoint, Q::ExtendedEdwardsPoint)
	c = P.Z // Q.Z
	return (P.X = c * Q.X) && (P.Y = c * Q.Y) && (P.U == Q.U) && (P.curve == Q.curve)
end

######################################################################
# Special points
######################################################################

function neutral(E::Edwards)
	K = base_ring(E)
	return EdwardsPoint(zero(K), one(K), one(K))
end

isneutral(P::EdwardsPoint) = P == neutral(P.curve)

extended_neutral(E::Edwards) = ExtendedEdwardsPoint(neutral(E))

isneutral(E::ExtendedEdwardsPoint) = isneutral(EdwardsPoint(P))

######################################################################
# Arithmetic of ordinary projective points
######################################################################

complete_formulas(E::Edwards) = issquare(E.a) && !issquare(E.d)

#For now addition formulae *assume* that the formulas are complete.

function +(P::EdwardsPoint, Q::EdwardsPoint)
	Ed = P.curve
	A = P.Z * Q.Z
	B = A^2
	C = P.X * Q.X
	D = P.Y * Q.Y
	E = Ed.d * C * D
	F = B - E
	G = B + E
	resx = A * F * ((P.X + Q.X) * (P.Y + Q.Y) - C - D)
	resy = A * G * (D - Ed.a * C)
	resz = F * G
	return EdwardsPoint(resx, resy, resz, Ed)
end

function double(P::EdwardsPoint)
	Ed = P.curve
	B = (P.X + P.Y)^2
	C = P.X^2
	D = P.Y^2
	E = Ed.a * C
	F = E + D
	H = P.Z^2
	J = F - 2*H
	resx = (B - C - D) * J
	resy = F * (E - D)
	resz = F * J
	return EdwardsPoint(resx, resy, resz, Ed)
end

function *(k::Nemo.Integer, P::EdwardsPoint)
	E = P.curve
	P = normalized(P)
	if k == 0
		return neutral(E)
	elseif k<0
		return (-k) * (-P)
	else
		if isneutral(P)
			return neutral(E)
		else
			return ladder(k, P)
		end
	end
end

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

function ladder(m::Integer, P::EdwardsPoint)
	p0 = P
	for b in bin(m)[2:end]
        if (b == '0')
        	p0 = double(p0)
        else
        	p0 = P + double(p0)
        end
    end
	return p0
end

 -(P::EdwardsPoint) = EdwardsPoint( - P.X, P.Y, P.Z, P.curve)

 -(P::EdwardsPoint, Q::EdwardsPoint) = P + (-Q)

######################################################################
# Arithmetic of extended projective points
######################################################################

function +(P::ExtendedEdwardsPoint, Q::ExtendedEdwardsPoint)
	Ed = P.curve
	A = P.X * Q.X
	B = P.Y * Q.Y
	C = Ed.d * P.U * Q.U
	D = P.Z * Q.Z
	E = (P.X + Q.X) * (P.Y + Q.Y) - A - B
	F = D - C
	G = D + C
	H = B - Ed.a * A
	resx = E * F
	resy = G * H
	resz = E * H
	resu = F * G
	return ExtendedEdwardsPoint(resx, resy, resz, resu, Ed)
end

function double(P::ExtendedEdwardsPoint)
	Ed = P.curve
	A = P.X^2
	B = P.Y^2
	C = 2 * P.Z^2
	D = Ed.a * A
	E = (P.X + P.Y)^2 - A - B
	G = D + B
	F = G - C
	H = D - B
	resx = E * F
	resy = G * H
	resz = F * G
	resu = E * H
	return ExtendedEdwardsPoint(resx, resy, resz, resu, Ed)
end

function *(k::Integer, P::ExtendedEdwardsPoint)
	E = P.curve
	P = normalized(P)
	if k == 0
		return neutral(E)
	elseif k<0
		return (-k) * (-P)
	else
		if isneutral(P)
			return neutral(E)
		else
			return ladder(k, P)
		end
	end
end

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

function ladder(m::Nemo.Integer, P::ExtendedEdwardsPoint)
	p0 = P
	for b in bin(m)[2:end]
        if (b == '0')
        	p0 = double(p0)
        else
        	p0 = P + double(p0)
        end
    end
	return p0
end

 -(P::ExtendedEdwardsPoint) = ExtendedEdwardsPoint( - P.X, P.Y, P.Z, -P.U, P.curve)

 -(P::ExtendedEdwardsPoint, Q::ExtendedEdwardsPoint) = P + (-Q)


text:
function plot_comp_static_firm(x, y; xlabel = "", ylabel = "", Filename::String, Folder::String)
    y_LU = y[1]
    y_BU = y[2]
    y_LS = y[3]
    y_BS = y[4]

    plot(x, y_LU ./ y_LU[1], xlabel = xlabel, ylabel = ylabel, label = "Small, unsophisticated", color = :black, ls = :dashdot)
    plot!(x, y_BU ./ y_BU[1], color = :blue, ls = :dot, label = "Big, unsophisticated")
    plot!(x, y_LS ./ y_LS[1], color = :green, ls = :dash, label = "Small, sophisticated")
    plot!(x, y_BS ./ y_BS[1], color = :red, label = "Big, sophisticated")

    savefig("../Graphs/v3/" *Folder* "/" * Filename * "_normalized.pdf")
    savefig("../Graphs/v3/" *Folder* "/" * Filename * "_normalized.png")

    plot(x, y_LU, color = :black, ls = :dashdot, xlabel = xlabel, ylabel = ylabel, label = "Small, unsophisticated")
    plot!(x, y_BU, color = :blue, ls = :dot, label = "Big, unsophisticated")
    plot!(x, y_LS, color = :green, ls = :dash, label = "Small, sophisticated")
    plot!(x, y_BS, color = :red, label = "Big, sophisticated")

    savefig("../Graphs/v3/" *Folder* "/" * Filename * ".pdf")
    savefig("../Graphs/v3/" *Folder* "/" * Filename * ".png")
end


text:
@info "Start Testing"
using PotentialMap
using SteadyStateDiffEq
using FindSteadyStates
using PyPlot
using Test

include("SampleDE.jl")

plt= PyPlot
sp = SampleDE

@info "Start Testing PotentialMap.LeastAction"
@testset "PotentialMap.jl" begin

    @testset "Generate diffusion matrix" begin
        include("diffusionM.jl")
    end

    @testset "Wrapper of FindSteadyStates" begin
        include("find_fixedpoints.jl")
    end

    @testset "gradient calculation compared to the reference" begin
        include("gradient.jl")
    end

    @info "Start Measuring LeastAction"
    @testset "LeastAction" begin
        #include("LeastAction.jl")
    end
end


text:
function default_pre(feas::Bool,X::Vector{Interval{T}},UBD::Float64,
                     k::Int64,pos::Int64,opt,LBDn::Float64,UBDn::Float64,x::BnBSolver,
                     y::BnBModel{Interval{T}}) where {T<:AbstractFloat}
   return feas,X
end

function default_pre(feas::Bool,X::Vector{MCInterval{T}},UBD::Float64,
                     k::Int64,pos::Int64,opt,LBDn::Float64,UBDn::Float64,x::BnBSolver,
                     y::BnBModel{MCInterval{T}}) where {T<:AbstractFloat}
   return feas,X
end

function default_post(feas::Bool,X::Vector{Interval{T}},k::Int64,
                     pos::Int64,opt,tempL,tempU,LBD,UBD) where {T<:AbstractFloat}
    return feas,X
end

function default_post(feas::Bool,X::Vector{MCInterval{T}},k::Int64,
                     pos::Int64,opt,tempL,tempU,LBD,UBD) where {T<:AbstractFloat}
    return feas,X
end


text:
# Because resize! will not do the initialization of new elements for us
function my_grow!(f, vec::AbstractVector, newlen::Integer)
    diff = newlen - length(vec)
    diff == 0 && return vec
    @assert diff > 0
    sizehint!(vec, newlen)
    for i in 1:diff
        push!(vec, f())
    end
    @assert length(vec) == newlen
    return vec
end

my_grow!(vec::AbstractVector{T}, newlen::Integer, value::T=zero(T)) where {T} =
    my_grow!(() -> value, vec, newlen)

mutable struct StateInfo
    # Currently unused, except in tests.  It can be useful to have
    # this to be able to test if two states are in the same
    # equivalence class, though.  In the future we may be able to
    # optimize some things by knowing this.
    #
    # This gets initialized to 0, which denotes uninitialized.  The
    # RepresentativeStateTable constructor should set it to its final
    # value.
    representative_index::Int

    # If the Hamiltonian contains multiple sectors that do not mix when the
    # Hamiltonian and translation/symmetry operators are applied, they will be
    # available separately here.  Usually there will only be one sector, unless
    # our basis is mistakenly too large.
    #
    # This gets initialized to 0, which denotes uninitialized.  The
    # RepresentativeStateTable constructor should set it to its final
    # value.
    sector_index::Int

    # The results of translating in each direction, as well as results of any
    # additional symmetry operations provided.  The length of this vector is
    # the number of dimensions plus the number of additional symmetries given.
    transformation_results::Vector{Tuple{Int, Rational{Int}}}

    StateInfo() = new(0, 0, Tuple{Int, Rational{Int}}[])
end

# It's the job of this class to know all the representative states (and what
# sector each is in) as well as what state is obtained by translating each
# state in each direction (and applying each additional symmetry operation).
#
# We may eventually decide to make special optimized versions of this
# class for certain models.
#
# By default (if given no transformation_exponent_v), it assumes we can
# translate in every direction (except those with open boundary conditions) and
# apply each provided "additional symmetry" operation (e.g. spin flip).
struct RepresentativeStateTable{HilbertSpaceType<:HilbertSpace}
    # The Hilbert space and the function that applies the Hamiltonian
    #
    # FIXME: Instead of the following two objects separately, maybe we
    # should have an object that represents HilbertSpace +
    # Hamiltonian.
    hs::HilbertSpaceType
    apply_hamiltonian::Function

    # How many times we should apply each translation (or additional symmetry)
    # operator to advance to the next state.  This is a vector of length `dd`
    # (i.e. the number of dimensions plus the number of additional symmetries).
    # For open boundary conditions in a direction, set to 0.  For periodic, set
    # to 1.
    #
    # Currently this is assumed to have elements that are all either
    # zero or one.  In the future, we could make it so that e.g. if it
    # has value 2 in a given direction, then it is translated by two
    # sites, not one, when implementing the symmetry.  This could be
    # useful for a Hamiltonian where alternating sites have different
    # properties (and has an even number of sites with PBC).
    transformation_exponent_v::Vector{Int}

    # StateInfo for each state in hs.indexer
    state_info_v::Vector{StateInfo}

    # A list of indices of each representative state
    representative_state_indices::Vector{Int}

    # Total number of sectors (see above in StateInfo for explanation)
    sector_count::Int

    # The period of each "additional symmetry operation" that was provided
    # (e.g. spin flip has period 2, since two applications brings us back to
    # the original state)
    additional_symmetry_periods::Vector{Int}

    function RepresentativeStateTable{HilbertSpaceType}(hs::HilbertSpaceType, apply_hamiltonian::Function,
                                                        additional_symmetries::Vector{Tuple{F,Int}}=Tuple{Function,Int}[],
                                                        transformation_exponent_v::Vector{Int}=Int[]) where {HilbertSpaceType<:HilbertSpace,F<:Function}
        for (symm_func, symm_period) in additional_symmetries
            @assert symm_period > 0
        end

        d = ndimensions(hs.lattice)
        dd = d + length(additional_symmetries)

        # FIXME: move this logic to the outer constructor
        if isempty(transformation_exponent_v)
            # transformation by a single step allowed in each direction that is
            # not OBC
            transformation_exponent_v = [Int(i > d || repeater(hs.lattice)[i,i] != 0) for i in 1:dd]
        end

        @assert length(transformation_exponent_v) == dd

        @assert length(hs.indexer) > 0 # otherwise we have no seeded states!

        for i in 1:d
            @assert repeater(hs.lattice)[i, i] != 0 || transformation_exponent_v[i] == 0
        end
        # (nothing equivalent to check for d+1:dd)

        for i in 1:dd
            # see XXX below
            @assert transformation_exponent_v[i] == 0 || transformation_exponent_v[i] == 1
        end

        state_info_v = StateInfo[]
        representative_state_indices = Int[]

        CacheType = LatticeTranslationCache{typeof(hs.lattice)}
        ltrc = [transformation_exponent_v[i] != 0 ? Nullable{CacheType}(LatticeTranslationCache(hs.lattice, i)) : Nullable{CacheType}() for i in 1:d]

        transformation_basis_queue = UniqueVector{Int}()
        hamiltonian_basis_queue = UniqueVector{Int}()

        sector_count = 0
        for z in 1:length(hs.indexer)
            if z <= length(state_info_v) && state_info_v[z].representative_index != 0
                # We've already visited this
                continue
            end

            # We've discovered a new sector!
            sector_count += 1
            @assert isempty(hamiltonian_basis_queue)
            push!(hamiltonian_basis_queue, z)

            while !isempty(hamiltonian_basis_queue)
                y = pop!(hamiltonian_basis_queue)

                if y <= length(state_info_v) && state_info_v[y].representative_index != 0
                    # We've already visited this
                    continue
                end

                # We found a new representative state in this sector!
                push!(representative_state_indices, y)

                # Now try to translate in each direction (and perform each
                # additional symmetry operation) as long as we are finding new
                # states
                @assert isempty(transformation_basis_queue)
                push!(transformation_basis_queue, y)

                while !isempty(transformation_basis_queue)
                    x = pop!(transformation_basis_queue)

                    # We should not have already visited this
                    @assert !(x <= length(state_info_v) && state_info_v[x].representative_index != 0)

                    # Make sure our array is big enough, in case the indexer grew
                    my_grow!(StateInfo, state_info_v, length(hs.indexer))
                    @assert 0 < x <= length(state_info_v)

                    state_info_v[x].representative_index = y
                    state_info_v[x].sector_index = sector_count
                    resize!(state_info_v[x].transformation_results, dd) # NOTE: elements are uninitialized here!

                    # Translate in each direction, and perform each additional
                    # symmetry operation on the current state to generate
                    # states in the same class.
                    for i in 1:dd
                        if transformation_exponent_v[i] == 0
                            state_info_v[x].transformation_results[i] = (x, 0//1)
                        else
                            w = x
                            η = 0//1
                            for j in 1:transformation_exponent_v[i] # is always 0 or 1, for now.
                                w, η_inc = (i <= d) ? translateη(hs, get(ltrc[i]), w) : additional_symmetries[i-d][1](hs, w)
                                η += η_inc
                            end
                            if w > length(state_info_v) || state_info_v[w].representative_index == 0
                                findfirst!(isequal(w), transformation_basis_queue)
                            end
                            state_info_v[x].transformation_results[i] = (w, η)
                        end
                    end
                end

                # Now apply the hamiltonian to the current representative
                # state to generate other states
                #
                # NOTE: if we wanted to strictly check translation
                # invariance we would apply the hamiltonian to all states,
                # not just the representative ones.  so in the end we
                # really need to check that it is an eigenstate.
                apply_hamiltonian(hs, y) do newidx, amplitude
                    if amplitude != 0
                        if newidx > length(state_info_v) || state_info_v[newidx].representative_index == 0
                            findfirst!(isequal(newidx), hamiltonian_basis_queue)
                        end
                    end
                end
            end

        end

        # FIXME: If being careful, we might as well check that
        # representative_state_indices includes precisely those states
        # that are the representative of themselves, and no repeats.
        # And while we are at it, check that no sector_index or
        # representative_index is 0 anymore.

        additional_symmetry_periods = [symm_period for (symm_func, symm_period) in additional_symmetries]
        return new(hs, apply_hamiltonian, transformation_exponent_v, state_info_v, representative_state_indices, sector_count, additional_symmetry_periods)
    end
end

RepresentativeStateTable(hs::HilbertSpaceType, apply_hamiltonian::Function, additional_symmetries::Vector{Tuple{F,Int}} where {F<:Function} = Tuple{Function,Int}[], transformation_exponent_v::Vector{Int}=Int[]) where {HilbertSpaceType<:HilbertSpace} =
    RepresentativeStateTable{HilbertSpaceType}(hs, apply_hamiltonian, additional_symmetries, transformation_exponent_v)

# At times we will want to be able to specify which states are used as the
# representative ones (e.g. if we are loading the results of a previous
# calculation).  This constructor handles this case by immediately modifying
# the representative states once the state table has been constructed.
function RepresentativeStateTable(hs::HilbertSpace{StateType}, apply_hamiltonian::Function, representative_states::AbstractVector{StateType}, additional_args...) where {StateType}
    state_table = RepresentativeStateTable(hs, apply_hamiltonian, additional_args...)

    # Now we just need to set the representative states to those given

    # Plan the mutation
    old2new = Dict{Int,Int}() # maps the old representative index to the new one
    for state in representative_states
        newidx = findfirst(isequal(state), state_table.hs.indexer)
        oldidx = state_table.state_info_v[newidx].representative_index
        haskey(old2new, oldidx) && throw(ArgumentError("Two states were given in the same equivalence class"))
        old2new[oldidx] = newidx
    end
    for oldidx in state_table.representative_state_indices
        haskey(old2new, oldidx) || throw(ArgumentError("The states provided do not account for all representative states"))
    end

    # Perform the mutation
    for stateinfo in state_table.state_info_v
        stateinfo.representative_index = old2new[stateinfo.representative_index]
    end
    for (i, oldidx) in enumerate(state_table.representative_state_indices)
        state_table.representative_state_indices[i] = old2new[oldidx]
    end

    return state_table
end

function iterate_transformations(f, state_table::RepresentativeStateTable, z::Integer, bounds::Vector{Int})
    dd = length(bounds) # i.e. number of dimensions plus number of additional symmetry operations

    # NOTE: iteration_helper[dd] will *only* ever be incremented in the dd
    # direction.  iteration_helper[dd - 1] will be incremented both in the dd
    # direction and the dd-1 direction.  iteration_helper[1] is the current
    # element we want.  The other elements of the array are simply stored so we
    # can "wrap around" any time an index is reset to zero.
    iteration_helper = [(z, 0//1) for i in 1:dd]
    iter = zeros(Int, dd)

    while true
        # First, do what we need.
        f(iter, iteration_helper[1])

        # Second, advance the counter.
        i = 1
        while true
            i > dd && return # we are done iterating
            iter[i] += 1
            iter[i] != bounds[i] && break # we have advanced the counter
            iter[i] = 0 # index i gets wrapped around to zero
            i += 1
        end

        # Perform the actual translation (or symmetry operation) associated
        # with the new counter state.
        idx, η = iteration_helper[i]
        idx, η_inc = state_table.state_info_v[idx].transformation_results[i]
        iteration_helper[i] = (idx, η + η_inc)

        # For each index that wrapped around to zero when
        # advancing the counter, reset it to be the new
        # untransformed state
        for j in 1:i-1
            iteration_helper[j] = iteration_helper[i]
        end
    end
end

# Takes a RepresentativeStateTable, a momentum, and a sector index.
#
# Builds up normalization for each representative state, an
# UniqueVector with all such states of nonzero norm in our sector, and
# a mapping from each such non-representative state to a
# representative state index (with phase).
struct DiagonalizationSector{HilbertSpaceType<:HilbertSpace}
    # NOTE: All fields remains constant after initialization.

    state_table::RepresentativeStateTable{HilbertSpaceType}
    sector_index::Int
    momentum_index::Int

    # Includes only nonzero-normed representative states, and refers to states
    # by their index in the main indexer.
    reduced_indexer::UniqueVector{Int}

    # Contains the (nonzero) norm of each element in reduced_indexer
    norm_v::Vector{Float64}

    # Allows us to get the representative state \ket{r} and its
    # coefficient given a state \ket{s} (some element of \ket{r_k}).
    representative_v::Vector{Tuple{Int, ComplexF64}}

    # Allows us to explicitly construct \ket{r_k} from \ket{r}
    coefficient_v::Vector{Tuple{Int, Int, ComplexF64}}

    function DiagonalizationSector{HilbertSpaceType}(state_table::RepresentativeStateTable{HilbertSpaceType},
                                                     sector_index::Int,
                                                     momentum_index::Int,
                                                     reduced_indexer::UniqueVector{Int},
                                                     additional_symmetry_indices::Vector{Int}=Int[]) where {HilbertSpaceType<:HilbertSpace}
        @assert 0 < sector_index <= state_table.sector_count
        @assert 0 < momentum_index <= nmomenta(state_table.hs.lattice)

        @assert length(additional_symmetry_indices) == length(state_table.additional_symmetry_periods)
        @assert all(0 .<= additional_symmetry_indices .< state_table.additional_symmetry_periods)

        norm_v = Float64[]
        representative_v = Tuple{Int, ComplexF64}[(0, 0.0im) for i in 1:length(state_table.hs.indexer)]
        coefficient_v = Tuple{Int, Int, ComplexF64}[]

        d = ndimensions(state_table.hs.lattice)
        dd = d + length(state_table.additional_symmetry_periods)

        # This will become the number of distinct combinations of
        # transformation (i.e. translation + additional symmetry) steps we may
        # take, including the identity operation.  It is impicitly assumed that
        # all such operations commute with one another.
        transformation_count = 1

        bounds = zeros(Int, dd)
        for i in 1:dd
            if state_table.transformation_exponent_v[i] != 0
                @assert state_table.transformation_exponent_v[i] == 1 # enforce assumption (see XXX below)
                bounds[i] = (i <= d) ? repeater(state_table.hs.lattice)[i,i] : state_table.additional_symmetry_periods[i-d]
                transformation_count *= bounds[i]
            else
                bounds[i] = 1
            end
        end

        additional_symmetry_mult = [x//y for (x, y) in zip(additional_symmetry_indices, state_table.additional_symmetry_periods)]

        n_discovered_indices = 0

        for z in state_table.representative_state_indices
            if state_table.state_info_v[z].sector_index != sector_index
                continue
            end

            current_terms = Dict{Int, ComplexF64}()

            # XXX FIXME: make sure all the states in our sector have the same number!
            total_charge = get_total_charge(state_table.hs, z)
            total_momentum = momentum(state_table.hs.lattice, momentum_index, total_charge)

            iterate_transformations(state_table, z, bounds) do iter, current_transformation
                # XXX: assumes transformation_exponent_v contains only zeroes and ones

                # fixme: we can precalculate each k_dot_r
                # above and store it in a vector somehow,
                # though it may be tricky since r is a
                # "vector".  then again, we can just have a
                # size_t counter alongside, or use a
                # multi-dimensional array.
                idx, η = current_transformation
                kdr = kdotr(state_table.hs.lattice, total_momentum, iter[1:d])
                η += dot(iter[d+1:dd], additional_symmetry_mult)
                oldval = get!(current_terms, idx, complex(0.0))
                current_terms[idx] = oldval + cis(kdr + 2π * η) / transformation_count
            end

            normsq = mapreduce(abs2, +, values(current_terms))
            #println("$(state_table.hs.indexer[z])\t$(momentum_index)\t$(length(current_terms))\t$(normsq)")

            if normsq < 1e-8 # XXX: numerical tuning
                # this state doesn't exist in this momentum sector
                continue
            end

            reduced_i = findfirst!(isequal(z), reduced_indexer)
            n_discovered_indices += 1

            my_grow!(norm_v, length(reduced_indexer))
            norm = sqrt(normsq)
            norm_v[reduced_i] = norm

            # Save for later some things we have already calculated
            #
            # FIXME: make this optional if we are doing a huge system
            for (idx, val) in current_terms
                # NOTE: if we save these, we probably don't even need
                # norm_v.  in fact, either way it may not make sense
                # to store norm_v.  ugh, looks like we need it after all!
                #
                # XXX: rename these, and write down some algebra equations
                alpha = val / norm
                representative_v[idx] = (reduced_i, alpha)
                push!(coefficient_v, (reduced_i, idx, alpha))
            end
        end

        n_discovered_indices == length(reduced_indexer) || throw(ArgumentError("The provided reduced_indexer contains states that do not exist in this DiagonalizationSector."))

        # This is necessary for us to be able to perform a binary search using
        # searchsorted().
        sort!(coefficient_v)

        return new(state_table, sector_index, momentum_index, reduced_indexer, norm_v, representative_v, coefficient_v)
    end
end

DiagonalizationSector(state_table::RepresentativeStateTable{HilbertSpaceType}, sector_index::Int, momentum_index::Int, additional_symmetry_indices::Vector{Int}=Int[]) where {HilbertSpaceType<:HilbertSpace} =
    DiagonalizationSector{HilbertSpaceType}(state_table, sector_index, momentum_index, UniqueVector{Int}(), additional_symmetry_indices)

function DiagonalizationSector(state_table::RepresentativeStateTable{HilbertSpaceType}, sector_index::Int, momentum_index::Int, provided_reduced_indexer::AbstractVector{StateType}, additional_symmetry_indices::Vector{Int}=Int[]) where {StateType,HilbertSpaceType<:HilbertSpace{StateType}}
    indexer = state_table.hs.indexer
    reduced_indexer = UniqueVector{Int}()
    #sizehint!(reduced_indexer, length(provided_reduced_indexer))
    for state in provided_reduced_indexer
        i = findfirst!(isequal(state), indexer)
        if state_table.state_info_v[i].representative_index != i
            throw(ArgumentError("reduced_indexer contains states that our state_table has not chosen to be representative"))
        end
        push!(reduced_indexer, i)
    end

    diagsect = DiagonalizationSector{HilbertSpaceType}(state_table, sector_index, momentum_index, reduced_indexer, additional_symmetry_indices)

    # Make sure we didn't pick up any additional states that weren't
    # in our indexer before
    @assert length(reduced_indexer) == length(provided_reduced_indexer)
    @assert all(x -> x != 0, diagsect.norm_v)

    return diagsect
end

length(diagsect::DiagonalizationSector) = length(diagsect.reduced_indexer)
checkbounds(diagsect::DiagonalizationSector, i::Integer) = 0 < i <= length(diagsect) || throw(BoundsError())

function apply_reduced_commuting_operator(f, diagsect::DiagonalizationSector, reduced_j::Integer, apply_operator)
    # This, of course, assumes that the operator commutes with any
    # translations allowed by the diagsect.
    checkbounds(diagsect, reduced_j)
    j = diagsect.reduced_indexer[reduced_j]
    reduced_j_norm = diagsect.norm_v[reduced_j]
    apply_operator(diagsect.state_table.hs, j) do i, amplitude
        reduced_i, alpha = diagsect.representative_v[i]
        # If reduced_i is zero, this term has zero norm and is
        # therefore not in our sector, so we ignore it.  This is
        # expected!
        if reduced_i != 0
            @assert 0 < reduced_i <= length(diagsect)
            f(reduced_i, amplitude * conj(alpha) / reduced_j_norm)
        end
    end
end

function construct_reduced_commuting_operator(diagsect::DiagonalizationSector, apply_operator)
    s = length(diagsect)
    rows = Int[]
    cols = Int[]
    vals = ComplexF64[]
    for j in 1:s
        apply_reduced_commuting_operator(diagsect, j, apply_operator) do i, amplitude
            push!(rows, i)
            push!(cols, j)
            push!(vals, amplitude)
        end
    end

    return sparse(rows, cols, vals, s, s)
end

function apply_reduced_hamiltonian(f, diagsect::DiagonalizationSector, reduced_j::Integer)
    # This, of course, assumes that the Hamiltonian commutes with any
    # translations allowed by the diagsect.
    apply_reduced_commuting_operator(f, diagsect, reduced_j, diagsect.state_table.apply_hamiltonian)
end

function construct_reduced_hamiltonian(diagsect::DiagonalizationSector)
    hmat = construct_reduced_commuting_operator(diagsect, diagsect.state_table.apply_hamiltonian)

    # XXX: (numerical) assert hermitian
    #inhermiticity = sum(abs(hmat - hmat'))
    #if inhermiticity > 1e-4 # XXX numerical tuning
    #    println("WARNING: matrix is not hermitian by a value of $(inhermiticity)")
    #end
    #@assert inhermiticity < 1 # XXX numerical assert

    # FIXME: possibly return a Hermitian object once Julia supports sparse,
    # Hermitian matrices

    return hmat
end

function apply_reduced_operator(f, diagsect::DiagonalizationSector, reduced_j::Integer, apply_operator, args...)
    # Unlike apply_reduced_hamiltonian, this does *not* assume that the
    # operator commutes with the translation operators, so is therefore slower
    # to execute.
    checkbounds(diagsect, reduced_j)
    for z in searchsorted(diagsect.coefficient_v, reduced_j, by=(a -> a[1]))
        reduced_j_, j, alpha_j = diagsect.coefficient_v[z]
        apply_operator(diagsect.state_table.hs, j, args...) do i, amplitude
            reduced_i, alpha_i = diagsect.representative_v[i]
            # If reduced_i is zero, this term has zero norm and is
            # therefore not in our sector, so we ignore it.  This is
            # expected!
            if reduced_i != 0
                @assert 0 < reduced_i <= length(diagsect)
                f(reduced_i, amplitude * conj(alpha_i) * alpha_j)
            end
        end
    end
end

function construct_reduced_operator(diagsect::DiagonalizationSector, apply_operator, args...)
    s = length(diagsect)
    rows = Int[]
    cols = Int[]
    vals = ComplexF64[]
    for j in 1:s
        apply_reduced_operator(diagsect, j, apply_operator, args...) do i, amplitude
            push!(rows, i)
            push!(cols, j)
            push!(vals, amplitude)
        end
    end

    return sparse(rows, cols, vals, s, s)
end

function construct_reduced_indexer(rst::RepresentativeStateTable)
    # The typical use case for this function is to save precisely which states
    # the RST had in it, so we can later reconstruct it.  This way, if the code
    # later changes such that *other*, equivalent states would have been
    # considered representative, we can still reconstruct the RST, diagsect,
    # etc.  For this use case, the exact order of the representative states
    # does not matter; it simply matters what the representative states look
    # like.
    original_indexer = rst.hs.indexer
    return UniqueVector{eltype(original_indexer)}([original_indexer[i] for i in rst.representative_state_indices])
end

function construct_reduced_indexer(diagsect::DiagonalizationSector)
    original_indexer = diagsect.state_table.hs.indexer
    return UniqueVector{eltype(original_indexer)}([original_indexer[i] for i in diagsect.reduced_indexer])
end

function get_full_psi!(full_psi::Vector{ComplexF64}, diagsect::DiagonalizationSector, reduced_psi::AbstractVector)
    length(reduced_psi) == length(diagsect) || throw(DimensionMismatch("reduced_psi has the wrong length"))
    length(full_psi) == length(diagsect.state_table.hs.indexer) || throw(DimensionMismatch("full_psi has the wrong length"))
    fill!(full_psi, zero(ComplexF64))
    for (reduced_i, i, alpha) in diagsect.coefficient_v
        full_psi[i] = reduced_psi[reduced_i] * alpha
    end
    return full_psi
end

get_full_psi(diagsect::DiagonalizationSector, reduced_psi::AbstractVector) =
    get_full_psi!(Array{ComplexF64}(undef, length(diagsect.state_table.hs.indexer)), diagsect, reduced_psi)

# FIXME: these diagonalization functions are overkill for systems without symmetry.
#   1) they force the use of ComplexType
#   2) they involve lots of infrastructure and unnecessary build up.
#
# The solution?  For now, simply don't use them unless we have symmetry!


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

# Core functions for the usage in `/scripts`-folder

using DynamicalSystems
using RecurrenceAnalysis
using DelayEmbeddings
using Distances
using Statistics
using LinearAlgebra

"""
Perform the RecurrenceAnalysis of some reconstruction trajectories `Y₁`, `Y₂`,
`Y₃`, `Y₄`. Specifically, compute the fraction of recurrence rates from the
"original"/reference trajectory `Y_ref` with the one from the JRP of the
original `Y₁`, `Y₂`, `Y₃`, `Y₄` together with the reconstructed trajectory. Also
compute the recurrence time entropy of the recurrence plots of `Y₁`, `Y₂`, `Y₃`,
`Y₄` and `Y_ref`.

Keyword arguments:
*`ε = 0.05`: The used threshold for constructing the recurrence plots
    The reconstruction method is fixed recurrence rate.
*`w = 1`: Theiler window used for all Datasets
*`lmin = 2`: Minimum used line length for digaonal line based RQA measures.
*`kNN = 1`: The number of nearest neighbors used for obtaining the mutual
    nearest neighbors measure
"""
function perform_recurrence_analysis(Y_ref::Dataset, Y₁::Dataset,
                        Y₂::Dataset, Y₃::Dataset, Y₄::Dataset;
                        ε::Real = 0.05, w::Int = 1, lmin::Int = 2, kNN::Int = 1)

    N1 = length(Y₁)
    N2 = length(Y₂)
    N3 = length(Y₃)
    N4 = length(Y₄)
    N = minimum(hcat(N1, N2, N3, N4))

    R_ref = RecurrenceMatrix(Y_ref[1:N,:], ε; metric = "euclidean", fixedrate = true)
    R1 = RecurrenceMatrix(Y₁[1:N,:], ε; fixedrate = true)
    R2 = RecurrenceMatrix(Y₂[1:N,:], ε; fixedrate = true)
    R3 = RecurrenceMatrix(Y₃[1:N,:], ε; fixedrate = true)
    R4 = RecurrenceMatrix(Y₄[1:N,:], ε; fixedrate = true)

    f1 = jrp_rr_frac(R_ref, R1)
    f2 = jrp_rr_frac(R_ref, R2)
    f3 = jrp_rr_frac(R_ref, R3)
    f4 = jrp_rr_frac(R_ref, R4)

    mfnn1 = mfnn(Y_ref[1:N,:], Y₁[1:N,:]; w = w, kNN = kNN)
    mfnn2 = mfnn(Y_ref[1:N,:], Y₂[1:N,:]; w = w, kNN = kNN)
    mfnn3 = mfnn(Y_ref[1:N,:], Y₃[1:N,:]; w = w, kNN = kNN)
    mfnn4 = mfnn(Y_ref[1:N,:], Y₄[1:N,:]; w = w, kNN = kNN)

    RQA_ref = rqa(R_ref; theiler = w, lmin = lmin)
    RQA1 = rqa(R1; theiler = w, lmin = lmin)
    RQA2 = rqa(R2; theiler = w, lmin = lmin)
    RQA3 = rqa(R3; theiler = w, lmin = lmin)
    RQA4 = rqa(R4; theiler = w, lmin = lmin)

    return mfnn1, mfnn2, mfnn3, mfnn4, f1, f2, f3, f4, RQA_ref, RQA1, RQA2,
                                            RQA3, RQA4, R_ref, R1, R2, R3, R4
end

"""
    standard_embedding_hegger(s::Vector; kwargs...) → `Y`, `τ`
Compute the reconstructed trajectory from a time series using the standard time
delay embedding. The delay `τ` is taken as the 1st minimum of the mutual
information [`estimate_dimension`](@ref) and the embedding dimension `m` is
estimated by using an FNN method from [^Hegger1999] [`fnn_uniform_hegger`](@ref).
Return the reconstructed trajectory `Y` and the delay `τ`.

Keyword arguments:

*`fnn_thres = 0.05`: a threshold defining at which fraction of FNNs the search
    should break.
* The `method` can be one of the following:
* `"ac_zero"` : first delay at which the auto-correlation function becomes <0.
* `"ac_min"` : delay of first minimum of the auto-correlation function.
* `"mi_min"` : delay of first minimum of mutual information of `s` with itself
  (shifted for various `τs`). <- Default

[^Hegger1999]: Hegger, Rainer and Kantz, Holger (1999). [Improved false nearest neighbor method to detect determinism in time series data. Physical Review E 60, 4970](https://doi.org/10.1103/PhysRevE.60.4970).
"""
function standard_embedding_hegger(s::Vector{T}; method::String = "mi_min",
                                fnn_thres::Real = 0.05, τs = 1:200) where {T}
    @assert method=="ac_zero" || method=="mi_min" || method=="ac_min"
    "The absolute correlation function has elements that are = 0. "*
    "We can't fit an exponential to it. Please choose another method."

    τ = estimate_delay(s, method, τs)
    _, _, Y = fnn_uniform_hegger(s, τ; fnn_thres = fnn_thres)
    return Y, τ
end

"""
    standard_embedding_cao(s::Vector; kwargs...) → `Y`, `τ`
Compute the reconstructed trajectory from a time series using the standard time
delay embedding. The delay `τ` is taken as the 1st minimum of the mutual
information [`estimate_dimension`](@ref) and the embedding dimension `m` is
estimated by using an FNN method from Cao [`estimate_dimension`](@ref), with the
threshold parameter `cao_thres`.
Return the reconstructed trajectory `Y` and the delay `τ`.

Keyword arguments:
*`cao_thres = 0.05`: This threshold determines the tolerable deviation of the
    proposed statistic from the optimal value of 1, for breaking the algorithm.
*`m_max = 10`: The maximum embedding dimension, which is encountered by the
    algorithm.
* The `method` can be one of the following:
* `"ac_zero"` : first delay at which the auto-correlation function becomes <0.
* `"ac_min"` : delay of first minimum of the auto-correlation function.
* `"mi_min"` : delay of first minimum of mutual information of `s` with itself
  (shifted for various `τs`). <- Default

[^Hegger1999]: Hegger, Rainer and Kantz, Holger (1999). [Improved false nearest neighbor method to detect determinism in time series data. Physical Review E 60, 4970](https://doi.org/10.1103/PhysRevE.60.4970).
"""
function standard_embedding_cao(s::Vector{T}; cao_thres::Real = 0.05,
             τs = 1:200, method::String = "mi_min", m_max::Int = 10) where {T}
    @assert method=="ac_zero" || method=="mi_min" || method=="ac_min"
    "The absolute correlation function has elements that are = 0. "*
    "We can't fit an exponential to it. Please choose another method."

    τ = estimate_delay(s, method, τs)
    rat = estimate_dimension(s, τ, 1:m_max, "afnn")
    for i = 1:m_max
        if abs(1-rat[i]) < cao_thres
            global m = i
            break
        end
    end
    try
        if m > 1
            global Y = embed(s, m, τ)
        else
            global Y = s
        end
    catch
        global Y = s
    end
    return Y, τ, rat
end

"""
    fnn_uniform_hegger(s::Vector, τ::Int; kwargs...) →  `m`, `FNNs`, `Y`
Compute and return the optimal embedding dimension `m` for the time series `s`
and a uniform time delay `τ` after [^Hegger1999]. Return the optimal `m` and the
corresponding reconstruction vector `Y` according to that `m` and the input `τ`.
The optimal `m` is chosen, when the fraction of `FNNs` falls below the threshold
`fnn_thres` or when fraction of FNN's increases.

Keyword argument:
*`fnn_thres = 0.05`: Threshold, which defines the tolerable fraction of FNN's
    for which the algorithm breaks.
*`max_dimension = 10`: The maximum dimension which is encountered by the
    algorithm and after which it breaks, if the breaking criterion has not been
    met yet.
*`r = 2`: Obligatory threshold, which determines the maximum tolerable spreading
    of trajectories in the reconstruction space.
*`metric = Euclidean`: The norm used for distance computations.
*`w = 1` = The Theiler window, which excludes temporally correlated points from
    the nearest neighbor search.

[^Hegger1999]: Hegger, Rainer and Kantz, Holger (1999). [Improved false nearest neighbor method to detect determinism in time series data. Physical Review E 60, 4970](https://doi.org/10.1103/PhysRevE.60.4970).
"""
function fnn_uniform_hegger(s::Vector{T}, τ::Int; max_dimension::Int = 10,
            r::Real = 2, w::Int = 1, fnn_thres::Real = 0.05, metric = Euclidean()) where {T}
    @assert max_dimension > 0
    s = (s .- mean(s)) ./ std(s)
    Y_act = s

    vtree = KDTree(Dataset(s), metric)
    _, NNdist_old = DelayEmbeddings.all_neighbors(vtree, Dataset(s), 1:length(s), 1, w)

    FNNs = zeros(max_dimension)
    for m = 2:max_dimension+1
        Y_act = DelayEmbeddings.hcat_lagged_values(Y_act, s, m*τ)
        Y_act = regularize(Y_act)
        vtree = KDTree(Y_act, metric)
        _, NNdist_new = DelayEmbeddings.all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

        FNNs[m-1] = DelayEmbeddings.fnn_embedding_cycle(view(NNdist_old,
                                            1:length(Y_act)), NNdist_new, r)

        flag = fnn_break_criterion(FNNs[1:m-1], fnn_thres)
        if flag
            global bm = m
            break
        else
            global bm = m
        end

        NNdist_old = NNdist_new
    end

    if bm>2
        Y_final = embed(s, bm-1, τ)
    else
        Y_final = s
    end
    return bm, FNNs[1:bm-1], Y_final
end

"""
Determines the break criterion for the Hegger-FNN-estimation
"""
function fnn_break_criterion(FNNs, fnn_thres)
    flag = false
    if FNNs[end] ≤ fnn_thres
        flag = true
        println("Algorithm stopped due to sufficiently small FNNs. "*
                "Valid embedding achieved ✓.")
    end
    if length(FNNs) > 1 && FNNs[end] > FNNs[end-1]
        flag = true
        println("Algorithm stopped due to rising FNNs. "*
                "Valid embedding achieved ✓.")
    end
    return flag
end

"""
Computes the similarity between recurrence plots `RP₁` and `RP₂`. Outputs the
fraction of recurrences rates gained from RP₁ and of the joint recurrence
plot `RP₁ .* RP₂`.
"""
function jrp_rr_frac(RP₁::RecurrenceMatrix, RP₂::RecurrenceMatrix)
    @assert size(RP₁) == size(RP₂)

    RR1 = sum(RP₁)/(size(RP₁,1)*size(RP₁,1))
    JRP = elementwise_product(RP₁, RP₂)
    RR2 = sum(JRP)/(size(JRP,1)*size(JRP,1))

    f = RR2 / RR1
    return f
end

"""
Computes the mututal false nearest neighbours (mfnn) for a reference trajectory
`Y_ref` and a reconstruction `Y_rec` after [^Rulkov1995].

Keyword arguments:

*`w = 1`: Theiler window for the surpression of serially correlated neighbors in
    the nearest neighbor-search
*`kNN = 1`: The number of considered nearest neighbours (in the paper always 1)

[^Rulkov1995]: Rulkov, Nikolai F. and Sushchik, Mikhail M. and Tsimring, Lev S. and Abarbanel, Henry D.I. (1995). [Generalized synchronization of chaos in directionally coupled chaotic systems. Physical Review E 51, 980](https://doi.org/10.1103/PhysRevE.51.980).
"""
function mfnn(Y_ref::Dataset, Y_rec::Dataset; w::Int = 1, kNN::Int = 1)

    @assert length(Y_ref) == length(Y_rec)
    @assert kNN > 0
    N = length(Y_ref)
    metric = Euclidean()

    # compute nearest neighbor distances for both trajectories
    vtree = KDTree(Y_ref, metric)
    allNNidxs_ref, _ = DelayEmbeddings.all_neighbors(vtree, Y_ref,
                                                        1:length(Y_ref), kNN, w)
    vtree = KDTree(Y_rec, metric)
    allNNidxs_rec, _ = DelayEmbeddings.all_neighbors(vtree, Y_rec,
                                                        1:length(Y_rec), kNN, w)

    F = zeros(N)
    factor1_nom = zeros(kNN)
    factor1_denom = zeros(kNN)
    factor2_nom = zeros(kNN)
    factor2_denom = zeros(kNN)
    for i = 1:N
        for j = 1:kNN
            factor1_nom[j] = evaluate(Euclidean(), Y_rec[i], Y_rec[allNNidxs_ref[i][j]])
            factor1_denom[j] = evaluate(Euclidean(), Y_ref[i], Y_ref[allNNidxs_ref[i][j]])
            factor2_nom[j] = evaluate(Euclidean(), Y_ref[i], Y_ref[allNNidxs_rec[i][j]])
            factor2_denom[j] = evaluate(Euclidean(), Y_rec[i], Y_rec[allNNidxs_rec[i][j]])
        end
        factor1 = sum(factor1_nom)/sum(factor1_denom)
        factor2 = sum(factor2_nom)/sum(factor2_denom)
        F[i] = factor1*factor2                                         # Eq.(27)
    end
    return mean(F)
end

"""
Compute elementwise product of two Recurrence Plots (==JRP)
"""
elementwise_product(RP₁, RP₂) = JointRecurrenceMatrix(RP₁, RP₂)

"""
Generate data from a AR(1) process for a initial value `u0`, a AR-coefficient
`α` and a white noise scaling parameter `p`. Return a time series of length `N`.
"""
function ar_process(u0::T, α::T, p::T, N::Int) where {T<:Real}
    x = zeros(T, N+10)
    x[1] = u0
    for i = 2:N+10
        x[i] = α*x[i-1] + p*randn()
    end
    return x[11:end]
end

"""
Generate data from a AR(2) process for a initial value `u0`, AR-coefficients
`α₁` and `α₂` and a white noise scaling parameter `p`. Return a time series of
length `N`.
"""
function ar_process2(u0::Vector, α₁::T, α₂::T, p::T, N::Int) where {T<:Real}
    @assert (length(u0) == 2)
    x = zeros(T, N+100)
    x[1] = u0[1]
    x[2] = u0[1]
    for i = 3:N+10
        x[i] = α₂*x[i-2] + α₁*x[i-1] + p*randn()
    end
    return x[101:end]
end

"""
Generate data from a AR(3) process for a initial value `u0`, AR-coefficients
`α₁`, `α₂` and `α₃` and a white noise scaling parameter `p`. Return a time series of
length `N`.
"""
function ar_process3(u0::Vector, α₁::T, α₂::T, α₃::T, p::T, N::Int) where {T<:Real}
    @assert (length(u0) == 3)
    x = zeros(T, N+100)
    x[1] = u0[1]
    x[2] = u0[2]
    x[3] = u0[3]
    for i = 4:N+10
        x[i] = α₃*x[i-3] + α₂*x[i-2] + α₁*x[i-1] + p*randn()
    end
    return x[101:end]
end

"""
Generate data from a AR(5) process for a initial value `u0`, AR-coefficients
`α₁`, `α₂`, `α₃`, `α₄`, `α₅` and a white noise scaling parameter
`p`. Return a time series of
length `N`.
"""
function ar_process5(u0::Vector, α₁::T, α₂::T, α₃::T, α₄::T, α₅::T, p::T, N::Int) where {T<:Real}
    @assert (length(u0) == 5)
    x = zeros(T, N+100)
    x[1] = u0[1]
    x[2] = u0[2]
    x[3] = u0[3]
    x[4] = u0[4]
    x[5] = u0[5]
    for i = 6:N+10
        x[i] = α₅*x[i-5] + α₄*x[i-4] + α₃*x[i-3] + α₂*x[i-2] + α₁*x[i-1] + p*randn()
    end
    return x[101:end]
end

"""
Generate data from a AR(7) process for a initial value `u0`, AR-coefficients
`α₁`, `α₂`, `α₃`, `α₄`, `α₅`, `α₆` and `α₇` and a white noise scaling parameter
`p`. Return a time series of
length `N`.
"""
function ar_process7(u0::Vector, α₁::T, α₂::T, α₃::T, α₄::T, α₅::T, α₆::T, α₇::T, p::T, N::Int) where {T<:Real}
    @assert (length(u0) == 7)
    x = zeros(T, N+100)
    x[1] = u0[1]
    x[2] = u0[2]
    x[3] = u0[3]
    x[4] = u0[4]
    x[5] = u0[5]
    x[6] = u0[6]
    x[6] = u0[7]
    for i = 8:N+10
        x[i] = α₇*x[i-7] + α₆*x[i-6] + α₅*x[i-5] + α₄*x[i-4] + α₃*x[i-3] + α₂*x[i-2] + α₁*x[i-1] + p*randn()
    end
    return x[101:end]
end

"""
compute_delta_L(s, τs, (js,) T_max; KNN = 3, w = 1, metric = Euclidean) → ΔL

    Compute the overall L-decrease `ΔL` of a given embedding of the time series
`s::Vector` with the delay values `τs` up to a maximum `T`-value `T_max`. We
respect the Theiler window `w`, the chosen `metric` and the number of considered
nearest neighbors `KNN`. It is also possible to compute `ΔL` for a multivariate
input `Y::Dataset`. Then one additionally needs to supply a vector `js`, which
lists the chosen time series corresponding to the given delay values in `τs`.
This is similar to the procedure in [`genembed`]@ref. The computations are based
on z-standardized input for ensuring comparability.
"""
function compute_delta_L(s::Vector{T}, τs::Vector{Int}, T_max::Int;
        KNN::Int = 3, w::Int = 1, metric = Euclidean()
    ) where {T}
    ts = regularize(Dataset(s))
    tss = vec(Matrix(Dataset(ts)))
    ΔL = 0
    Y_act = ts
    for i = 1 : (length(τs)-1)
        τ = τs[i+1]
        # embedding one cycle
        Y_next = DelayEmbeddings.hcat_lagged_values(Y_act, tss, τ)
        # compute ΔL for this cycle
        ΔL += uzal_cost_pecuzal(Y_act, Y_next, T_max; K = KNN, w = w, metric = metric)
        Y_act = Y_next
    end
    return ΔL
end

function compute_delta_L(Y::Dataset{D, T}, τs::Vector{Int}, js::Vector{Int},
        T_max::Int; KNN::Int = 3, w::Int = 1, metric = Euclidean()) where {D, T}
    @assert length(τs) == length(js)
    ts = regularize(Y)
    ΔL = 0
    Y_act = Dataset(ts[:,js[1]])
    for i = 1: (length(τs)-1)
        τ = τs[i+1]
        # embedding one cycle
        Y_next = DelayEmbeddings.hcat_lagged_values(Y_act, ts[:,js[i+1]], τ)
        # compute ΔL for this cycle
        ΔL += uzal_cost_pecuzal(Y_act, Y_next, T_max; K = KNN, w = w,
        metric = metric)
        Y_act = Y_next
    end
    return ΔL
end


text:
module MICE
using GLM, DataFrames
import Random
import StatsModels
using Statistics
import GLM: predict

struct MICEmodel
   varorder
   p1
   regmodels
end

struct PredictSameModel
   value
end

predict(p::PredictSameModel, x) = repeat([p.value], size(x,1))

function formula_n(n, onevaluevariables)
   diversevariables = .!(onevaluevariables[1:(n-1)])
   StatsModels.FormulaTerm(Term(Symbol("x" * string(n))),
         Tuple(map( i -> Term(Symbol("x" * string(i))), (1:(n-1))[diversevariables])))
end

function fitmice(x::Matrix{Float64})
   nsamples, nvariables = size(x)

   varorder = Random.shuffle(1:nvariables)
   df = DataFrame(x[:, varorder])
   onevaluevariables = map(i -> reduce(==, df[:,i]), 1:nvariables)
   colinearvariables = map(i -> any(map(j -> df[:,i] == df[:,j], 1:(i-1))), 1:nvariables)
   p1 = mean(df[:,1])
   regmodels = Array{Any, 1}(undef, nvariables-1)
   for i in 2:nvariables
      if onevaluevariables[i-1] == true
         regmodels[i-1] = PredictSameModel(df[1,i-1])
      else
         formula = formula_n(i, onevaluevariables .| colinearvariables)
         try
            regmodels[i-1] = glm(formula, df, Bernoulli(), LogitLink())
         catch ex
            colinearvariables[i-1] = true
            formula = formula_n(i, onevaluevariables .| colinearvariables)
            regmodels[i-1] = glm(formula, df, Bernoulli(), LogitLink())
         end
      end
   end
   MICEmodel(varorder, p1, regmodels)
end

function genmice(m::MICEmodel, nsamples)
   x1 = float.(rand(nsamples) .< m.p1)
   df = DataFrame(x1 = x1)
   nvariables = length(m.varorder)
   for i in 2:nvariables
      df[:,i] = round.(rand(nsamples) .< predict(m.regmodels[i-1], df[:,1:(i-1)]))
   end
   df = df[:,invperm(m.varorder)]
   convert(Array, df)
end

end

max_stars_repo_name: stefan-m-lenz/dist-gen-comp

max_issues_repo_name: stefan-m-lenz/dist-gen-comp

max_forks_repo_name: stefan-m-lenz/dist-gen-comp


text:
@testset "807.max-increase-to-keep-city-skyline.jl" begin
    @test max_increase_keeping_skyline(hcat([[3,0,8,4],[2,4,5,7],[9,2,6,3],[0,3,1,0]]...)) == 35
end


text:
test_gage_principal_axes(
    [1. 1. 0.;
     0. 0. 0.;],
    [sqrt(2)/2. -sqrt(2)/2. 0.0;
     sqrt(2)/2. sqrt(2)/2. 0.0;
     0.0 0.0 1.0], "vector in the X-Y plane")


text:
"""
Calls type inference on a function type.
"""
function Base.code_typed{T<:Function}(functype::Type{T}, argtypes...)
    # Extract a type-inferred AST:
    m = Base._methods_by_ftype(Tuple{functype, argtypes...}, -1)[1]
    codeinfo = Core.Inference.typeinf_uncached(m[3], m[1], m[2], true)[1]
    ast = Base.uncompressed_ast(m[3], codeinfo)
end

function sexpr(obj)
    io = IOBuffer()
    Base.Meta.show_sexpr(io, obj)
    takebuf_string(io)
end

""" Depth-first tree traversal and transformation. """
function depthfirst
end
function depthfirst(func::Function, e::Expr)
    new_args = [depthfirst(func, arg) for arg in e.args]
    new_expr = Expr(e.head, new_args...)
    func(new_expr)
end
function depthfirst(func::Function, code::CodeInfo)
    newcode = deepcopy(code)
    newcode.code = Any[depthfirst(func, stmt) for stmt in code.code]
    # must be Any[] otherwise the CodeInfo will show as <compressed>.
    func(newcode)
end
depthfirst(func::Function, o::Any) = func(o)


text:
"""
    Module Basis2DQuad

Includes DG core functions.
"""

module Basis2DQuad

export vandermonde_2D, grad_vandermonde_2D
export nodes_2D, equi_nodes_2D, quad_nodes_2D

using Basis1D
using CommonUtils
using LinearAlgebra

"""
    vandermonde_2D(N, r)

Initialize the 2D Vandermonde matrix of order N "Legendre" polynomials at
nodes (r,s)

# Examples
```jldoctest
"""
function vandermonde_2D(N, r, s)

    Np = convert(Int,(N+1)*(N+1))
    sk = 1
    V = zeros(length(r), Np);
    for i=0:N
        for j=0:N
            V[:,sk] = jacobiP(r, 0, 0, i).*jacobiP(s, 0, 0, j)
            sk += 1
        end
    end
    return V
end

"""
    gradV_2D(N, Np, r, s)

Initilize 2D gradient Vandermonde matrix of derivatives 2D "Legendre" polynomials
of order N at (r,s)

# Examples
```jldoctest
"""
function grad_vandermonde_2D(N, r, s)

    Np = convert(Int,(N+1)*(N+1))
    sk = 1
    V2Dr = zeros(length(r), Np);
    V2Ds = zeros(length(r), Np);
    for i=0:N
        for j=0:N
            V2Dr[:,sk] = grad_jacobiP(r, 0, 0, i).*jacobiP(s, 0, 0, j)
            V2Ds[:,sk] = jacobiP(r, 0, 0, i).*grad_jacobiP(s, 0, 0, j)
            sk += 1
        end
    end

    return V2Dr, V2Ds
end

# ===================================================

"""
    nodes_2D(N)

Compute optimized interpolation nodes using blend & warp method on equilateral
triangles for polynomial of order N, with Np points

# Examples
```jldoctest
"""

function nodes_2D(N)
    r1D,w1D = gauss_lobatto_quad(0,0,N)
    s,r = meshgrid(r1D)
    return r[:], s[:]
end

"""
    equi_nodes_2D(N)

Compute optimized interpolation nodes using blend & warp method on equilateral
triangles for polynomial of order N, with Np points

# Examples
```jldoctest
"""

function equi_nodes_2D(N)
    r1D = LinRange(-1,1,N+1)
    s,r = meshgrid(r1D)

    return r[:], s[:]
end

"""
    quad_nodes_2D(N)

Compute optimized interpolation nodes using blend & warp method on equilateral
triangles for polynomial of order N, with Np points

# Examples
```jldoctest
"""

function quad_nodes_2D(N)
    r1D,w1D = gauss_quad(0,0,N)
    s,r = meshgrid(r1D)
    ws,wr = meshgrid(w1D)
    w = @. wr*ws
    return r[:], s[:], w[:]
end

end

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
abstract type AbstractCovariate end

struct Linear <: AbstractCovariate end
covariate(::Linear, X) = [ones(1, size(X, 2)); X]

struct Constant <: AbstractCovariate end
covariate(::Constant, X) = ones(eltype(X), 1, size(X, 2))

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# These are functions and structs taken from Luxor.jl: https://github.com/JuliaGraphics/Luxor.jl
# We couldn't use Luxor.jl directly since it doesn't compile on Arm64 so far.

struct Point
    x::AbstractFloat
    y::AbstractFloat
end

mutable struct BoundingBox
    corner1::Point
    corner2::Point
 end

function BoundingBox(pointlist::Array{Point, 1})
    lowx, lowy = pointlist[1].x, pointlist[1].y
    highx, highy = pointlist[end].x, pointlist[end].y
    for p in pointlist
        p.x < lowx  && (lowx  = p.x)
        p.y < lowy  && (lowy  = p.y)
        p.x > highx && (highx = p.x)
        p.y > highy && (highy = p.y)
    end
    return [Point(lowx, lowy), Point(highx, highy)]
end

function det3p(q1::Point, q2::Point, p::Point)
    (q1.x - p.x) * (q2.y - p.y) - (q2.x - p.x) * (q1.y - p.y)
end

function isinside(p::Point, pointlist::Array{Point, 1};
    allowonedge::Bool=false)
c = false
@inbounds for counter in eachindex(pointlist)
    q1 = pointlist[counter]
    # if reached last point, set "next point" to first point
    if counter == length(pointlist)
        q2 = pointlist[1]
    else
        q2 = pointlist[counter + 1]
    end
    if q1 == p
        allowonedge || error("isinside(): VertexException a")
        continue
    end
    if q2.y == p.y
        if q2.x == p.x
            allowonedge || error("isinside(): VertexException b")
            continue
        elseif (q1.y == p.y) && ((q2.x > p.x) == (q1.x < p.x))
            allowonedge || error("isinside(): EdgeException")
            continue
        end
    end
    if (q1.y < p.y) != (q2.y < p.y) # crossing
        if q1.x >= p.x
            if q2.x > p.x
                c = !c
            elseif ((det3p(q1, q2, p) > 0) == (q2.y > q1.y))
                c = !c
            end
        elseif q2.x > p.x
            if ((det3p(q1, q2, p) > 0) == (q2.y > q1.y))
                c = !c
            end
        end
    end
end
return c
end


text:
using Test

@testset "example1" begin
    include("example1.jl")
end


text:
using DataFrames
using NveData
using IcanProj
using VannModels
using StatsBase

import IcanProj.read_vic_forcing

struct SimRes

    swe_error
    swe_fine
    swe_coarse
    elev
    xcoord
    ycoord

end

function get_simulation_domain(res, opt)

    # Indices of coarse resolution grid boxes

    senorge_ind, xcoord, ycoord = senorge_info()

    nrow, ncol = size(senorge_ind)

    boxes_ind = resolution_ind(res, nrow, ncol)

    df_boxes = DataFrame(senorge_ind = senorge_ind[:],
                         boxes_ind = boxes_ind[:],
                         xcoord = xcoord[:],
                         ycoord = ycoord[:])

    # Information from soil parameter file

    soil_param = read_soil_params(opt["base_folder"])

    df_soil = DataFrame(senorge_ind = soil_param[:gridcel] + 1,
                        lat = soil_param[:lat],
                        lon = soil_param[:lon],
                        elev = soil_param[:elev])

    # Keep cells for simulations

    df = join(df_soil, df_boxes, on = :senorge_ind, kind = :left)

    iremove = find(df[:boxes_ind] .== -9999)

    nkeep = floor(0.8 * res^2)

    for ibox in 1:maximum(df[:boxes_ind])

        isel = find(df[:boxes_ind] .== ibox)

        if length(isel) <= nkeep
            append!(iremove, isel)
        elseif length(isel) > nkeep
            nremove = Int(length(isel) - nkeep)
            append!(iremove, sample(isel, nremove, replace=false))
        end

    end

    deleterows!(df, iremove)

    return df

end

function read_vic_forcing(path_forcing, lat, lon)

    lat = @sprintf("%0.5f", lat)

    lon = @sprintf("%0.5f", lon)

    file_src = joinpath(path_forcing, "forcing/data_$(lat)_$(lon)")

    prec, tmin, tmax, wind = read_vic_forcing(file_src)

    return prec, tmin, tmax, wind

end

function snow_sim!(swe, prec, tair)

    mdata = TinSnow(24.0, DateTime(2000,1,1), DataFrame(dummy=1))

    for i in 1:length(prec)

        mdata.p_in[1] = prec[i]
        mdata.tair[1] = tair[i]

        run_timestep(mdata)

        swe[i] += mdata.swe[1,1]

    end

    return nothing

end

function run_single_box(ibox, path_forcing, df)

    # Initilize output arrays

    prec, tmin, tmax, wind = read_vic_forcing(path_forcing, df[:lat][1], df[:lon][1])

    swe_fine = fill(0.0, size(prec))

    swe_coarse = fill(0.0, size(prec))

    tmean_coarse = fill(0.0, size(prec))

    pmean_coarse = fill(0.0, size(prec))

    # Loop over cells in one box

    for icell in ibox

        # Load forcings

        prec, tmin, tmax, wind = read_vic_forcing(path_forcing, df[:lat][icell], df[:lon][icell])

        tmean = (tmin + tmax)/2

        tmean_coarse += tmean
        pmean_coarse += prec

        # Snow simulation for fine resolution

        snow_sim!(swe_fine, prec, tmean)

    end

    # Average fine scale results

    swe_fine /= length(ibox)
    tmean_coarse /= length(ibox)
    pmean_coarse /= length(ibox)

    # Snow simulation for coarse resolution

    snow_sim!(swe_coarse, pmean_coarse, tmean_coarse)

    # Compute error between fine and coarse resolution

    swe_error = rmsd(swe_coarse, swe_fine) / mean(swe_fine)

    # Collect results

    elev = df[:elev][ibox]
    xcoord = mean(df[:xcoord][ibox])
    ycoord = mean(df[:ycoord][ibox])

    res = SimRes(swe_error,
                 swe_fine,
                 swe_coarse,
                 elev,
                 xcoord,
                 ycoord)

    return res

end

# Get information about simulation domain

res = 50

opt = get_options()

df = get_simulation_domain(res, opt)

# Snow simulations on different resolutions

path_forcing = joinpath(opt["base_folder"])

boxes_unique = unique(df[:boxes_ind])

res = SimRes[]

for i in boxes_unique

    @show i

    ibox = find(df[:boxes_ind] .== i)

    push!(res, run_single_box(ibox, path_forcing, df))

end

# Save results