text:
""" Low-level wrapper of C library Serd.
"""
module CSerd
export SerdException, SerdNode, SerdStatement, SerdStatementFlags, SerdStyles,
  SerdReader, serd_reader_new, serd_reader_free,
  serd_reader_set_error_sink, serd_reader_set_strict,
  serd_reader_add_blank_prefix, serd_reader_set_default_graph,
  serd_reader_read_file, serd_reader_read_string,
  SerdWriter, serd_writer_new, serd_writer_free,
  serd_writer_set_error_sink, serd_writer_chop_blank_prefix,
  serd_writer_set_base_uri, serd_writer_set_root_uri,
  serd_writer_set_prefix, serd_writer_write_statement, serd_writer_finish

using AutoHashEquals

import Serd_jll
const serd = Serd_jll.libserd

const Cbool = UInt8

""" Export an enum and all its values.
"""
macro export_enum(name::Symbol)
  expr = :(eval(Expr(:export, $(QuoteNode(name)),
                     (Symbol(inst) for inst in instances($name))...)))
  esc(Expr(:toplevel, expr))
end

# Data types
############

@enum SerdStatus::Cint begin
  SERD_SUCCESS
  SERD_FAILURE
  SERD_ERR_UNKNOWN
  SERD_ERR_BAD_SYNTAX
  SERD_ERR_BAD_ARG
  SERD_ERR_NOT_FOUND
  SERD_ERR_ID_CLASH
  SERD_ERR_BAD_CURIE
  SERD_ERR_INTERNAL
end
@export_enum SerdStatus

@enum SerdSyntax::Cint begin
  SERD_TURTLE   = 1  # Turtle - Terse RDF Triple Language
  SERD_NTRIPLES = 2  # NTriples - Line-based RDF triples
  SERD_NQUADS   = 3  # NQuads - Line-based RDF quads
  SERD_TRIG     = 4  # TRiG - Terse RDF quads
end
@export_enum SerdSyntax

@enum SerdStatementFlag::Cint begin
  SERD_EMPTY_S      = 1 << 1  # Empty blank node subject
  SERD_EMPTY_O      = 1 << 2  # Empty blank node object
  SERD_ANON_S_BEGIN = 1 << 3  # Start of anonymous subject
  SERD_ANON_O_BEGIN = 1 << 4  # Start of anonymous object
  SERD_ANON_CONT    = 1 << 5  # Continuation of anonymous node
  SERD_LIST_S_BEGIN = 1 << 6  # Start of list subject
  SERD_LIST_O_BEGIN = 1 << 7  # Start of list object
  SERD_LIST_CONT    = 1 << 8  # Continuation of list
end
@export_enum SerdStatementFlag

const SerdStatementFlags = Cint
Base.convert(::Type{SerdStatementFlags}, x::SerdStatementFlag) = Cint(x)

@enum SerdType::Cint begin
  SERD_NOTHING = 0  # The type of a nonexistent node.
  SERD_LITERAL = 1  # Literal value
  SERD_URI     = 2  # URI (absolute or relative)
  SERD_CURIE   = 3  # CURIE, a shortened URI
  SERD_BLANK   = 4  # A blank node
end
@export_enum SerdType

@enum SerdStyle::Cint begin
  SERD_STYLE_ABBREVIATED = 1       # Abbreviate triples when possible
  SERD_STYLE_ASCII       = 1 << 1  # Escape all non-ASCII characters
  SERD_STYLE_RESOLVED    = 1 << 2  # Resolve URIs against base URI
  SERD_STYLE_CURIED      = 1 << 3  # Shorten URIs into CURIEs
  SERD_STYLE_BULK        = 1 << 4  # Write output in pages
end
@export_enum SerdStyle

const SerdStyles = Cint
Base.convert(::Type{SerdStyles}, x::SerdStyle) = Cint(x)

struct SerdException <: Exception
  status::SerdStatus
end

@auto_hash_equals struct SerdNode
  value::String
  typ::SerdType
end

@auto_hash_equals struct SerdStatement
  flags::SerdStatementFlags
  graph::Union{SerdNode,Nothing}
  subject::SerdNode
  predicate::SerdNode
  object::SerdNode
  object_datatype::Union{SerdNode,Nothing}
  object_lang::Union{SerdNode,Nothing}
end

mutable struct SerdReader
  ptr::Ptr{Cvoid}
  base_sink::Union{Function,Nothing}
  prefix_sink::Union{Function,Nothing}
  statement_sink::Union{Function,Nothing}
  end_sink::Union{Function,Nothing}
  error_sink::Union{Function,Nothing}
end

mutable struct SerdWriter
  ptr::Ptr{Cvoid}
  env::Ptr{Cvoid}
  sink::Function
  error_sink::Union{Function,Nothing}
end

mutable struct CSerdNode
  buf::Ptr{UInt8}
  n_bytes::Csize_t
  n_chars::Csize_t
  flags::Cint
  typ::Cint
end

mutable struct CSerdError
  status::Cint
  filename::Ptr{UInt8}
  line::Cuint
  col::Cuint
  char::Ptr{Cchar}
  args::Ptr{Cvoid} # va_list *
end

# Node
######

""" Convert Serd node from Julia struct to C struct.
"""
function c_serd_node(node::SerdNode)::CSerdNode
  ccall((:serd_node_from_string, serd), CSerdNode, (Cint, Cstring),
        Cint(node.typ), node.value)
end
function c_maybe_serd_node(node::Union{SerdNode,Nothing})::Ptr{CSerdNode}
  isnothing(node) ? C_NULL : pointer_from_objref(c_serd_node(node))
end

""" Convert Serd node from C struct to Julia struct.
"""
function unsafe_serd_node(ptr::Ptr{CSerdNode})::Union{SerdNode,Nothing}
  if ptr == C_NULL
    nothing
  else
    typ = unsafe_load(Ptr{Cint}(ptr + fieldoffset(CSerdNode,5)))
    if typ == SERD_NOTHING
      nothing
    else
      n_bytes = unsafe_load(Ptr{Csize_t}(ptr + fieldoffset(CSerdNode,2)))
      value_ptr = unsafe_load(Ptr{Ptr{UInt8}}(ptr))
      value = unsafe_string(value_ptr, n_bytes)
      SerdNode(value, SerdType(typ))
    end
  end
end

serd_status(status) = Cint(isa(status, SerdStatus) ? status : SERD_SUCCESS)

function check_serd_status(status)
  status = SerdStatus(status)
  if status != SERD_SUCCESS
    throw(SerdException(status))
  end
end

# Reader
########

""" Create a new RDF reader.
"""
function serd_reader_new(syntax::SerdSyntax, base_sink, prefix_sink,
                         statement_sink, end_sink)::SerdReader
  serd_base_sink_ptr = @cfunction(
    $(CSerd.serd_base_sink), Cint, (Ptr{Cvoid}, Ptr{CSerdNode}))
  serd_prefix_sink_ptr = @cfunction(
    $(CSerd.serd_prefix_sink), Cint, (Ptr{Cvoid}, Ptr{CSerdNode}, Ptr{CSerdNode}))
  serd_statement_sink_ptr = @cfunction(
    $(CSerd.serd_statement_sink),
    Cint,
    (Ptr{Cvoid}, Cint, Ptr{CSerdNode}, Ptr{CSerdNode}, Ptr{CSerdNode},
     Ptr{CSerdNode}, Ptr{CSerdNode}, Ptr{CSerdNode}))
  serd_end_sink_ptr = @cfunction(
    $(CSerd.serd_end_sink), Cint, (Ptr{Cvoid}, Ptr{CSerdNode}))

  reader = SerdReader(
    C_NULL, base_sink, prefix_sink, statement_sink, end_sink, nothing)
  reader.ptr = ccall(
    (:serd_reader_new, serd),
    Ptr{Cvoid},
    (Cint, Any, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
    syntax, reader, C_NULL, serd_base_sink_ptr, serd_prefix_sink_ptr,
    serd_statement_sink_ptr, serd_end_sink_ptr)
  finalizer(serd_reader_free, reader)
  return reader
end
function serd_base_sink(handle::Ptr{Cvoid}, uri::Ptr{CSerdNode})
  reader = unsafe_pointer_to_objref(handle)::SerdReader
  serd_status(if !isnothing(reader.base_sink)
    reader.base_sink(unsafe_serd_node(uri))
  end)
end
function serd_prefix_sink(handle::Ptr{Cvoid}, name::Ptr{CSerdNode}, uri::Ptr{CSerdNode})
  reader = unsafe_pointer_to_objref(handle)::SerdReader
  serd_status(if !isnothing(reader.prefix_sink)
    reader.prefix_sink(unsafe_serd_node(name), unsafe_serd_node(uri))
  end)
end
function serd_statement_sink(
    handle::Ptr{Cvoid}, flags::Cint, graph::Ptr{CSerdNode},
    subject::Ptr{CSerdNode}, predicate::Ptr{CSerdNode}, object::Ptr{CSerdNode},
    object_datatype::Ptr{CSerdNode}, object_lang::Ptr{CSerdNode}
  )
  reader = unsafe_pointer_to_objref(handle)::SerdReader
  serd_status(if !isnothing(reader.statement_sink)
    reader.statement_sink(SerdStatement(
      flags,
      unsafe_serd_node(graph),
      unsafe_serd_node(subject),
      unsafe_serd_node(predicate),
      unsafe_serd_node(object),
      unsafe_serd_node(object_datatype),
      unsafe_serd_node(object_lang)))
  end)
end
function serd_end_sink(handle::Ptr{Cvoid}, node::Ptr{CSerdNode})
  reader = unsafe_pointer_to_objref(handle)::SerdReader
  serd_status(if !isnothing(reader.end_sink)
    reader.end_sink(unsafe_serd_node(node))
  end)
end

""" Set a function to be called when errors occur during reading.

If no error function is set, errors are printed to stderr in GCC style.
"""
function serd_reader_set_error_sink(reader::SerdReader, error_sink)
  reader.error_sink = error_sink
  serd_error_sink_ptr = isnothing(reader.error_sink) ? C_NULL :
    @cfunction($(CSerd.serd_error_sink), Cint, (Ptr{Cvoid}, Ptr{CSerdError}))
  ccall(
    (:serd_reader_set_error_sink, serd),
    Cvoid,
    (Ptr{Cvoid}, Ptr{Cvoid}, Any),
    reader.ptr, serd_error_sink_ptr, reader)
end
function serd_error_sink(handle::Ptr{Cvoid}, error::Ptr{CSerdError})
  obj = unsafe_pointer_to_objref(handle)::Union{SerdReader,SerdWriter}
  serd_status(if !isnothing(obj.error_sink)
    # FIXME: Include error information besides status code.
    status = unsafe_load(Ptr{Cint}(error))
    obj.error_sink(SerdStatus(status))
  end)
end

""" Enable or disable strict parsing.

The reader is non-strict (lax) by default, which will tolerate URIs with
invalid characters. Setting strict will fail when parsing such files. An error
is printed for invalid input in either case.
"""
function serd_reader_set_strict(reader::SerdReader, strict::Bool)
  ccall((:serd_reader_set_strict, serd), Cvoid, (Ptr{Cvoid}, Cbool),
        reader.ptr, strict)
end

""" Set a prefix to be added to all blank node identifiers.

This is useful when multiple files are to be parsed into the same output
(e.g. a store, or other files). Since Serd preserves blank node IDs, this could
cause conflicts where two non-equivalent blank nodes are merged, resulting in
corrupt data. By setting a unique blank node prefix for each parsed file, this
can be avoided, while preserving blank node names.
"""
function serd_reader_add_blank_prefix(reader::SerdReader, prefix::String)
  ccall((:serd_reader_add_blank_prefix, serd), Cvoid, (Ptr{Cvoid}, Cstring),
        reader.ptr, prefix)
end

""" Set the URI of the default graph.

If this is set, the reader will emit quads with the graph set to the given node
for any statements that are not in a named graph.
"""
function serd_reader_set_default_graph(reader::SerdReader, graph::SerdNode)
  ccall((:serd_reader_set_default_graph, serd), Cvoid, (Ptr{Cvoid}, Ref{CSerdNode}),
        reader.ptr, c_serd_node(graph))
end

""" Read a file at a given URI.
"""
function serd_reader_read_file(reader::SerdReader, uri::String)
  check_serd_status(ccall(
    (:serd_reader_read_file, serd),
    Cint,
    (Ptr{Cvoid}, Cstring),
    reader.ptr, uri
  ))
end

""" Read from UTF8 string.
"""
function serd_reader_read_string(reader::SerdReader, str::String)
  check_serd_status(ccall(
    (:serd_reader_read_string, serd),
    Cint,
    (Ptr{Cvoid}, Cstring),
    reader.ptr, str
  ))
end

""" Free RDF reader.

This function will be called automatically when the Julia Serd reader is
garbage collected.
"""
function serd_reader_free(reader::SerdReader)::Nothing
  if reader.ptr != C_NULL
    ccall((:serd_reader_free, serd), Cvoid, (Ptr{Cvoid},), reader.ptr)
    reader.ptr = C_NULL
  end
  nothing
end

# Writer
########

""" Create a new RDF writer.
"""
function serd_writer_new(syntax::SerdSyntax, style::SerdStyles, io::IO)::SerdWriter
  sink = text -> write(io, text)
  serd_writer_new(syntax, style, sink)
end
function serd_writer_new(syntax::SerdSyntax, style::SerdStyles, sink::Function)::SerdWriter
  serd_sink_ptr = @cfunction(
    $(CSerd.serd_writer_sink), Csize_t, (Ptr{Cvoid}, Csize_t, Ptr{Cvoid}))
  env_ptr = ccall((:serd_env_new, serd), Ptr{Cvoid}, (Ptr{CSerdNode},), C_NULL)
  writer = SerdWriter(C_NULL, env_ptr, sink, nothing)
  writer.ptr = ccall(
    (:serd_writer_new, serd),
    Ptr{Cvoid},
    (Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Any),
    syntax, style, env_ptr, C_NULL, serd_sink_ptr, writer)
  finalizer(serd_writer_free, writer)
  return writer
end
function serd_writer_sink(buf::Ptr{Cvoid}, len::Csize_t, handle::Ptr{Cvoid})
  writer = unsafe_pointer_to_objref(handle)::SerdWriter
  writer.sink(unsafe_string(Ptr{UInt8}(buf), len))
  return len
end

""" Set a function to be called when errors occur during writing.

The error_sink will be called with handle as its first argument. If no error
function is set, errors are printed to stderr.
"""
function serd_writer_set_error_sink(writer::SerdWriter, error_sink)
  writer.error_sink = error_sink
  serd_error_sink_ptr = isnothing(writer.error_sink) ? C_NULL :
    @cfunction($(CSerd.serd_error_sink), Cint, (Ptr{Cvoid}, Ptr{CSerdError}))
  ccall(
    (:serd_writer_set_error_sink, serd),
    Cvoid,
    (Ptr{Cvoid}, Ptr{Cvoid}, Any),
    writer.ptr, serd_error_sink_ptr, writer)
end

""" Set a prefix to be removed from matching blank node identifiers.
"""
function serd_writer_chop_blank_prefix(writer::SerdWriter, prefix::String)
  ccall(
    (:serd_writer_chop_blank_prefix, serd),
    Cvoid,
    (Ptr{Cvoid}, Cstring),
    writer.ptr, prefix)
end

""" Set the current output base URI.
"""
function serd_writer_set_base_uri(writer::SerdWriter, uri::SerdNode)
  check_serd_status(ccall(
    (:serd_writer_set_base_uri, serd),
    Cint,
    (Ptr{Cvoid}, Ref{CSerdNode}),
    writer.ptr, c_serd_node(uri)
  ))
end

""" Set the current root URI.

The root URI should be a prefix of the base URI. The path of the root URI is
the highest path any relative up-reference can refer to.
"""
function serd_writer_set_root_uri(writer::SerdWriter, uri::SerdNode)
  check_serd_status(ccall(
    (:serd_writer_set_root_uri, serd),
    Cint,
    (Ptr{Cvoid}, Ref{CSerdNode}),
    writer.ptr, c_serd_node(uri)
  ))
end

""" Set a namespace prefix.
"""
function serd_writer_set_prefix(writer::SerdWriter, name::SerdNode, uri::SerdNode)
  check_serd_status(ccall(
    (:serd_writer_set_prefix, serd),
    Cint,
    (Ptr{Cvoid}, Ref{CSerdNode}, Ref{CSerdNode}),
    writer.ptr, c_serd_node(name), c_serd_node(uri)
  ))
end

""" Write a statement (RDF triple or quad).
"""
function serd_writer_write_statement(writer::SerdWriter, stmt::SerdStatement)
  check_serd_status(ccall(
    (:serd_writer_write_statement, serd),
    Cint,
    (Ptr{Cvoid}, Cint, Ptr{CSerdNode}, Ref{CSerdNode}, Ref{CSerdNode},
     Ref{CSerdNode}, Ptr{CSerdNode}, Ptr{CSerdNode}),
    writer.ptr,
    stmt.flags,
    c_maybe_serd_node(stmt.graph),
    c_serd_node(stmt.subject),
    c_serd_node(stmt.predicate),
    c_serd_node(stmt.object),
    c_maybe_serd_node(stmt.object_datatype),
    c_maybe_serd_node(stmt.object_lang)
  ))
end

""" Mark the end of an anonymous node's description.
"""
function serd_writer_end_anon(writer::SerdWriter, node::SerdNode)
  check_serd_status(ccall(
    (:serd_writer_end_anon, serd),
    Cint,
    (Ptr{Cvoid}, Ref{CSerdNode}),
    writer.ptr, c_serd_node(node)
  ))
end

""" Finish a write.
"""
function serd_writer_finish(writer::SerdWriter)
  ccall((:serd_writer_finish, serd), Cvoid, (Ptr{Cvoid},), writer.ptr)
end
Base.close(writer::SerdWriter) = serd_writer_finish(writer)

""" Free RDF writer.

This function will be called automatically when the Julia Serd writer is
garbage collected.
"""
function serd_writer_free(writer::SerdWriter)::Nothing
  if writer.ptr != C_NULL
    ccall((:serd_writer_free, serd), Cvoid, (Ptr{Cvoid},), writer.ptr)
    writer.ptr = C_NULL
  end
  if writer.env != C_NULL
    ccall((:serd_env_free, serd), Cvoid, (Ptr{Cvoid},), writer.env)
    writer.env = C_NULL
  end
  nothing
end

end


text:
function normal_order(x; rewriter = serial_normal_order_simplifier)
    !istree(x) && return x
    return simplify(x, rewriter = rewriter)
end

needs_hilbert_order_sorting(f) = x -> isterm_op(f)(x) && !issorted(arguments(x), lt = <ₑ)

function sort_hilbert_args(f, t)
    args = arguments(t)
    if length(args) < 2
        return similarterm(t, f, args)
    elseif length(args) == 2
        x, y = args
        return similarterm(t, f, x <ₑ y ? [x,y] : [y,x])
    end
    args = args isa Tuple ? [args...] : args
    similarterm(t, f, sort(args, lt=<ₑ))
end

<ₑ(a::Real,    b::Real) = abs(a) < abs(b)
<ₑ(a::Complex, b::Complex) = (abs(real(a)), abs(imag(a))) < (abs(real(b)), abs(imag(b)))
<ₑ(a::Real,    b::Complex) = true
<ₑ(a::Complex, b::Real) = false

<ₑ(a::Symbolic, b::Number) = false
<ₑ(a::Number,   b::Symbolic) = true
<ₑ(a::Sym{Number},   b::Sym{Number}) = a.name < b.name
<ₑ(a::Sym{Number},   b::Sym{QOperator}) = true
<ₑ(a::Sym{QOperator}, b::Sym{Number}) = false
<ₑ(a::Sym{QOperator}, b::Sym{QOperator}) = getmetadata(a, QOperatorMeta).h_idx < getmetadata(b, QOperatorMeta).h_idx

arglength(a) = length(arguments(a))
function <ₑ(a, b)
    # isequal(a, b) && return false
    if istree(a) && (b isa Symbolic && !istree(b))
        if symtype(b) != QOperator
            return false
        else
            if symtype(a) != QOperator
                return true
            elseif _isthere_h_idx(a, b)
                return false ## Needs to return false, following that for any real number r < r = false
            else
                return getmetadata(b, QOperatorMeta).h_idx > _min_h_idx(a)
            end
        end
    elseif istree(b) && (a isa Symbolic && !istree(a))
        if symtype(a) != QOperator
            return true
        else
            if symtype(b) != QOperator
                return false
            elseif _isthere_h_idx(b, a)
                return false
            else
                return getmetadata(a, QOperatorMeta).h_idx < _min_h_idx(b)
            end
        end
    elseif istree(a) && istree(b)
        if symtype(b) != QOperator
            if symtype(a) != QOperator
                return cmp_term_term(a,b)
            else
                return false
            end
        else
            return false ## Needs to return false, following that for any real number r < r = false
        end
    else
        return a <ₑ b
    end
end

function <ₑ(a::Symbol, b::Symbol)
    # Enforce the order [+,-,\,/,^,*]
    if b === :*
        a in (:^, :/, :\, :-, :+)
    elseif b === :^
        a in (:/, :\, :-, :+) && return true
    elseif b === :/
        a in (:\, :-, :+) && return true
    elseif b === :\
        a in (:-, :+) && return true
    elseif b === :-
        a === :+ && return true
    elseif a in (:*, :^, :/, :-, :+)
        false
    else
        a < b
    end
end

<ₑ(a::Function, b::Function) = nameof(a) <ₑ nameof(b)

<ₑ(a::Type, b::Type) = nameof(a) <ₑ nameof(b)

function cmp_term_term(a, b)
    la = arglength(a)
    lb = arglength(b)

    if la == 0 && lb == 0
        return operation(a) <ₑ operation(b)
    elseif la === 0
        return operation(a) <ₑ b
    elseif lb === 0
        return a <ₑ operation(b)
    end

    na = operation(a)
    nb = operation(b)

    if 0 < arglength(a) <= 2 && 0 < arglength(b) <= 2
        # e.g. a < sin(a) < b ^ 2 < b
        @goto compare_args
    end

    if na !== nb
        return na <ₑ nb
    elseif arglength(a) != arglength(b)
        return arglength(a) < arglength(b)
    else
        @label compare_args
        aa, ab = arguments(a), arguments(b)
        if length(aa) !== length(ab)
            return length(aa) < length(ab)
        else
            terms = zip(Iterators.filter(!is_literal_number, aa), Iterators.filter(!is_literal_number, ab))

            for (x,y) in terms
                if x <ₑ y
                    return true
                elseif y <ₑ x
                    return false
                end
            end

            # compare the numbers
            nums = zip(Iterators.filter(is_literal_number, aa),
                       Iterators.filter(is_literal_number, ab))

            for (x,y) in nums
                if x <ₑ y
                    return true
                elseif y <ₑ x
                    return false
                end
            end

        end
        return na <ₑ nb # all args are equal, compare the name
    end
end

function flatten_term(⋆, x)
    args = arguments(x)
    # flatten nested ⋆
    flattened_args = []
    for t in args
        if istree(t) && operation(t) === (⋆)
            if isnotflat(⋆)(t)
                append!(flattened_args, arguments(flatten_term(⋆, t)))
            else
                append!(flattened_args, arguments(t))
            end
        else
            push!(flattened_args, t)
        end
    end
    similarterm(x, ⋆, flattened_args)
end

function isnotflat(⋆)
    function (x)
        if istree(x)
            if operation(x) === (⋆)
                args = arguments(x)
                for t in args
                    if istree(t) && operation(t) === (⋆)
                        return true
                    end
                end
            end
        end
        return false
    end
end

function _isthere_h_idx(my_term, my_op)
    metadt = getmetadata(my_op, QOperatorMeta)
    h_ind = metadt.h_idx
    op_type = metadt.type
    if symtype(my_term) != QOperator
        return false
    else
        if !istree(my_term)
            return getmetadata(my_term, QOperatorMeta).h_idx == h_ind# && getmetadata(my_term, QOperatorMeta).type == op_type
        else
            args = arguments(my_term)
            for arg in args
                _isthere_h_idx(arg, my_op) && return true
            end
        end
    end
    return false
end

function _min_h_idx(my_term::Term{QOperator})
    idxs = _list_h_idxs(my_term)
    return minimum(idxs)
end

function _list_h_idxs(my_term::Term{QOperator})
    args = arguments(my_term)
    idxs = []
    for arg in args
        if istree(arg)
            if symtype(arg) == QOperator
                h_idxs = _list_h_idxs(arg)
            else
                h_idxs = []
            end
        else
            if symtype(arg) == QOperator
                h_idxs = [getmetadata(arg, QOperatorMeta).h_idx]
            else
                h_idxs = []
            end
        end
        for h_idx in h_idxs
            if !(h_idx in idxs)
                push!(idxs, h_idx)
            end
        end
    end
    return sort(idxs)
end

get_qop_order(x::Symbolic{QOperator}) = 1
get_qop_order(x::Symbolic{Number}) = 0
get_qop_order(x::Number) = 0

function get_qop_order(x::Term{QOperator})
    args = arguments(x)
    if isterm_op(*)(x)
        return +([get_qop_order(arg) for arg in args]...)
    elseif isterm_op(+)(x)
        return maximum([get_qop_order(arg) for arg in args])
    elseif isterm_op(/)(x)
        return get_qop_order(args[1]) - get_qop_order(args[2])
    elseif isterm_op(^)(x)
        if args[2] isa Number
            return get_qop_order(args[1]) * args[2]
        else
            return 0
        end
    else
        return 0
    end
end

function needs_sum_ordering(x)
    !isterm_op(+)(x) && return false
    args = arguments(x)
    orders = [get_qop_order(arg) for arg in args]
    return !issorted(orders)
end

function order_sum(x)
    args = arguments(x)
    orders = [get_qop_order(arg) for arg in args]
    return +(arguments(x)[sortperm(orders)]...)
end

ext: jl


text:
import DomainSets: Interval, Ball, infimum, supremum

@deprecate IntervalDomain(a, b) Interval(a, b)
@deprecate CircleDomain() Ball()

# type piracy on Interval for downstream compatibility to be reverted once upgrade is complete
function Base.getproperty(domain::Interval, sym::Symbol)
    if sym === :lower
        @warn "domain.lower is deprecated, use infimum(domain) instead"
        return infimum(domain)
    elseif sym === :upper
        @warn "domain.upper is deprecated, use supremum(domain) instead"
        return supremum(domain)
    else
        return getfield(domain, sym)
    end
end


text:
module Yota

export grad, update!, @primitive, @grad

include("core.jl")

version() = v"0.1.0-2"

end


text:
module TestForeach

using Test
using ThreadsX
using Referenceables

simdflags = [false, true, :ivdep]

function addavg!(c, a, b)
    c[] += (a + b) / 2
end

@testset for simd in [simdflags; Val.(simdflags)], n in [10, 100, 1000]
    basesize = n^2 ÷ 4

    @testset "linear" begin
        A = randn(n, n)
        B = randn(n, n)
        C0 = randn(n, n)
        C1 = copy(C0)
        foreach(addavg!, referenceable(C0), A, B)
        ThreadsX.foreach(addavg!, referenceable(C1), A, B; simd = simd, basesize = basesize)
        @test C0 == C1
    end

    @testset "cartesian" begin
        A = randn(n, n)
        C0 = randn(n, n)
        C1 = copy(C0)
        foreach(addavg!, referenceable(C0), A, A')
        ThreadsX.foreach(
            addavg!,
            referenceable(C1),
            A,
            A';
            simd = simd,
            basesize = basesize,
        )
        @test C0 == C1
    end

    @testset "foreach(f, referenceable(C), product(A, B))" begin
        A = 1:3
        B = 1:2
        C = zeros(3, 2)

        ThreadsX.foreach(referenceable(C), Iterators.product(A, B)) do c, (a, b)
            c[] = a * b
        end

        @test C == A .* reshape(B, 1, :)
    end
end

@testset "argument validation" begin
    @test_throws(
        ArgumentError("Invalid `simd` option: :invalid"),
        ThreadsX.foreach(identity, 1:0; simd = :invalid),
    )
    @test_throws(Exception, ThreadsX.foreach(identity, 1:0; simd = Val(:invalid)),)
end

end  # module


text:
using Highlights

using Test
using InteractiveUtils

#
# Utilities.
#

function tokentest(lexer, source, expects...)
    tokens = Highlights.Compiler.lex(source, lexer).tokens
    @test length(tokens) == length(expects)
    @test join([s for (n, s) in expects]) == source
    for (token, (name, str)) in zip(tokens, expects)
        @test token.value == name
        @test SubString(source, token.first, token.last) == str
    end
end

function print_all(lexer, file)
    source = read(joinpath(@__DIR__, "samples", file), String)
    buffer = IOBuffer()
    for m in ["html", "latex"]
        mime = MIME("text/$m")
        for theme in subtypes(Themes.AbstractTheme)
            stylesheet(buffer, mime, theme)
            highlight(buffer, mime, source, lexer, theme)
        end
    end
    return buffer
end

#
# Setup.
#

using Highlights.Tokens, Highlights.Themes, Highlights.Lexers

# Error reporting for broken lexers.
abstract type BrokenLexer <: Highlights.AbstractLexer end
@lexer BrokenLexer Dict(
    :tokens => Dict(:root => [(r"\w+", TEXT, (:b, :a))], :a => [], :b => []),
)

# Lexer inheritance.
abstract type ParentLexer <: Highlights.AbstractLexer end
abstract type ChildLexer <: ParentLexer end
@lexer ParentLexer Dict(
    :tokens => Dict(:root => [(r"#.*$"m, COMMENT)])
)
@lexer ChildLexer Dict(
    :tokens => Dict(:root => [:__inherit__, (r"\d+", NUMBER), (r".", TEXT)])
)

abstract type NumberLexer <: Highlights.AbstractLexer end
@lexer NumberLexer Dict(
    :tokens => Dict(
        :root => [
            (r"0b[0-1]+", NUMBER_BIN),
            (r"0o[0-7]+", NUMBER_OCT),
            (r"0x[0-9a-f]+", NUMBER_HEX),
        ],
        :integers => [
            (r"\d+", NUMBER_INTEGER),
        ],
    ),
)

abstract type SelfLexer <: Highlights.AbstractLexer end
@lexer SelfLexer Dict(
    :tokens => Dict(
        :root => [
            (r"(#)( )(.+)(;)$"m, (COMMENT, WHITESPACE, :root, PUNCTUATION)),
            (r"(!)( )(.+)(;)$"m, (COMMENT, WHITESPACE, :string, PUNCTUATION)),
            (r"\d+", NUMBER),
            (r"\w+", NAME),
            (r" ", WHITESPACE),
        ],
        :string => [
            (r"'[^']*'", STRING),
            (r"0[box][\da-f]+", NumberLexer),
            (r"\d+", NumberLexer => :integers),
            (r"\s", WHITESPACE),
        ],
    ),
)

# Custom output format. A DOM-like node "buffer".

struct Node
    name::Symbol
    text::String
    class::String
end

struct NodeBuffer <: IO
    nodes::Vector{Node}
end

function Format.render(io::IO, ::MIME"ast/dom", tokens::Format.TokenIterator)
    for (str, id, style) in tokens
        push!(io.nodes, Node(:span, str, "hljl-$id"))
    end
end

#
# Testsets.
#

@testset "Highlights" begin
    @testset "Lexers" begin
        @testset for file in readdir(joinpath(@__DIR__, "lexers"))
            include("lexers/$file")
        end
        @testset "Self-referencing" begin
            tokentest(
                SelfLexer,
                "# word;",
                COMMENT => "#",
                WHITESPACE => " ",
                NAME => "word",
                PUNCTUATION => ";",
            )
            tokentest(
                SelfLexer,
                "# 1 word;",
                COMMENT => "#",
                WHITESPACE => " ",
                NUMBER => "1",
                WHITESPACE => " ",
                NAME => "word",
                PUNCTUATION => ";",
            )
            tokentest(
                SelfLexer,
                "! '...';",
                COMMENT => "!",
                WHITESPACE => " ",
                STRING => "'...'",
                PUNCTUATION => ";",
            )
            tokentest(
                SelfLexer,
                "! 0b1010;",
                COMMENT => "!",
                WHITESPACE => " ",
                NUMBER_BIN => "0b1010",
                PUNCTUATION => ";",
            )
            tokentest(
                SelfLexer,
                "! 0b1010 0xacd1;",
                COMMENT => "!",
                WHITESPACE => " ",
                NUMBER_BIN => "0b1010",
                WHITESPACE => " ",
                NUMBER_HEX => "0xacd1",
                PUNCTUATION => ";",
            )
            tokentest(
                SelfLexer,
                "! 1234 0b01;",
                COMMENT => "!",
                WHITESPACE => " ",
                NUMBER_INTEGER => "1234",
                WHITESPACE => " ",
                NUMBER_BIN => "0b01",
                PUNCTUATION => ";",
            )
        end
        @testset "Inheritance" begin
            tokentest(
                ParentLexer,
                "# ...",
                COMMENT => "# ...",
            )
            tokentest(
                ParentLexer,
                "1 # ...",
                ERROR => "1",
                ERROR => " ",
                COMMENT => "# ...",
            )
            tokentest(
                ChildLexer,
                "1 # ...",
                NUMBER => "1",
                TEXT => " ",
                COMMENT => "# ...",
            )
        end
        @testset "Utilities" begin
            let w = Lexers.words(["if", "else"]; prefix = "\\b", suffix = "\\b")
                @test occursin(w, "if")
                @test !occursin(w, "for")
                @test !occursin(w, "ifelse")
            end
            let c = Highlights.Compiler.Context("@lexer CustomLexer dict(")
                @test Lexers.julia_is_macro_identifier(c) == 1:6
                @test Lexers.julia_is_iterp_identifier(c) == 0:0
            end
            let c = Highlights.Compiler.Context("raw\"\"\"...\"\"\"")
                @test Lexers.julia_is_triple_string_macro(c) == 1:6
            end
        end
        @testset "Errors" begin
            tokentest(BrokenLexer, " ", ERROR => " ")
        end
    end
    @testset "Themes" begin
        let s = Themes.Style("fg: 111")
            @test s.fg.r == 0x11
            @test s.fg.g == 0x11
            @test s.fg.b == 0x11
            @test s.fg.active == true

            @test s.bg.r == 0x00
            @test s.bg.g == 0x00
            @test s.bg.b == 0x00
            @test s.bg.active == false

            @test !s.bold
            @test !s.italic
            @test !s.underline
        end
        let s = Themes.Style("bg: f8c; italic; bold")
            @test s.fg.r == 0x00
            @test s.fg.g == 0x00
            @test s.fg.b == 0x00
            @test s.fg.active == false

            @test s.bg.r == 0xff
            @test s.bg.g == 0x88
            @test s.bg.b == 0xcc
            @test s.bg.active == true
            @test s.bold
            @test s.italic
            @test !s.underline
        end
        let t = Themes.maketheme(Themes.DefaultTheme),
            m = Themes.metadata(Themes.DefaultTheme)
            @test t.base == m[:style]
            @test t.styles[1] == m[:tokens][TEXT]
        end
    end
    @testset "Format" begin
        render = function(mime, style)
            buffer = IOBuffer()
            Highlights.Format.render(buffer, mime, style)
            return Highlights.takebuf_str(buffer)
        end
        @testset "CSS" begin
            let mime = MIME("text/css")
                @test render(mime, Themes.Style("fg: 111"))   == "color: rgb(17,17,17); "
                @test render(mime, Themes.Style("bg: 111"))   == "background-color: rgb(17,17,17); "
                @test render(mime, Themes.Style("bold"))      == "font-weight: bold; "
                @test render(mime, Themes.Style("italic"))    == "font-style: italic; "
                @test render(mime, Themes.Style("underline")) == "text-decoration: underline; "
            end
        end
        @testset "LaTeX" begin
            escape = function(mime,str;charescape=false)
                buffer = IOBuffer()
                Highlights.Format.escape(buffer,mime,str,charescape=charescape)
                return Highlights.takebuf_str(buffer)
            end
            render_nonwhitespace = function(mime,str)
                buffer = IOBuffer()
                Highlights.Format.render_nonwhitespace(buffer,mime,str,0,nothing)
                return Highlights.takebuf_str(buffer)
            end
            escapeinside(str) = "(*@{"*str*"}@*)"
            let mime = MIME("text/latex")
                @test render(mime, Themes.Style("fg: 111"))   == "[1]{\\textcolor[RGB]{17,17,17}{#1}}"
                @test render(mime, Themes.Style("bg: 111"))   == "[1]{\\colorbox[RGB]{17,17,17}{#1}}"
                @test render(mime, Themes.Style("bold"))      == "[1]{\\textbf{#1}}"
                @test render(mime, Themes.Style("italic"))    == "[1]{\\textit{#1}}"
                @test render(mime, Themes.Style("underline")) == "[1]{\\underline{#1}}"
                mime = MIME("text/latex")
                ebrace_L = escape(mime,"{")
                ebrace_R = escape(mime,"}")
                # This is placed in a verbatim env. so there must be an additional charescape
                @test escape(mime,"{}",charescape=true) == escapeinside(ebrace_L)*escapeinside(ebrace_R)
                # This is already escaped so now charescape
                @test escape(mime,"{}",charescape=false) == ebrace_L*ebrace_R
                # This is placed in a verbatim so whitespace is already handled by the verbatim env.
                # Only special characters need to be escaped
                @test escape(mime,"   some text\n next line{}",charescape=true) ==
                    "   some text\n next line"*escapeinside(ebrace_L)*escapeinside(ebrace_R)
                # This is escaped so LaTeX automatically removes whitespace characters.
                # Also no additional escapeinside for special characters needed.
                s = escape(mime," ",charescape=false)
                @test escape(mime,"   some text\n next line{}",charescape=false) ==
                    "$(s)$(s)$(s)some$(s)text\n$(s)next$(s)line"*ebrace_L*ebrace_R
                r(x) = "(*@\\HLJL0{$x}@*)"
                @test render_nonwhitespace(MIME("text/latex"),"\tfoo  \nbar\t") == "\t"*r("foo")*"  \n"*r("bar")*"\t"
                @test render_nonwhitespace(MIME("text/latex"),"foo") == r("foo")
                @test render_nonwhitespace(MIME("text/latex"),"\t") == "\t"
                @test render_nonwhitespace(MIME("text/latex"),"α") == r("α")
                @test render_nonwhitespace(MIME("text/latex")," α\tβfoo ") == " "*r("α")*"\t"*r("βfoo")*" "
                @test render_nonwhitespace(MIME("text/latex")," α ") == " "*r("α")*" "
            end
        end
        @testset "Custom Nodes" begin
            let buf = NodeBuffer([]),
                eq = (a, b) -> a.name == b.name && a.text == b.text && a.class == b.class
                highlight(buf, MIME("ast/dom"), "x = 1", Lexers.JuliaLexer)
                @test eq(buf.nodes[1], Node(:span, "x", "hljl-n"))
                @test eq(buf.nodes[2], Node(:span, " ", "hljl-t"))
                @test eq(buf.nodes[3], Node(:span, "=", "hljl-oB"))
                @test eq(buf.nodes[4], Node(:span, " ", "hljl-t"))
                @test eq(buf.nodes[5], Node(:span, "1", "hljl-ni"))
            end
        end
    end
    @testset "Compiler" begin
        let buf = IOBuffer()
            Highlights.Compiler.debug(buf, "x", Highlights.Lexers.JuliaLexer)
            @test Highlights.takebuf_str(buf) == "<NAME> := \"x\"\n"
            # Should print nothing to STDOUT.
            Highlights.Compiler.debug("", Highlights.Lexers.JuliaLexer)
        end
    end
    @testset "Tokens" begin
        @test Highlights.Tokens.TokenValue(:TEXT).value == 1
        @test Highlights.Tokens.parent(:TEXT) == :TEXT
        @test Highlights.Tokens.parent(:COMMENT) == :TEXT
        @test Highlights.Tokens.parent(:COMMENT_SINGLE) == :COMMENT
    end
    @testset "Miscellaneous" begin
        @test Highlights.lexer("julia") == Lexers.JuliaLexer
        @test Highlights.lexer("jl") == Lexers.JuliaLexer
        @test_throws ArgumentError Highlights.lexer("???")
    end
end

# Hack to make the test in External actually usable...
using Pkg
Pkg.activate(joinpath(@__DIR__, "External"))
Pkg.develop("Highlights")
using External


text:
export Transform

struct Transform{ℱ, ℬ}
    forward::ℱ
    backward::ℬ
end

function Transform(𝒢::FourierGrid; arraytype=Array)
    grid_size = length.(𝒢.grid)
    f = arraytype(randn(grid_size...) .+ 0im)
    FFTW.set_num_threads(Threads.nthreads())
    P = plan_fft(f) # typeof(iP) <: AbstractFFTs.ScaledPlan
    iP = plan_ifft(f)
    return Transform(P, iP)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
# commenting out till we decide to start deprecating things again

# import Base: @deprecate

# TODO: remove these definitions in year 2021
by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))

aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "ComputabilityTheory.jl"))

using .ComputabilityTheory
using Test

@time @testset "ComputabilityTheory.jl" begin
### CODING
	@test 1 ∸ 2 == 0
	@test 2 ∸ 1 == 1

    @test pair_tuple(5,7) == 83
    @test pair_tuple(5,7,20) == 5439
    @test pair_tuple([1,2,3,4,5,6,7,8,9]...) == 131504586847961235687181874578063117114329409897615188504091716162522225834932122128288032336298131

    @test π(83, 2, algebraic = false) == (5,7)
    @test π(83, 2, 0, algebraic = false) == 5
    @test π(83, 3, algebraic = false) == (2, 0, 7)
    @test π(1023, 2, 1, algebraic = false) == 11
    @test π(big(1315045868479612356871818745780631171143), 1, algebraic = false) == 1315045868479612356871818745780631171143
    @test π(5987349857934, 0, algebraic = false) == nothing
    @test π(83, 3) == (2, 0, 7)
    @test π(10001, 10) == (0, 0, 0, 0, 0, 0, 1, 3, 4, 9)
    @test π(big(1315045868479612356871818745780631171143), 1) == 1315045868479612356871818745780631171143
    @test π(5987349857934, 0) == nothing
    small_random = abs(rand(1:100))
    large_random = abs(rand(1:10000))
    @test π(small_random, 3) == π(small_random, 3)
    @test π(large_random, 2) == π(large_random, 2)
    @test π(972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535) == (7, 44097457325828774284791367377726710860393112311607949541420228252579855118277296197320694764257474575167892495444390648066046785424108252105160)
    # @test π(rand(1:10^1000), algebraic)

    @test cℤ([0,-1,1,-2,2,-3,3,-4,4]) == [0, 1, 2, 3, 4, 5, 6, 7, 8]
    @test cℤ((-1,2)) == 16
    @test cℤ((1,1)) == 12
    @test cℤ(3,4) == (6, 8)

    @test cℤ⁻¹([0, 1, 2, 3, 4, 5, 6, 7, 8]) == [0, -1, 1, -2, 2, -3, 3, -4, 4]
    @test cℤ(cℤ⁻¹(79)) == 79
    @test cℤ⁻¹(cℤ(-40)) == -40
    @test cℤ⁻¹(10029) == -5015
	int_random1 = rand(Int)
	int_random2 = rand(Int)
	@test cℤ(cℤ⁻¹(abs(int_random1))) == abs(int_random1)
	@test cℤ⁻¹(cℤ(int_random2)) == int_random2

### PROGRAMME
	@test GoToProgramme(121).programme_length == 1
	@test GoToProgramme(972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535).programme_length == 7

	io = IOBuffer()
	show_programme(io, 121)
	@test occursin(r"0[[:space:]]*halt", String(take!(io)))

	io = IOBuffer()
	show_programme(io, 972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535)
	@test occursin(r"6[[:space:]]*halt", String(take!(io)))

### INSTRUCTIONS
	@test Instruction(7).instruction == (1, 2)
	@test Instruction(7).first == 1
	@test Instruction(7).second == 2
	@test Instruction(7).third == nothing
	@test Instruction(58).second == 1
	@test Instruction(58).third == 2
	@test Instruction((3, (1, 2))).I == 58
	@test Instruction((3, 1, 2)).I == 58
	@test Instruction(3, 1, 2).I == 58
	@test Instruction(4, 0).I == 14

	@test run_goto_programme(121) == tuple(0)
	@test run_goto_programme(363183787614755732766753446033240, RegisterMachine(1, 2, 3)) == tuple(2, 1, 3)
	@test run_goto_programme(972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535, (3, 1, 1)) == tuple(4, 0, 1)
	@test run_goto_programme(223927124717252638681796127568111496742698512180714163023526568831822132049774230195118493600338381851704911289089017787857413597404988541338847992210485499035130781598573348126975813033599979584118468386578559370872281223796019975651453204160222340161200928955988080303016030336930443304245145601448154569733221264381538305168568931449480758528673246332402357557389885850502823987502500615923705666157574953869234044687656342114443578201036399996781757409020565065471713601342262949677684562201881239545438799475865994538241042725394500589505167297761345293434717297622500387774051925726286270828011593028154888482643933182720230237091007716467114407940246910876104242298811097230924566373379040422996695018496814405564315525869126758556165426204979060785448299907433924693776629082547220185023385651497716000521725878783901307004251127081977484444374187001323436144140037580903329365965139055311902469340438083050750301104320244848087790925135435440549816132311216931423038414958594054194345955328449597351606723033103611762253425583323756613052352836016212463541320408142893340287911691519461708609068873926465524275772411714118028019844718904105500330674632838771428681137130136100546871186789563992080600787879878511406660580151025959186841418554827311891320184308989618516409485943704466045877998837588192158884659614365477933514118459878126838662936672543012162679280965785850875755241288227892173373302073145637378088093219495972690064022128760997033867243024673762202229065063901594672092113076799336225728489852391219624477091216808040617586878843148705687256752473769380170895833165256985381066579262059980097322456692517960764426357076378502591620840598428480692392558582556549972917361099374457224379447165159150577431975140314355645616431652416541261437750090494682945620970384916782903782096424825827978169475334019068642656156261072845169384253772103615589882161531092319774765872915674167557806826727743079891261233740545808889798017329652560892724042630851029267797497619411669453984450602449597986084188806823478595595515128790674171852156626442453849830854084257473658507273475076486871236759558604076795418209173234381881185600648269038607854092549797177142124163026021172162187852101251819011578114092027496262838478199601979508935066069585104555431366780222617619635159394639488403252338359673576720653539411491577093278404645966204761032796579698188013739359885785188427988372869573092916550596059755937168888105513627640203861038587598828520438054944864945385514172661237373005642977222691283341092758484095459330166556676864303459850601868153088510594207374543818165331580579808638399043608636650159746931808355843721671440780553546753959760350975906249580764849993052921544867678610217996851589828713473915270585295877597236593571404899502354444405479601753716517182781778676830983428681476532848100787997408528892345104613523387370533475767447648602910646760355660958540505065155007961785878896275786367616723562594625648464444056163755337577182731373365691826878791057479596578563844746371661) == tuple(1, 0)

	# Basically all we can test with rand is that the output has a halting instruction in the correct place
	io = IOBuffer()
	show_programme(io, rand(GoToProgramme, 3))
	@test occursin(r"[[:digit:]][[:space:]]*halt", String(take!(io)))
	io = IOBuffer()
	show_programme(io, rand(GoToProgramme, 11))
	@test occursin(r"[[:digit:]][[:digit:]][[:space:]]*halt", String(take!(io)))
end # end testset


text:
###############################################################################
#
# Calculate diffusive thermo and SIF
#
###############################################################################
"""
    diffusive_S(τ_dd::Array{FT},
                ρ_dd::Array{FT},
                S⁻::Array{FT},
                S⁺::Array{FT},
                boundary_top::Array{FT},
                boundary_bottom::Array{FT},
                rsoil::Array{FT}
    ) where {FT<:AbstractFloat}

Computes 2-stream diffusive radiation transport (used for thermal and SIF) given:
- `τ_dd` A 2D Array with layer reflectances
- `ρ_dd` A 2D Array with layer transmissions
- `S⁻` A 2D Array with layer source terms in the downwelling direction
- `S⁺` A 2D Array with layer source terms in the upwelling direction
- `boundary_top` A 1D array with downwelling radiation at the top (top of canopy)
- `boundary_bottom` A 1D array with upwnwelling radiation at the bottom (soil)
- `rsoil` A 1D array with soil reflectance
"""
function diffusive_S(
            τ_dd::Array{FT},
            ρ_dd::Array{FT},
            S⁻::Array{FT},
            S⁺::Array{FT},
            boundary_top::Array{FT},
            boundary_bottom::Array{FT},
            rsoil::Array{FT}
) where {FT<:AbstractFloat}
    # Get dimensions (1st is wavelength, 2nd is layers), for Stefab Boltzmann
    # just one effective wavelength
    nWL,nl = size(τ_dd);
    Xdd    = similar(τ_dd);
    Rdd    = zeros(FT,nWL,nl+1)
    # Create Y,U matrices (see mSCOPE paper, eqs 30-39)
    Y  = similar(τ_dd);
    U  = similar(Rdd);
    E⁻ = similar(Rdd);
    E⁺ = similar(Rdd);
    net_diffuse = similar(τ_dd);

    Rdd[:,end] = rsoil

    # Source at TOC (0 for SIF, downwelling LW from atmosphere for thermal)
    E⁻[:,1]  = boundary_top;
    # Source at bottom layer (0 for SIF, emitted from soil for thermal)
    U[:,end].= boundary_bottom;

    @inbounds for j=nl:-1:1   # bottom to top
        dnorm    = 1 .-ρ_dd[:,j].*Rdd[:,j+1];
        Xdd[:,j] = τ_dd[:,j]./dnorm;
        Rdd[:,j] = ρ_dd[:,j]+τ_dd[:,j].*Rdd[:,j+1].*Xdd[:,j];
        Y[:,j]   = (ρ_dd[:,j].*U[:,j+1]+S⁻[:,j])./dnorm;
        U[:,j]   = τ_dd[:,j].*(Rdd[:,j+1].*Y[:,j]+U[:,j+1])+S⁺[:,j];
    end
    @inbounds for j=1:nl      # from top to bottom
        E⁻[:,j+1]= Xdd[:,j].*E⁻[:,j]+Y[:,j];
        E⁺[:,j]  = Rdd[:,j].*E⁻[:,j]+U[:,j];
    end
    E⁺[:,end] = Rdd[:,end].*E⁻[:,end]+U[:,end];
    # compute net diffuse radiation per layer:
    @inbounds for j=1:nl
        net_diffuse[:,j]= (E⁻[:,j] + E⁺[:,j+1]).*(1 .-(τ_dd[:,j]+ρ_dd[:,j]))
    end

    return E⁻,E⁺,net_diffuse
end

"""
    diffusive_S!(
                sf_con::SFCache{FT},
                soil::SoilOpticals{FT},
                rt_dim::RTDimensions
    ) where {FT<:AbstractFloat}

Computes 2-stream diffusive radiation transport (used for thermal and SIF),
    given
- `sf_con` [`SFCache`](@ref) type cache
- `soil` [`SoilOpticals`](@ref) type struct
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function diffusive_S!(
            sf_con::SFCache{FT},
            soil::SoilOpticals{FT},
            rt_dim::RTDimensions
) where {FT<:AbstractFloat}
    # 1. unpack values from sf_con
    @unpack dnorm,  F⁻, F⁺, net_diffuse, Rdd, S⁻, S⁺, U, Xdd, Y, zeroB, ρ_dd,
            τ_dd = sf_con;
    @unpack albedo_SW_SIF = soil;
    @unpack nLayer, nLevel = rt_dim;

    # Get dimensions (1st is wavelength, 2nd is layers), for Stefan Boltzmann
    # just one effective wavelength
    dnorm .= 1 .- view(ρ_dd, :, 1);

    Rdd[:,end] .= albedo_SW_SIF;

    # Source at TOC (0 for SIF, downwelling LW from atmosphere for thermal)
    F⁻[:,1]  .= zeroB;
    # Source at bottom layer (0 for SIF, emitted from soil for thermal)
    U[:,end] .= zeroB;

    @inbounds for j=nLayer:-1:1   # bottom to top
        dnorm    .= 1 .- view(ρ_dd, :, j) .* view(Rdd, :, j+1);
        Xdd[:,j] .= view(τ_dd, :, j) ./ dnorm;
        Rdd[:,j] .= view(ρ_dd, :, j) .+ view(τ_dd, :, j) .* view(Rdd, :, j+1) .* view(Xdd, :, j);
        Y[:,j]   .= (view(ρ_dd, :, j) .* view(U, :, j+1) .+ view(S⁻, :, j)) ./ dnorm;
        U[:,j]   .= view(τ_dd, :, j) .* (view(Rdd, :, j+1) .* view(Y, :, j) .+ view(U, :, j+1)) .+ view(S⁺, :, j);
    end
    @inbounds for j=1:nLayer      # from top to bottom
        F⁻[:,j+1] .= view(Xdd, :, j) .* view(F⁻, :, j) .+ view(Y, :, j);
        F⁺[:,j  ] .= view(Rdd, :, j) .* view(F⁻, :, j) .+ view(U, :, j);
    end
    F⁺[:,end] .= view(Rdd, :, nLevel) .* view(F⁻, :, nLevel) .+ view(U, :, nLevel);

    # compute net diffuse radiation per layer:
    @inbounds for j=1:nLayer
        net_diffuse[:,j] .= (view(F⁻, :, j) .+ view(F⁺, :, j+1)) .* (1 .- (view(τ_dd, :, j) .+ view(ρ_dd, :, j)));
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 8a8be09a-49e9-4563-9afc-ae72bd05d26d
using BenchmarkTools, Colors, Manopt, Manifolds, PlutoUI, Random

# ╔═╡ f58f28fe-c3d8-11ec-0acb-e948605dd63f
md"""
# Stochastic Gradient Descent

This tutorial illustrates how to use the [`stochastic_gradient_descent`](https://manoptjl.org/stable/solvers/stochastic_gradient_descent.html)
solver and different [`DirectionUpdateRule`](https://manoptjl.org/stable/solvers/gradient_descent.html#Direction-Update-Rules-1)s in order to introduce
the average or momentum variant, see [Stochastic Gradient Descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent).

Computationally we look at a very simple but large scale problem,
the Riemannian Center of Mass or [Fréchet mean](https://en.wikipedia.org/wiki/Fréchet_mean):
For given points ``p_i ∈\mathcal M``, ``i=1,…,N`` this optimization problem reads

```math
\operatorname*{arg\,min}_{x∈\mathcal M} \frac{1}{2}\sum_{i=1}^{N}
  \operatorname{d}^2_{\mathcal M}(x,p_i),
```

which of course can be (and is) solved by a gradient descent, see the introductionary tutorial or the [Statistics in Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/stable/features/statistics.html).
If ``N`` is very large it might be quite expensive to evaluate the complete gradient.
A remedy is, to evaluate only one of the terms at a time and choose a random order for these.

We first initialize the packages
"""

# ╔═╡ 84528b83-f11d-42de-9e86-fc9f3d451500
md"and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)"

# ╔═╡ c104cad0-c489-42f8-b896-986811eebfef
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733") # Start
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB") # a path
    TolVibrantTeal = RGBA{Float64}(colorant"#009988") # points
end;

# ╔═╡ 11ffe7e1-22f5-4098-9830-d56ff59745ee
md"We next generate a (little) large(r) data set"

# ╔═╡ f3a29d0e-91a5-437f-8681-5b72072b9cdf
begin
    n = 5000
    σ = π / 12
    M = Sphere(2)
    x = 1 / sqrt(2) * [1.0, 0.0, 1.0]
    Random.seed!(42)
    data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # files folder
    image_prefix = localpath * "/stochastic_gradient_descent"
    @info image_prefix
    render_asy = false # on CI or when you do not have asymptote, this should be false
end

# ╔═╡ 4f527a26-dc55-4eb9-ab1b-62a4712182d5
render_asy && asymptote_export_S2_signals(
    image_prefix * "/center_and_large_data.asy";
    points=[[x], data],
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
    dot_sizes=[2.5, 1.0],
    camera_position=(1.0, 0.5, 0.5),
)

# ╔═╡ 8baf553d-8d22-4f13-9741-43531b080998
render_asy && render_asymptote(image_prefix * "/center_and_large_data.asy"; render=2);

# ╔═╡ daf0dc47-cba2-4264-84d9-47137924addd
PlutoUI.LocalResource(image_prefix * "/center_and_large_data.png")

# ╔═╡ 30f0a209-afde-40c4-a8ab-3f40fd374d3c
md"""
Note that due to the construction of the points as zero mean tangent vectors, the mean should
be very close to our initial point `x`.

In order to use the stochastic gradient, we now need a function that returns the vector of gradients.
There are two ways to define it in `Manopt.jl`: as one function, that returns a vector or a vector of funtions.

The first variant is of course easier to define, but the second is more efficient when only evaluating one of the gradients.

For the mean we have as a gradient

```math
 gradF(x) = \sum_{i=1}^N \operatorname{grad}f_i(x) \quad \text{where} \operatorname{grad}f_i(x) = -\log_x p_i
```

Which we define in `Manopt.jl` in two different ways:
Either as one function returning all gradients as a vector (see `gradF`) or – maybe more fitting for a large scale problem, as a vector of small gradient functions (see `gradf`)
"""

# ╔═╡ 12a60ec6-cd18-4eea-a0af-0cabd43686b8
F(M, x) = 1 / (2 * n) * sum(map(p -> distance(M, x, p)^2, data))

# ╔═╡ d76aa8ef-b234-40d6-940b-cfc5db38cee8
gradF(M, x) = [grad_distance(M, p, x) for p in data]

# ╔═╡ d41271a4-cbbe-46a3-b350-fce79d860671
gradf = [(M, x) -> grad_distance(M, p, x) for p in data];

# ╔═╡ 2d977c61-c556-4d58-8ca6-de0f7977cd30
md"""
The calls are only slightly different, but notice that accessing the 2nd gradient element
requires evaluating all logs in the first function.
So while you can use both `gradF` and `gradf` in the following call, the second one is (much) faster:
"""

# ╔═╡ 2b1ed0d6-d64c-41a2-b396-c3bb243af79c
x_opt1 = stochastic_gradient_descent(M, gradF, x)

# ╔═╡ a6d7707d-8d8b-4136-b058-60c1a33b182d
@benchmark stochastic_gradient_descent($M, $gradF, $x)

# ╔═╡ 201d9ab2-2848-472e-abd6-826301c328e3
x_opt2 = stochastic_gradient_descent(M, gradf, x)

# ╔═╡ ce5e6f7b-f101-401f-9470-74ae94b87096
@benchmark stochastic_gradient_descent($M, $gradf, $x)

# ╔═╡ b463184b-3b66-4c1d-810c-bb398d918e2a
x_opt2

# ╔═╡ f2093d02-a6e0-49f3-ad39-55d9d88e75d0
md"""This result is reasonably close. But we can improve it by using a `DirectionUpdateRule`, namely:

On the one hand [`MomentumGradient`](https://manoptjl.org/stable/solvers/gradient_descent.html#Manopt.MomentumGradient), which requires both the manifold and the initial value,  in order to keep track of the iterate and parallel transport the last direction to the current iterate.
You can also set a `vector_transport_method`, if `ParallelTransport()` is not
available on your manifold. Here we simply do
"""

# ╔═╡ 22c66c8a-3db1-4584-8cb1-93fbbeb34d4e
x_opt3 = stochastic_gradient_descent(
    M, gradf, x; direction=MomentumGradient(M, x, StochasticGradient(zero_vector(M, x)))
)

# ╔═╡ dbefa2af-d68f-4918-8ab0-4af257425686
MG = MomentumGradient(M, x, StochasticGradient(zero_vector(M, x)));

# ╔═╡ 5d5122b6-ba71-48ef-b459-d676ad4cca0d
@benchmark stochastic_gradient_descent($M, $gradf, $x; direction=$MG)

# ╔═╡ d7526835-73d3-40e6-8f8d-1711519ee57a
md"""
And on the other hand the [`AverageGradient`](https://manoptjl.org/stable/solvers/gradient_descent.html#Manopt.AverageGradient) computes an average of the last `n` gradients, i.e.
"""

# ╔═╡ a23ee4d6-6e9b-413a-a3ee-75367c9ea943
x_opt4 = stochastic_gradient_descent(
    M, gradf, x; direction=AverageGradient(M, x, 10, StochasticGradient(zero_vector(M, x)))
);

# ╔═╡ b0daf769-f8ea-4d7d-9bd4-48748bc02346
AG = AverageGradient(M, x, 10, StochasticGradient(zero_vector(M, x)));

# ╔═╡ 8b7cbfc8-7ec5-4297-b8f1-d3691725578e
@benchmark stochastic_gradient_descent($M, $gradf, $x; direction=$AG)

# ╔═╡ 5a3cb657-b25b-4622-a5a1-c4e93a17204f
md"""
Note that the default `StoppingCriterion` is a fixed number of iterations which helps the comparison here.

For both update rules we have to internally specify that we are still in the Stochastic setting (since both rules can also be used with the `IdentityUpdateRule` within [`gradient_descent`](file:///Users/ronny/Repositories/Julia/Manopt.jl/docs/build/solvers/gradient_descent.html),

For this not-that-large-scale example we can of course also use a gradient descent with `ArmijoLinesearch`, but it will be a little slower usually
"""

# ╔═╡ b60b65a4-a46f-4f12-a416-48a716fa7386
fullGradF(M, x) = sum(grad_distance(M, p, x) for p in data)

# ╔═╡ ed9d3864-9a24-4c9a-938b-c2d3c59bb97e
x_opt5 = gradient_descent(M, F, fullGradF, x; stepsize=ArmijoLinesearch())

# ╔═╡ 534bbae1-c777-460b-bdc3-9d4b78132644
AL = ArmijoLinesearch();

# ╔═╡ dba33908-3a1b-491a-abd8-497a00d85230
@benchmark gradient_descent($M, $F, $fullGradF, $x; stepsize=$AL)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
BenchmarkTools = "~1.3.1"
Colors = "~0.12.8"
Manifolds = "~0.7.8"
Manopt = "~0.3.23"
PlutoUI = "~0.7.38"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "c933ce606f6535a7c7b98e1d86d5d1014f730596"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.7"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.3.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.10.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "a970d55c2ad8084ca317a4658ba6ce99b7523571"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.12"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "cc5e34ecd061bf9bedd62a8b8a28abe96c3b7c6e"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.8"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "b8518cbf0f06f1375c95479759b266bd75d89ad9"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.13"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "a8a292a422b0a6924f46d4e022fc8263721120c5"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.23"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "3b429f37de37f1fc603cc1de4a799dc7fbe4c0b6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "bfe14f127f3e7def02a6c2b1940b39d0dabaa3ef"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.26.3"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.3.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─f58f28fe-c3d8-11ec-0acb-e948605dd63f
# ╠═8a8be09a-49e9-4563-9afc-ae72bd05d26d
# ╟─84528b83-f11d-42de-9e86-fc9f3d451500
# ╠═c104cad0-c489-42f8-b896-986811eebfef
# ╟─11ffe7e1-22f5-4098-9830-d56ff59745ee
# ╠═f3a29d0e-91a5-437f-8681-5b72072b9cdf
# ╠═4f527a26-dc55-4eb9-ab1b-62a4712182d5
# ╠═8baf553d-8d22-4f13-9741-43531b080998
# ╠═daf0dc47-cba2-4264-84d9-47137924addd
# ╟─30f0a209-afde-40c4-a8ab-3f40fd374d3c
# ╠═12a60ec6-cd18-4eea-a0af-0cabd43686b8
# ╠═d76aa8ef-b234-40d6-940b-cfc5db38cee8
# ╠═d41271a4-cbbe-46a3-b350-fce79d860671
# ╟─2d977c61-c556-4d58-8ca6-de0f7977cd30
# ╠═2b1ed0d6-d64c-41a2-b396-c3bb243af79c
# ╠═a6d7707d-8d8b-4136-b058-60c1a33b182d
# ╠═201d9ab2-2848-472e-abd6-826301c328e3
# ╠═ce5e6f7b-f101-401f-9470-74ae94b87096
# ╠═b463184b-3b66-4c1d-810c-bb398d918e2a
# ╟─f2093d02-a6e0-49f3-ad39-55d9d88e75d0
# ╠═22c66c8a-3db1-4584-8cb1-93fbbeb34d4e
# ╠═dbefa2af-d68f-4918-8ab0-4af257425686
# ╠═5d5122b6-ba71-48ef-b459-d676ad4cca0d
# ╟─d7526835-73d3-40e6-8f8d-1711519ee57a
# ╠═a23ee4d6-6e9b-413a-a3ee-75367c9ea943
# ╠═b0daf769-f8ea-4d7d-9bd4-48748bc02346
# ╠═8b7cbfc8-7ec5-4297-b8f1-d3691725578e
# ╟─5a3cb657-b25b-4622-a5a1-c4e93a17204f
# ╠═b60b65a4-a46f-4f12-a416-48a716fa7386
# ╠═ed9d3864-9a24-4c9a-938b-c2d3c59bb97e
# ╠═534bbae1-c777-460b-bdc3-9d4b78132644
# ╠═dba33908-3a1b-491a-abd8-497a00d85230
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    PedestrianLatLonAccel
Pedestrian walking action. Acceleration in the Frenet frame, along with
desired lane after crossing the street.

# Fields
- `a_lat::Float64` lateral acceleration [m/s^2]
- `a_lon::Float64` longitudinal accelaration [m/s^2]
- `lane_des::Lane` desired lane to move to
"""
struct PedestrianLatLonAccel
    a_lat::Float64
    a_lon::Float64
    lane_des::Lane
end

function AutomotiveDrivingModels.propagate(ped::Entity{VehicleState, D, I},
                                           action::PedestrianLatLonAccel,
                                           roadway::Roadway,
                                           ΔT::Float64=0.1
                                           ) where {D,I}
    state = propagate(ped, LatLonAccel(action.a_lat, action.a_lon), roadway, ΔT)
    roadproj = proj(posg(state), roadway[action.lane_des.tag], roadway, move_along_curves=false)
    retval = VehicleState(Frenet(roadproj, roadway), roadway, vel(state))
    return retval
end


text:
# ---
# title: 583. Delete Operation for Two Strings
# id: problem583
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/delete-operation-for-two-strings/description/>
# hidden: true
# ---
#
# Given two words _word1_ and _word2_ , find the minimum number of steps
# required to make _word1_ and _word2_ the same, where in each step you can
# delete one character in either string.
#
# **Example 1:**
#
#
#
#     Input: "sea", "eat"
#     Output: 2
#     Explanation: You need one step to make "sea" to "ea" and another step to make "eat" to "ea".
#
#
# **Note:**
#
#   1. The length of given words won't exceed 500.
#   2. Characters in given words can only be lower-case letters.
#
#
## @lc code=start
using LeetCode

min_distance583(word1::String, word2::String) =
    length(word2) + length(word2) - 2longest_common_subsequence(word1, word2)
## @lc code=end


text:
struct ECE{B<:AbstractBinningAlgorithm,D} <: CalibrationErrorEstimator
    """Binning algorithm."""
    binning::B
    """Distance function."""
    distance::D
end

"""
    ECE(binning[, distance = TotalVariation()])

Estimator of the expected calibration error (ECE) for a classification model with
respect to the given `distance` function using the `binning` algorithm.

For classification models, the predictions ``P_{X_i}`` and targets ``Y_i`` are identified
with vectors in the probability simplex. The estimator of the ECE is defined as
```math
\\frac{1}{B} \\sum_{i=1}^B d\\big(\\overline{P}_i, \\overline{Y}_i\\big),
```
where ``B`` is the number of non-empty bins, ``d`` is the distance function, and
``\\overline{P}_i`` and ``\\overline{Y}_i`` are the average vector of the predictions and
the average vector of targets in the ``i``th bin. By default, the total variation distance
is used.

The `distance` has to be a function of the form
```julia
distance(pbar::Vector{<:Real}, ybar::Vector{<:Real}).
```
In particular, distance measures of the package
[Distances.jl](https://github.com/JuliaStats/Distances.jl) are supported.
"""
ECE(binning::AbstractBinningAlgorithm) = ECE(binning, TotalVariation())

# estimate ECE
function (ece::ECE)(
    predictions::AbstractVector{<:AbstractVector{<:Real}},
    targets::AbstractVector{<:Integer},
)
    @unpack binning, distance = ece

    # check number of samples
    nsamples = check_nsamples(predictions, targets)

    # bin predictions and labels
    bins = perform(binning, predictions, targets)
    nbins = length(bins)
    nbins > 0 || error("there must exist at least one bin")

    # compute the weighted mean of the distances in each bin
    # use West's algorithm for numerical stability

    # evaluate the distance in the first bin
    @inbounds begin
        bin = bins[1]
        x = distance(bin.mean_predictions, bin.proportions_targets)

        # initialize the estimate
        estimate = x / 1

        # for all other bins
        n = bin.nsamples
        for i in 2:nbins
            # evaluate the distance
            bin = bins[i]
            x = distance(bin.mean_predictions, bin.proportions_targets)

            # update the estimate
            m = bin.nsamples
            n += m
            estimate += (m / n) * (x - estimate)
        end
    end

    return estimate
end


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

using FEMBase, MortarContact2D, Test
using MortarContact2D: create_contact_segmentation, calculate_normals

X = Dict(1 => [0.0, 0.0],
         2 => [1.0, 0.0],
         3 => [0.0, 1.0],
         4 => [1.0, 1.0])

slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
elements = [slave, master]
update!(elements, "geometry", X)
update!(slave, "displacement", 0.0 => ([0.0,0.0], [0.0,0.0]))
update!(master, "displacement", 0.0 => ([0.0,0.0], [0.0,0.0]))
update!(slave, "displacement", 1.0 => ([0.0,0.0], [0.0,0.0]))
update!(master, "displacement", 1.0 => ([2.0,0.0], [2.0,0.0]))
problem = Problem(Contact2D, "test problem", 2, "displacement")
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
normals = calculate_normals([slave], 0.0)
update!([slave], "normal", 0.0 => normals)
seg1 = create_contact_segmentation(problem, slave, [master], 0.0; deformed=true)
seg2 = create_contact_segmentation(problem, slave, [master], 1.0; deformed=true)
problem.properties.max_distance = 0.0
seg3 = create_contact_segmentation(problem, slave, [master], 0.0; deformed=false)
@test length(seg1) == 1
@test length(seg2) == 0
@test length(seg3) == 0

slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
elements = [slave, master]
update!(elements, "geometry", X)
problem = Problem(Contact2D, "test problem", 2, "displacement")
problem.assembly.la = zeros(8)
problem.properties.store_fields = [
    "contact area",
    "contact error",
    "weighted gap",
    "contact pressure",
    "complementarity condition",
    "active nodes",
    "inactive nodes",
    "stick nodes",
    "slip nodes"]
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
assemble!(problem, 0.0)

X[3] += [2.0, 0.0]
X[4] += [2.0, 0.0]
slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
elements = [slave, master]
update!(elements, "geometry", X)
problem = Problem(Contact2D, "test problem", 2, "displacement")
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
assemble!(problem, 0.0)

X = Dict(1 => [0.0, 0.0],
         2 => [1.0, 0.0],
         3 => [0.0, 0.0],
         4 => [1.0, 0.0])

slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
elements = [slave, master]
update!(elements, "geometry", X)
problem = Problem(Contact2D, "test problem", 2, "displacement")
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
assemble!(problem, 0.0)
empty!(problem.assembly)
problem.properties.iteration = 0
problem.properties.contact_state_in_first_iteration = :INACTIVE
assemble!(problem, 0.0)


text:
import Base: length, write, read, readall, readuntil, readbytes

const PROTOCOL_VERSION = uint32(0x0003_0000)
const CANCEL_CODE      = uint32(0x04d2_162e)
const SSL_CODE         = uint32(0x04d2_162f)

write_be(io, x) = write(io, hton(x))
read_be(io, T) = ntoh(read(io, T))

immutable MSG{I}
  data :: IOBuffer

  MSG() = new(PipeBuffer())
  MSG(data :: IOBuffer) = new(data)
end

# Definition taken from
# https://github.com/JuliaLang/julia/blob/04893a165aed19d898a9ed2f9dc2202553906256/base/iobuffer.jl#L106
function length(msg :: MSG)
  io = data(msg)
  (io.seekable ? io.size : nb_available(io))
end
data(msg :: MSG) = msg.data
ident{I}(:: MSG{I}) = I

write(m :: MSG, x) = write(data(m), x)
read(m :: MSG, x) = read(data(m), x)
read(m :: MSG, x, l) = read(data(m), x, l)
readall(m :: MSG) = readall(data(m))
readuntil(m :: MSG, c) = readuntil(data(m), c)
readbytes(m :: MSG) = readbytes(data(m))

readstring(m :: MSG) = readuntil(m, '\0')[1:end-1]

pg_msg(i :: Symbol, args...) = pg_msg(MSG{i}(), args...)
pg_msg(i :: MSG, args...) = error("Could not construct message for ident: $(ident(i))")

function pg_msg(msg :: MSG{:start}, options :: Dict{ByteString, ByteString}, version = PROTOCOL_VERSION)
  write_be(msg, version)

  for (key, value) in options
    write(msg, key)
    write(msg, '\0')
    write(msg, value)
    write(msg, '\0')
  end

  write(msg, '\0')
  msg
end

function writemsg(conn, msg :: MSG)
  out = PipeBuffer()
  id = ident(msg)
  if id != :start
    write(out, id)
  end

  write_be(out, int32(length(msg) + sizeof(Int32))) # Length - big endian
  write(out, data(msg))
  write(conn, takebuf_array(out))
end

function readmsg(conn)
  ident = symbol(read(conn, Char))
  length = read_be(conn, Int32) - sizeof(Int32)
  data = IOBuffer(read(conn, UInt8, length))
  MSG{ident}(data)
end

##
# Function: parsemsg
# Parses a postgres Backend message after:
# http://www.postgresql.org/docs/9.4/static/protocol-message-formats.html
#
# Each message is identified by a ident char
##
function parsemsg(msg :: MSG)
  error("Can't parse postgres message with ident $(ident(msg))")
end

function parsemsg(msg :: MSG{:R})
  request = read_be(msg, Int32)

  if request == 0
    return (:AuthenticationOK, )
  elseif request == 2
    return (:AuthenticationKerberosV5, )
  elseif request == 3
    return (:AuthenticationClearTextPassword, )
  elseif request == 5
    salt = read(msg, UInt8, 4)
    return (:AuthenticationMD5Password, salt)
  elseif request == 6
    return (:AuthenticationSCMCredential, )
  elseif request == 7
    return (:AuthenticationGSS, )
  elseif request == 9
    return (:AuthenticationSSPI, )
  elseif request == 8
    data = read(msg, UInt8, length(msg))
    return (:AuthenticationGSSContinue, data)
  else
    error("Does not support AuthenticationRequest $request")
  end
end

function parsemsg(msg :: MSG{:K})
  id  = read_be(msg, Int32)
  key = read_be(msg, Int32)
  :BackendKeyData, id, key
end

#TODO does this work
parsemsg(:: MSG{symbol('2')}) = (:BindComplete, )

#TODO does this work
parsemsg(:: MSG{symbol('3')}) = (:CloseComplete, )

parsemsg(msg :: MSG{:C}) = (:CommandComplete, readstring(msg))

function parsemsg(msg :: MSG{:d})
  payload = read(msg, UInt8, length(msg))

  :CopyData, payload
end

parsemsg(msg :: MSG{:c}) = (:CopyDone, )

function parsemsg(msg :: MSG{:G})
  format = read_be(msg, Int8) == 0 ? :textual : :binary
  n = read_be(msg, Int16)
  columns = Array(Int16, n)
  for i in 1:n
    columns[i] = read_be(msg, Int16)
  end

  :CopyInResponse, format, columns
end

function parsemsg(msg :: MSG{:H})
  format = read_be(msg, Int8) == 0 ? :textual : :binary
  n = read_be(msg, Int16)
  columns = Array(Int16, n)
  for i in 1:n
    columns[i] = read_be(msg, Int16)
  end

  :CopyOutResponse, format, columns
end

function parsemsg(msg :: MSG{:W})
  format = read_be(msg, Int8) == 0 ? :textual : :binary
  n = read_be(msg, Int16)
  columns = Array(Int16, n)
  for i in 1:n
    columns[i] = read_be(msg, Int16)
  end

  :CopyBothResponse, format, columns
end

function parsemsg(msg :: MSG{:D})
  n = read_be(msg, Int16)
  columns = Array(Nullable{Vector{UInt8}}, n)

  for i in 1:n
    c = read_be(msg, Int32)

    if c == -1
      val = Nullable{Vector{UInt8}}()
    else
      val = Nullable(read(msg, UInt8, c))
    end

    columns[i] = val
  end

  :DataRow, columns
end

parsemsg(msg :: MSG{:I}) = (:EmptyQueryResponse, )

parsemsg(msg :: MSG{:E}) = (:ErrorResponse, DBError(msg))

function parsemsg(msg :: MSG{:V})
  lr = read_be(msg, Int32)
  if lr == -1
    val = Nullable{Vector{UInt8}}()
  else
    val = Nullable(read(msg, UInt8, lr))
  end

  :FunctionCallResponse, val
end

parsemsg(msg :: MSG{:n}) = (:NoData, )

parsemsg(msg :: MSG{:N}) = (:NoticeResponse, DBError(msg))

function parsemsg(msg :: MSG{:A})
  id = read_be(msg, Int32)
  name = readstring(msg)
  payload = readstring(msg)

  :NotificationResponse, id, name, payload
end

function parsemsg(msg :: MSG{:t})
  n = read_be(msg, Int16)
  params = Array(Int32, n)
  for i in 1:n
    params[i] = read_be(msg, Int32)
  end

  :ParameterDescription, params
end

function parsemsg(msg :: MSG{:S})
  name = readstring(msg)
  value = readstring(msg)

  :ParameterStatus, name, value
end

#TODO: Does this work
parsemsg(msg :: MSG{symbol('1')}) = (:ParseComplete, )

parsemsg(msg :: MSG{:s}) = :PortalSuspended

function parsemsg(msg :: MSG{:Z})
  status = read(msg, Char)
  :ReadyForQuery, status
end

function parsemsg(msg :: MSG{:B})
  n = read_be(msg, Int16)
  fields = Array{Field, n}

  for i in 1:n
    name          = readstring(msg)
    tableID       = read_be(msg, Int32)
    columnID      = read_be(msg, Int16)
    datatypeID    = read_be(msg, Int32)
    datatypeSize  = read_be(msg, Int16)
    typeModifier  = read_be(msg, Int32)
    format        = read_be(msg, Int16) == 0 ? :textual : :binary

    fields[i] = Field{format}(name, tableID, columnID, datatypeID,
                              datatypeSize, typeModifier)
  end

  :RowDescription, fields
end


text:
tests = ("constructors.jl",)

for test in tests
    try
        include(test)
    catch err
        println("Test $test failed.")
        throw(err)
    end
end


text:
"""
Render an inline R script, substituting invalid "\$" signs for Julia symbols
"""
function render(script::Compat.String)
    symdict = OrderedDict{Compat.String,Any}()
    local status
    local msg = ""
    local k = 0
    while true
        status = parseVector(sexp(script))[2]
        if status != 1
            msg = getParseErrorMsg()
        end

        # break if not parse error (status = 3)
        (status != 3) && break

        # due to a bug in the R parser https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=16524
        # parse and getParseData do not work with unicode
        # R_ParseContextLast and R_ParseContext are not documentated, but they seem to work

        # the position of the error byte
        b = Int(unsafe_load(cglobal((:R_ParseContextLast, libR), Cint)))
        # it is unsafe to use `unsafe_string` if the last byte is a part of a unicode,
        c = Char(unsafe_load(cglobal((:R_ParseContext, libR), Cchar)+b))
        c != '\$' && break

        ast,i = parse(script,b+1,greedy=false,raise=false)

        if isa(ast,Symbol)
            sym = "$ast"
        elseif isa(ast, Expr) && !(ast.head == :error || ast.head == :continue || ast.head == :incomplete)
            sym = "($ast)"
            # if an expression has already appeared, we generate a new symbol so it will be evaluated twice (e.g. `R"$(rand(10)) == $(rand(10))"`)
            if haskey(symdict, sym)
                sym *= "##$k"
                k = k + 1
            end
        elseif isa(ast, Expr) && (ast.head == :error || ast.head == :continue)
            status = 3
            msg = ast.args[1]
            break
        elseif isa(ast, Expr) && ast.head == :incomplete
            status = 2
            msg = ast.args[1]
            break
        else
            status = 3
            msg = "unknown render error"
            break
        end
        symdict[sym] = ast
        script = string(script[1:b-1],"`#JL`\$`",sym,'`',script[i:end])
    end

    return script, symdict, status, msg
end