text:
# Biryani -- A conversion and validation toolbox
# By: Emmanuel Raviart <emmanuel@raviart.com>
#
# Copyright (C) 2015 Emmanuel Raviart
# https://github.com/eraviart/Biryani.jl
#
# This file is part of Biryani.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

abstract Context

type Convertible
  value
  context::Context
  error
end

Convertible(value, context::Context) = Convertible(value, context, nothing)
Convertible(value) = Convertible(value, EmptyContext())

type EmptyContext <: Context
end

_(context::EmptyContext, message::AbstractString) = message

N_(message) = message

function call(func::Function; handle_nothing = false)
  """Return a converter that applies a function to value and returns a new value.

  .. note:: Like most converters, by default a ``nothing`` value is not converted (ie function is not
     called). Set ``handle_nothing`` to ``true`` to call function when value is ``nothing``.

  .. note:: When your function doesn't modify value but may generate an error, use a :func:`test` instead.

  .. note:: When your function modifies value and may generate an error, write a full converter instead of a function.

  See :doc:`how-to-create-converter` for more information.
  """
  return function run_call(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing && !handle_nothing
      return convertible
    end
    return Convertible(func(convertible.value), convertible.context)
  end
end

function condition(converters::Function...)
  """When first converter succeeds (ie no error), then apply the second converter, otherwise test the third converter...

  When the number of converters is odd, the last converter is applied when all tests have failed (ie: its is an "else"
  converter).

  .. note:: See also :func:`first_match`.
  """
  return function run_condition(convertible::Convertible)
    last_index = length(converters)
    ok = false
    for (index, converter) in enumerate(converters)
      if index & 1 == 0
        if ok
          return converter(convertible)
        end
      elseif index == last_index
        # Else converter
        return converter(convertible)
      elseif converter(convertible).error === nothing
        # Test converter is OK.
        ok = true
      end
    end
    # Every test has failed and there is no else clause. Return the input convertible.
    return convertible
  end
end

function default(value)
  """Return a converter that replace a ``nothing`` value with given one.

  .. note:: See converter :func:`from_value` to replace a non-``default`` value.
  """
  return function run_default(convertible::Convertible)
    if convertible.error !== nothing || convertible.value !== nothing
      return convertible
    end
    return Convertible(value, convertible.context)
  end
end

function embed_error(value, error)
  """Embed errors inside value.

  .. note:: value is modified. No copy is made.
  """
  if error === nothing
    return nothing
  end
  if isa(value, Union(Dict, OrderedDict))
    if isa(error, Union(Dict, OrderedDict))
      for (child_key, child_error) in error
        child_error = embed_error(get(value, child_key, nothing), child_error)
        if child_error !== nothing
          errors = get!(value, "errors", Dict{typeof(child_key), Any}())
          errors[child_key] = child_error
        end
      end
    else
      value["errors"] = error
    end
    return nothing
  end
  if isa(value, Array) && isa(error, Union(Dict, OrderedDict))
    if all(key -> isa(key, Integer) && 1 <= key <= size(value, 1), error)
      for (child_index, child_error) in error
        child_error = embed_error(value[child_index], child_error)
        @assert(child_error === nothing, child_error)
        # if child_error !== nothing
        #   return error
        # end
      end
      return nothing
    end
    if all(key -> isa(key, AbstractString) && isdigit(key) && 1 <= parse(Int, key) <= size(value, 1), error)
      for (child_key, child_error) in error
        child_error = embed_error(value[parse(Int, child_key)], child_error)
        @assert(child_error === nothing, child_error)
        # if child_error !== nothing
        #   return error
        # end
      end
      return nothing
    end
  end
  return error
end

function empty_to_nothing(convertible::Convertible)
  """When value is an empty collection, replace it with ``nothing`` else keep it as is."""
  if convertible.error !== nothing || convertible.value === nothing || !isempty(convertible.value)
    return convertible
  end
  return Convertible(nothing, convertible.context)
end

eval_error(context::Context, error_by_key::Dict) = [
  key => eval_error(context, value)
  for (key, value) in error_by_key
]

function eval_error(context::Context, error_by_key::OrderedDict)
  evaluated_error_by_key = similar(error_by_key)
  for (key, value) in error_by_key
    evaluated_error_by_key[key] = eval_error(context, value)
  end
  return evaluated_error_by_key
end

eval_error(context::Context, func::Function) = func(context)

eval_error(context::Context, ::Void) = nothing

eval_error(context::Context, message::AbstractString) = _(context, message)

extract_when_singleton(convertible::Convertible) = condition(
  test(value -> length(value) == 1),
  call(value -> value[1]),
)(convertible)
"""Extract first item of sequence when it is a singleton, otherwise keep it unchanged."""

function fail(convertible::Convertible; error = nothing)
  """Always return an error."""
  if convertible.error !== nothing
    return convertible
  end
  return Convertible(
    convertible.value,
    convertible.context,
    error === nothing ? N_("An error occured.") : error,
  )
end

fail(error::AbstractString) = convertible::Convertible -> fail(convertible, error = error)

fail(; error::AbstractString = nothing) = convertible::Convertible -> fail(convertible, error = error)

function first_match(converters::Function...; error = nothing)
  """Try each converter successively until one succeeds. When every converter fail, return an error.

  .. note:: See also :func:`condition`.
  """
  return function run_first_match(convertible::Convertible)
    if isempty(converters)
      return convertible
    end
    converted = convertible
    for converter in converters
      converted = converter(convertible)
      if converted.error === nothing
        return converted
      end
    end
    return error === nothing ? converted : Convertible(convertible.value, convertible.context, error)
  end
end

function from_value(value)
  """Return a converter that replace any value with given one."""
  return function run_from_value(convertible::Convertible)
    if convertible.error !== nothing
      return convertible
    end
    return Convertible(value, convertible.context)
  end
end

function guess_bool(convertible::Convertible)
  """Convert the content of a string (or a number) to a boolean. Do nothing when input value is already a boolean.

  This converter accepts usual values for ``true`` and ``false``: "0", "f", "false", "n", etc.

  .. warning:: Like most converters, a ``nothing`` value is not converted.

      When you want ``nothing`` to be converted to ``false``, use::

        pipe(guess_bool, default(False))
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  if !isa(convertible.value, AbstractString)
    return to_bool(convertible)
  end
  lower_value = strip(convertible.value) |> lowercase
  if isempty(lower_value)
    return Convertible(nothing, convertible.context)
  end
  if lower_value in ("0", "f", "false", "n", "no", "off")
    value = false
  elseif lower_value in ("1", "on", "t", "true", "y", "yes")
    value = true
  else
    return Convertible(lower_value, convertible.context, N_("Value must be a boolean."))
  end
  return Convertible(value, convertible.context)
end

input_to_bool(convertible::Convertible) = pipe(strip, to_bool)(convertible)
"""Convert a string to a boolean.

.. warning:: Like most converters, a ``nothing`` value is not converted.

  When you want ``nothing`` to be converted to ``false``, use::

    pipe(input_to_bool, default(false))
"""

input_to_email(convertible::Convertible; accept_ip_address = false) = pipe(
  strip,
  string_to_email(accept_ip_address = accept_ip_address),
)(convertible)
"""Convert a string to an email address."""

input_to_email(; accept_ip_address = false) = convertible::Convertible -> input_to_email(convertible,
  accept_ip_address = accept_ip_address)

input_to_float(convertible::Convertible; accept_expression = false) = pipe(
  strip,
  to_float(accept_expression = accept_expression),
)(convertible)
"""Convert a string to a float."""

function input_to_float(; accept_expression = false)
  return convertible::Convertible -> input_to_float(convertible; accept_expression = accept_expression)
end

input_to_int(convertible::Convertible; accept_expression = false) = pipe(
  strip,
  to_int(accept_expression = accept_expression),
)(convertible)
"""Convert a string to an integer number."""

function input_to_int(; accept_expression = false)
  return convertible::Convertible -> input_to_int(convertible; accept_expression = accept_expression)
end

input_to_ipv4(convertible::Convertible) = pipe(
  strip,
  string_to_ipv4,
)(convertible)
"""Convert a string to an IPv4 address."""

function input_to_url_name(convertible::Convertible; separator = '_', stripmark = false)
  """Normalize a string to allow its use in an URL (or file system) path or a query parameter.

  .. note:: For a converter that keeps only letters, digits and separator, see :func:`input_to_slug`.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  value = normalize_string(convertible.value, casefold = true, compat = true, decompose = true, stable = true,
    stripcc = true, stripignore = true, stripmark = stripmark)
  # Replace unsafe characters (for URLs and file-systems).
  value = replace(value, ['\n', '\r', '\\', '/', ';', ':', '"', '#', '*', '?', '&', '<', '>', '|', '.'], ' ')
  value = join(split(value), separator)
  return Convertible(isempty(value) ? nothing : value, convertible.context)
end

function input_to_url_name(; separator = '_', stripmark = false)
  return convertible::Convertible -> input_to_url_name(convertible, separator = separator, stripmark = stripmark)
end

is_valid(convertible::Convertible) = convertible.error === nothing
"""Check a conversion and return true when it has no error, false otherwise."""

function item_or_sequence(converters::Function...; drop_nothing = false, item_type = nothing, sequence_type = Array)
  """Return a converter that accepts either an item or a sequence of items and applies the given converter to them."""
  return convertible::Convertible -> condition(
    test_isa(sequence_type),
    pipe(
      uniform_sequence(converters..., drop_nothing = drop_nothing, item_type = item_type,
        sequence_type = sequence_type),
      extract_when_singleton,
    ),
    pipe(converters...),
  )(convertible)
end

function item_to_singleton(convertible::Convertible; sequence_type = Array):
  """Convert an item to a singleton, but keep a sequence of items unchanged."""
  if convertible.error !== nothing || convertible.value === nothing || isa(convertible.value, sequence_type)
    return convertible
  end
  # TODO: Use sequence_type to build result value.
  return Convertible([convertible.value], convertible.context)
end

function log_info(convertible::Convertible, messages...)
  """Display value using "info" function."""
  info([string(message) for message in messages]..., "value: ", string(convertible.value), ", type: ",
    string(typeof(convertible.value)), ", error: ", string(convertible.error))
  return convertible
end

log_info(messages...) = convertible::Convertible -> log_info(convertible, messages...)

function log_warning(convertible::Convertible, messages...)
  """Display value using "warn" function."""
  warn([string(message) for message in messages]..., "value: ", string(convertible.value), ", error: ",
    string(convertible.error))
  return convertible
end

log_warning(messages...) = convertible::Convertible -> log_warning(convertible, messages...)

noop(convertible::Convertible) = convertible

function pipe(converters::Function...)
  """Return a compound converter that applies each of its converters till the end or an error occurs."""
  return function run_pipe(convertible::Convertible)
    for converter in converters
      convertible = converter(convertible)
      @assert(typeof(convertible) <: Convertible)
    end
    return convertible
  end
end

require(convertible::Convertible) = convertible.error === nothing && convertible.value === nothing ?
  Convertible(convertible.value, convertible.context, N_("Missing value.")) :
  convertible
"""Return an error when value is ``nothing``."""

function string_to_email(convertible::Convertible; accept_ip_address = false)
  """Convert a clean string to an email address.

  .. note:: For a converter that doesn't require a clean string, see :func:`input_to_email`.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  value = lowercase(convertible.value)
  if startswith(value, "mailto:")
    value = replace(value, "mailto:", "", 1)
  end
  split_value = split(value, '@')
  if length(split_value) != 2
    return Convertible(value, convertible.context, N_("""An email must contain exactly one "@"."""))
  end
  username, domain = split_value
  if !ismatch(r"[^ \t\n\r@<>()]+$"i, username)
    return Convertible(value, convertible.context, N_("Invalid username."))
  end
  if domain != "localhost" &&
      !ismatch(r"""
        ^(?:[a-z0-9][a-z0-9\-]{0,62}\.)+ # (sub)domain - alpha followed by 62max chars (63 total)
        [a-z]{2,}$                      # TLD
      """ix, domain) &&
      (!accept_ip_address || !ismatch(r"^\d+\.\d+\.\d+\.\d+$", domain))
    return Convertible(value, convertible.context, N_("Invalid domain name."))
  end
  return Convertible(value, convertible.context)
end

string_to_email(; accept_ip_address = false) = convertible::Convertible -> string_to_email(convertible,
  accept_ip_address = accept_ip_address)

function string_to_ipv4(convertible::Convertible)
  """Convert a clean string to an IPv4 address.

  .. note:: For a converter that doesn't require a clean string, see :func:`input_to_ipv4`.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  if !ismatch(r"^\d+\.\d+\.\d+\.\d+$", convertible.value)
    return Convertible(convertible.value, convertible.context, N_("Invalid IPv4 address."))
  end
  bytes = map(value -> parse(Int, value), split(convertible.value, '.'))
  if !all(byte -> 0 <= byte <= 255, bytes)
    return Convertible(convertible.value, convertible.context, N_("Invalid IPv4 address."))
  end
  return Convertible(IPv4(bytes...), convertible.context)
end

strip(convertible::Convertible) = pipe(call(strip), empty_to_nothing)(convertible)
"""Strip spaces from a string and remove it when empty."""

function struct(converters::Dict; default = nothing, drop_missing = false, drop_nothing = false)
  """Return a converter that maps a dictionary of converters to a dictionary of values."""
  return function run_struct(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing
      return convertible
    end
    if default == "drop"
      converter_by_key = converters
    else
      converter_by_key = copy(converters)
      for key in keys(convertible.value)
        if !haskey(converter_by_key, key)
          converter_by_key[key] = default === nothing ? fail(N_("Unexpected item.")) : default
        end
      end
    end
    # TODO Add a {T} to the function to parametrize the keys of the 2 below dicts?
    converted_value_by_key = Dict()
    error_by_key = Dict()
    for (key, converter) in converter_by_key
      converted = converter(Convertible(get(convertible.value, key, nothing), convertible.context))
      if converted.value !== nothing || (!drop_nothing && (!drop_missing || key in convertible.value))
        converted_value_by_key[key] = converted.value
      end
      if converted.error !== nothing
        error_by_key[key] = converted.error
      end
    end
    return Convertible(converted_value_by_key, convertible.context, isempty(error_by_key) ? nothing : error_by_key)
  end
end

function struct(converters::Tuple; default = nothing)
  """Return a converter that maps a tuple of converters to a tuple (or array) of values."""
  return function run_struct(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing
      return convertible
    end
    values = convertible.value
    if default == "drop"
      values_converter = converters
    else
      values_converter = converters
      while length(values) > length(values_converter)
        values_converter = tuple(
          values_converter...,
          default === nothing ? fail(N_("Unexpected item.")) : default,
        )
      end
    end
    while length(values) < length(values_converter)
      values = tuple(values..., nothing)
    end
    converted_values = cell(length(values_converter))
    error_by_index = Dict{Int,Any}()
    for (index, (converter, value)) in enumerate(zip(values_converter, values))
      converted = converter(Convertible(value, convertible.context))
      converted_values[index] = converted.value
      if converted.error !== nothing
        error_by_index[index] = converted.error
      end
    end
    return Convertible(
      tuple(converted_values...),
      convertible.context,
      isempty(error_by_index) ? nothing : error_by_index,
    )
  end
end

function test(func; error = nothing, handle_nothing = false)
  """Return a converter that applies a test function to a value and returns an error when test fails.

  ``test`` always returns the initial value, even when test fails.

   See :doc:`how-to-create-converter` for more information.
  """
  return function run_test(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing && !handle_nothing || func(convertible.value)
      return convertible
    end
    return Convertible(
      convertible.value,
      convertible.context,
      error === nothing ? N_("Test failed.") : error,
    )
  end
end

function test_between(min_value, max_value; error = nothing)
  """Return a converter that accepts only values between the two given bounds (included).

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  return test(
    value -> min_value <= value <= max_value,
    error = error === nothing ? context -> _(context, "Value must be between $min_value and $max_value.") : error,
  )
end

function test_equal(comparison_value; error = nothing, handle_nothing = false)
  """Return a converter that accepts only values equals to given value.

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  return test(
    value -> value == comparison_value,
    error = error === nothing ? context -> _(context, "Value must be equal to $comparison_value.") : error,
    handle_nothing = handle_nothing,
  )
end

function test_greater_or_equal(min_value; error = nothing)
  """Return a converter that accepts only values greater than or equal to given value.

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  return test(
    value -> min_value <= value,
    error = error === nothing ? context -> _(context, "Value must be greater than or equal to $min_value.") : error,
  )
end

function test_in(values::Dict; error = nothing, handle_nothing = false)
  """Return a converter that accepts only values belonging to the keys of a given dictionary.

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  values_key = keys(values)
  return test(
    value -> value in values_key,
    error = error === nothing ? context -> _(
      context,
      string(
        "Value must belong to ",
        length(values_key) > 5 ? string(join(collect(values_key)[1:5], ", "), "..") : values_key,
        ".",
      ),
    ) : error,
    handle_nothing = handle_nothing,
  )
end

function test_in(values::Set; error = nothing, handle_nothing = false)
  """Return a converter that accepts only values belonging to a given set.

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  return test(
    value -> value in values,
    error = error === nothing ? context -> _(
      context,
      string(
        "Value must belong to ",
        length(values) > 5 ? string(join(collect(values)[1:5], ", "), "..") : collect(values),
        ".",
      ),
    ) : error,
    handle_nothing = handle_nothing,
  )
end

function test_in(values; error = nothing, handle_nothing = false)
  """Return a converter that accepts only values belonging to a given set (or array or...).

  .. warning:: Like most converters, a ``nothing`` value is not compared.
  """
  return test(
    value -> value in values,
    error = error === nothing ? context -> _(
      context,
      string(
        "Value must belong to ",
        length(values) > 5 ? string(join(values[1:5], ", "), "..") : values,
        ".",
      ),
    ) : error,
    handle_nothing = handle_nothing,
  )
end

function test_isa(data_type::Union{DataType, Union}; error = nothing, handle_nothing = false)
  """Return a converter that accepts only an instance of given type."""
  return test(
    value -> isa(value, data_type),
    error = error === nothing ? context -> _(context, "Value must be an instance of $data_type.") : error,
    handle_nothing = handle_nothing,
  )
end

function test_nothing(; error = nothing)
    """Return a converters that signals an error when value is not ``None``."""
  return function run_test_nothing(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing
      return convertible
    end
    error = error === nothing ? context -> _(context, "Unexpected value") : error
    return Convertible(convertible.value, convertible.context, error)
  end
end

test_nothing(convertible::Convertible) = test_nothing()(convertible)

function to_bool(convertible::Convertible)
  """Convert a julia data to a boolean.

  .. note:: For a converter that doesn't require a clean string, see :func:`input_to_bool`.

  .. note:: For a converter that accepts special strings like "f", "off", "no", etc, see :func:`guess_bool`.

  .. warning:: Like most converters, a ``nothing`` value is not converted.

    When you want ``nothing`` to be converted to ``false``, use::

      pipe(str_to_bool, default(False))
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  return to_bool(convertible.value, convertible.context)
end

function to_bool(value, context::Context)
  try
    return Convertible(convert(Int, value) != 0, context)
  catch
    return Convertible(value, context, N_("Value must be a boolean."))
  end
end

function to_bool(value::AbstractString, context::Context)
  try
    return Convertible(parse(Int, value) != 0, context)
  catch
    return Convertible(value, context, N_("Value must be a boolean."))
  end
end

function to_float(convertible::Convertible; accept_expression = false)
  """Convert a Julia data to a float number.

  .. warning:: Like most converters, a ``nothing`` value is not converted.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  value = convertible.value
  if accept_expression && isa(value, AbstractString)
    value = strip(value)
    if !ismatch(r"^[ \t\n\r\d.+\-*/()]+$", value)
      return Convertible(value, convertible.context, N_("Value must be a valid floating point expression."))
    end
    try
      value = eval(parse(value))
    catch
      return Convertible(value, convertible.context, N_("Value must be a valid floating point expression."))
    end
  end
  return to_float(value, convertible.context)
end

function to_float(; accept_expression = false)
  return convertible::Convertible -> to_float(convertible; accept_expression = accept_expression)
end

function to_float(value::AbstractString, context::Context)
  try
    return Convertible(float(value), context)
  catch
    return Convertible(value, context, N_("Value must be a float number."))
  end
end

function to_float(value, context::Context)
  try
    return Convertible(convert(Float64, value), context)
  catch
    return Convertible(value, context, N_("Value must be a float number."))
  end
end

function to_int(convertible::Convertible; accept_expression = false)
  """Convert a Julia data to an integer number.

  .. warning:: Like most converters, a ``nothing`` value is not converted.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  value = convertible.value
  if accept_expression && isa(value, AbstractString)
    value = strip(value)
    if !ismatch(r"^[ \t\n\r\d.+\-*/()]+$", value)
      return Convertible(value, convertible.context, N_("Value must be a valid integer expression."))
    end
    try
      value = eval(parse(value))
    catch
      return Convertible(value, convertible.context, N_("Value must be a valid integer expression."))
    end
  end
  return to_int(value, convertible.context)
end

function to_int(; accept_expression = false)
  return convertible::Convertible -> to_int(convertible; accept_expression = accept_expression)
end

function to_int(value::AbstractString, context::Context)
  try
    return Convertible(parse(Int, value), context)
  catch
    return Convertible(value, context, N_("Value must be an integer number."))
  end
end

function to_int(value, context::Context)
  try
    return Convertible(convert(Int, value), context)
  catch
    return Convertible(value, context, N_("Value must be an integer number."))
  end
end

function to_string(convertible::Convertible)
  """Convert a Julia data to a string.

  .. warning:: Like most converters, a ``nothing`` value is not converted.
  """
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  return Convertible(string(convertible.value), convertible.context)
end

function to_value(convertible::Convertible)
  """Check a conversion and either return its value or raise an *ErrorException* exception."""
  if convertible.error !== nothing
    error(eval_error(convertible.context, convertible.error), "\nValue: ", convertible.value)
  end
  return convertible.value
end

to_value_error(convertible::Convertible) = (convertible.value, eval_error(convertible.context, convertible.error))

function uniform_mapping(key_converter::Function, value_converters::Function...; drop_nothing = false,
    drop_nothing_keys = false, key_type = nothing, value_type = nothing)
  """Return a converter that applies a unique converter to each key and another one to each value of a mapping."""
  # TODO: Handle constructor.
  return function run_uniform_mapping(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing
      return convertible
    end
    error_by_key = Dict{Any, Any}()
    value_by_key = OrderedDict{Any, Any}()
    for (key, value) in convertible.value
      key_converted = key_converter(Convertible(key, convertible.context))
      if key_converted.error !== nothing
        error_by_key[key] = key_converted.error
      end
      if drop_nothing_keys && key_converted.value === nothing
        continue
      end
      value_converted = pipe(value_converters...)(Convertible(value, convertible.context))
      if !drop_nothing || value_converted.value !== nothing
        value_by_key[key_converted.value] = value_converted.value
      end
      if value_converted.error !== nothing
        error_by_key[key_converted.value] = value_converted.error
      end
    end

    typed_value_by_key = OrderedDict{
      key_type === nothing ?
        (isempty(value_by_key) ? Any : mapreduce(typeof, promote_type, Void, keys(value_by_key))) :
        key_type,
      value_type === nothing ?
        (isempty(value_by_key) ? Any : mapreduce(typeof, promote_type, Void, values(value_by_key))) :
        value_type,
    }()
    for (key, value) in value_by_key
      typed_value_by_key[key] = value
    end

    if isempty(error_by_key)
      return Convertible(typed_value_by_key, convertible.context)
    else
      typed_error_by_key = OrderedDict{
        key_type === nothing ? mapreduce(typeof, promote_type, Void, keys(error_by_key)) : key_type,
        value_type === nothing ? mapreduce(typeof, promote_type, Void, values(error_by_key)) : value_type,
      }()
      for (key, error) in error_by_key
        typed_error_by_key[key] = error
      end
      return Convertible(typed_value_by_key, convertible.context, typed_error_by_key)
    end
  end
end

function uniform_sequence(converters::Function...; drop_nothing = false, item_type = nothing, sequence_type = Array)
  """Return a converter that applies the same converter to each value of an array."""
  # TODO: Handle constructor or sequence_type.
  return function run_uniform_sequence(convertible::Convertible)
    if convertible.error !== nothing || convertible.value === nothing
      return convertible
    end
    error_by_index = Dict{Any, Any}()
    values = Any[]
    for (index, value) in enumerate(convertible.value)
      converted = pipe(converters...)(Convertible(value, convertible.context))
      if !drop_nothing || converted.value !== nothing
        push!(values, converted.value)
      end
      if converted.error !== nothing
        error_by_index[index] = converted.error
      end
    end
    return Convertible(
      collect(item_type === nothing ?
        (isempty(values) ? Any : mapreduce(typeof, promote_type, Void, values)) :
        item_type, values),
      convertible.context,
      isempty(error_by_index) ? nothing : error_by_index,
    )
  end
end


text:
immutable NURBSBasis <: Basis1D
    bs::BSplineBasis
    weights::Vector{Float64}
    deriv::Int

    function NURBSBasis(bs::BSplineBasis, weights, deriv=0)
        @assert_ex(bs.deriv == 0,
                   ArgumentError("Underlying B-Spline basis must not be differentiated"))
        @assert_ex(length(bs) == length(weights),
                   ArgumentError("Number of weights must be equal to number of basis functions"))
        @assert_ex(minimum(weights) > 0.0, ArgumentError("Weights must be positive"))
        @assert_ex(deriv <= min(nderivs(bs), 2), ArgumentError("Differentiation order not supported"))
        new(bs, weights, deriv)
    end

    NURBSBasis(b::BSplineBasis, deriv::Int) = NURBSBasis(b, ones(length(b)), deriv)
    NURBSBasis(b::BSplineBasis) = NURBSBasis(b, ones(length(b)))
end

typealias NURBS BasisFunction1D{NURBSBasis}

# Inherit functions from B-Spline bases
Base.length(b::NURBSBasis) = length(b.bs)
Base.size(b::NURBSBasis) = size(b.bs)
Base.getindex(b::NURBSBasis, i) = NURBS(b, i, b.deriv)

nderivs(b::NURBSBasis) = min(2, nderivs(b.bs))

domain(b::NURBSBasis) = domain(b.bs)
domain(b::NURBS) = domain(b.basis[b.index])

order(b::NURBSBasis) = throw(ArgumentError("NURBS are not polynomial"))
order(b::NURBS) = throw(ArgumentError("NURBS are not polynomial"))

deriv(b::NURBSBasis) = NURBSBasis(b.bs, b.weights, b.deriv + 1)
deriv(b::NURBSBasis, order) = NURBSBasis(b.bs, b.weights, b.deriv + order)

supported(b::NURBSBasis, pts) = supported(b.bs, pts)

function evaluate_raw{T<:Real}(b::NURBSBasis, pts::Vector{T}, deriv::Int, rng::UnitRange{Int})
    const bwts = b.weights[rng]

    bvals = evaluate_raw(b.bs, pts, 0, rng) .* bwts
    wts = sum(bvals, 1)

    if deriv == 0
        return bvals ./ wts
    end

    bvals1 = evaluate_raw(b.bs, pts, 1, rng) .* bwts
    wts1 = sum(bvals1, 1)
    d1 = (bvals1 .* wts - bvals .* wts1) ./ (wts .^ 2)

    if deriv == 1
        return d1
    end

    bvals2 = evaluate_raw(b.bs, pts, 2, rng) .* bwts
    wts2 = sum(bvals2, 1)
    d2 = (bvals2 .* wts - bvals .* wts2) ./ (wts .^ 2)

    if deriv == 2
        return d2 - 2 * d1 .* wts1 ./ wts
    end

    throw(ArgumentError("Third order derivatives or higher are not supported by NURBS bases"))
end

ext: jl


text:
import ClassImbalance
import DataFrames
import StatsBase
import Test

df = DataFrames.DataFrame()
num_rows = 1_000
df[:x1] = randn(num_rows)
df[:x2] = StatsBase.sample(["a", "b"], num_rows)
df[:x3] = randn(num_rows)
df[:x4] = StatsBase.sample(["d", "e", "f"], num_rows)
df[:x5] = randn(num_rows)
df[:x6] = StatsBase.sample(["g", "h", "i", "j"], num_rows)
df[:y] = StatsBase.sample([0,1], num_rows)

df_smote_obs_result = ClassImbalance.smote_obs(df,0.01,0)
df_smote_obs_result = ClassImbalance.smote_obs(df,300,0)

df_1 = DataFrames.DataFrame(x = ["a", "b", "c"])
df_2 = ClassImbalance.dataframe_to_matrix(df_1, [1], 3, 1)
df_3 = ClassImbalance.matrix_to_dataframe(df_2, df_1, [1])
Test.@test(all(df_1[:, 1] .== df_3[:, 1]))

num_rows = 1_000
x1 = randn(num_rows)
x2 = randn(num_rows)
x3 = randn(num_rows)
X = hcat(x1, x2, x3)
X_smote_obs_result = ClassImbalance.smote_obs(X, 0.01, 0)
X_smote_obs_result = ClassImbalance.smote_obs(X, 300, 0)

w1 = [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Test.@test(
    ClassImbalance.cases_needed(w1) == 3
    )
Test.@test(
    ClassImbalance.pct_needed(w1) == 100.0
    )


text:
@testset "Array interface" begin

    # Test aliases
    @test QuantizedArray(rand(10), k=2, method=:sample) isa QuantizedVector
    @test QuantizedArray(rand(2, 10), k=2, method=:sample) isa QuantizedMatrix

    # Test outer constructor checks
    @test_throws AssertionError QuantizedArray(rand(2, 2, 2), k=2, m=1, method=:sample)
    @test_throws AssertionError QuantizedArray(rand(2, 10), k=-1, m=1, method=:sample)
    @test_throws AssertionError QuantizedArray(rand(2, 10), k=2, m=-1, method=:sample)
    @test_throws AssertionError QuantizedArray(rand(3, 10), k=2, m=2, method=:sample)
    for method in [:pq, :opq, :rvq]
        @test_throws AssertionError QuantizedArray(rand([1,2,3], 3, 10), k=2, m=2, method=method)
    end

    # Test quantize function
    T = Float32
    X = rand(T, 3, 10)
    qa = QuantizedArray(X, k=2, m=3, method=:sample)
    @test QuantizedArrays.quantizer(qa) === qa.quantizer
    @test quantize(QuantizedArrays.quantizer(qa), X) == qa

    # Test other interface parts
    @test size(qa) == qa.quantizer.dims
    @test eltype(qa) == T
    @test QuantizedArrays.nvars(rand(10)) == 1
    @test QuantizedArrays.nvars(rand(2, 10)) == 2

    # Test setindex!
    @test_throws ErrorException qa[1] = one(T)
    @test_throws ErrorException qa[1,1] = one(T)

    # Test orthogonal reconstruction: getindex

    # vector
    truevector = [0, 1, 0, 2, 0]
    codes = [0x00, 0x01, 0x02]
    vectors = [0 1 2]
    data = [0x00 0x01 0x00 0x02 0x00]
    cb = CodeBook(codes, vectors)
    q = QuantizedArrays.OrthogonalQuantization()
    rot = diagm(0=>ones(eltype(truevector), length(truevector)))
    aq = ArrayQuantizer(q, size(truevector), [cb], length(codes), Distances.SqEuclidean(), rot)
    qa = QuantizedArray(aq, data)
    @test qa[:] == truevector
    @test all(qa .== truevector)

    # matrix
    truematrix = [0 1 0 2 0 3;
                  0 1 0 2 0 3]
    codes = [0x00, 0x01, 0x02, 0x03]
    vectors = [0 1 2 3;
               0 1 2 3]
    data = [0x00 0x01 0x00 0x02 0x00 0x03]
    cb = CodeBook(codes, vectors)
    q = QuantizedArrays.OrthogonalQuantization()
    rot = diagm(0=>ones(eltype(truematrix), size(truematrix, 1)))
    aq = ArrayQuantizer(q, size(truematrix), [cb], length(codes), Distances.SqEuclidean(), rot)
    qa = QuantizedArray(aq, data)
    @test qa[:] == truematrix[:]
    @test qa[:,:] == truematrix
    @test all(qa .== truematrix)

    # Test aditive reconstruction: getindex

    # vector
    truevector = [0, 1, 0, 2, 0] .+ [-1, -1, 1, 1, 2]
    codes = [0x00, 0x01, 0x02]
    vectors1 = [0 1 2]
    vectors2 = [-1 1 2]
    data = [0x00 0x01 0x00 0x02 0x00;
            0x00 0x00 0x01 0x01 0x02]
    cbs = [CodeBook(codes, vectors1), CodeBook(codes, vectors2)]
    q = QuantizedArrays.AdditiveQuantization()
    rot = diagm(0=>ones(eltype(truevector), length(truevector)))
    aq = ArrayQuantizer(q, size(truevector), cbs, length(codes), Distances.SqEuclidean(), rot)
    qa = QuantizedArray(aq, data)
    @test qa[:] == truevector
    @test all(qa .== truevector)

    # matrix
    truematrix = [0 1 0 2 0 3;
                  0 1 0 2 0 3] .+
                 [ 2  3  1  4  4  1;
                  -2 -3 -1 -4 -4 -1]
    codes = [0x00, 0x01, 0x02, 0x03]
    vectors1 = [0 1 2 3;
                0 1 2 3]
    vectors2 = [1  2  3  4
                -1 -2 -3 -4]
    data = [0x00 0x01 0x00 0x02 0x00 0x03;
            0x01 0x02 0x00 0x03 0x03 0x00]
    cbs = [CodeBook(codes, vectors1), CodeBook(codes, vectors2)]
    q = QuantizedArrays.AdditiveQuantization()
    rot = diagm(0=>ones(eltype(truematrix), size(truematrix, 1)))
    aq = ArrayQuantizer(q, size(truematrix), cbs, length(codes), Distances.SqEuclidean(), rot)
    qa = QuantizedArray(aq, data)
    @test qa[:] == truematrix[:]
    @test qa[:,:] == truematrix
    @test all(qa .== truematrix)

    # Setindex tests
    @test_throws ErrorException qa[1] = 1
    @test_throws ErrorException qa[1, 1] = 1

end


text:
function formulation(
    i::SpanningTreeInstance{Int, Maximise},
    f::DefaultLinearFormulation
)
    # Input graph supposed to be undirected.
    n = nv(i.graph)
    graph = DiGraph(n)

    # Build the equivalent directed graph.
    for e in edges(i.graph)
        add_edge!(graph, src(e), dst(e))
        add_edge!(graph, dst(e), src(e))
    end

    # Helper methods.
    inedges(g, v) = (edgetype(g)(x, v) for x in inneighbors(g, v))
    outedges(g, v) = (edgetype(g)(v, x) for x in outneighbors(g, v))

    # Choose a source arbitrarily.
    source = first(vertices(graph))
    other_nodes = filter(v -> v != source, vertices(graph))

    # Start building the optimisation model. Keep integer variables, because
    # this formulation is not tight and, more importantly, there might be other
    # constraints later on that would destroy any interesting structure of the
    # problem.
    # Based on Magnanti, T.L.; Wolsey, L. Optimal Trees, section Flow
    # formulation (p. 38).
    m = Model(f.solver)
    set_silent(m)

    x = @variable(m, [e in edges(graph)], binary = true)
    flow =
        @variable(m, [e in edges(graph)], lower_bound = 0, upper_bound = n - 1)

    function inflow(v)
        return length(inedges(graph, v)) == 0 ? 0 :
               sum(flow[e] for e in inedges(graph, v))
    end
    function outflow(v)
        return length(outedges(graph, v)) == 0 ? 0 :
               sum(flow[e] for e in outedges(graph, v))
    end

    function inx(v)
        return length(inedges(graph, v)) == 0 ? 0 :
               sum(x[e] for e in inedges(graph, v))
    end
    function outx(v)
        return length(outedges(graph, v)) == 0 ? 0 :
               sum(x[e] for e in outedges(graph, v))
    end

    # Set the names, due to the use of anonymous variables with JuMP.
    for e in edges(graph)
        set_name(x[e], "x_$(src(e))_$(dst(e))")
        set_name(flow[e], "flow_$(src(e))_$(dst(e))")
    end

    # - Flow towards the source: zero. Flow from the source: one unit for each
    #   vertex in the graph.
    @constraint(m, inflow(source) == 0)
    @constraint(m, outflow(source) == n - 1)
    @constraint(m, inx(source) == 0)
    @constraint(m, outx(source) >= 1)
    @constraint(m, outx(source) <= 4)

    # - When the flow exits one of the other nodes, the node keeps one unit
    #   for itself and redistributes the rest.
    for v in other_nodes
        @constraint(m, inflow(v) - outflow(v) == 1)
    end

    # - (n - 1) edges must be taken for n vertices.
    @constraint(m, sum(x) == n - 1)

    # - Flow can flow in an edge only if the edge belongs to the tree.
    for e in edges(graph)
        @constraint(m, flow[e] <= (n - 1) * x[e])
    end

    # - Among direct and reverse edges, only one can be chosen.
    for e in edges(graph)
        @constraint(m, x[e] + x[reverse(e)] <= 1)
    end

    # Finally, the objective function.
    @objective(
        m,
        Max,
        sum(i.rewards[e] * (x[e] + x[reverse(e)]) for e in keys(i.rewards))
    )

    # Don't return flows, these are not really decision variables.
    return m, x
end

function solve(
    i::SpanningTreeInstance{Int, Maximise},
    f::DefaultLinearFormulation
)
    m, x = formulation(i, f)
    optimize!(m)
    return SpanningTreeSolution(i, _extract_lp_solution(i, x))
end


text:
iv = Interval(0,1)
@test minimum(iv) == 0
@test maximum(iv) == 1
@test 0.5 in iv
@test !(-0.5 in iv)
@test (1.0 ∈ iv)

iv = Interval(-0.2,1)
@test minimum(iv) == -0.2
@test maximum(iv) == 1
@test -0.2 in iv
@test !(1.5 in iv)
@test (1.0 ∈ iv)

@test issubset(Interval(0,1),Interval(0,2))
@test !issubset(Interval(0,1),Interval(0.5,2))
@test !issubset(Interval(0,2),Interval(0.5,1))
@test !issubset(Interval(0,1),Interval(-0.5,0.5))
@test Interval(0,1) ⊆ Interval(0,2)
@test !(Interval(0,1) ⊆ Interval(0.5,2))
@test Interval(0,1) ⊈ Interval(0.5,2)
@test !(Interval(0.5,1) ⊈ Interval(0,2))

@test isempty(Interval(1,-1))
@test isempty(Interval(1,-0.5))
@test ~isempty(Interval(1,1))

@test intersect(Interval(1,-0.5), Interval(1,2)) == Interval(1,-1)
@test intersect(Interval(0,1), Interval(1,2)) == Interval(1,1)
@test intersect(Interval(0,1), Interval(0.5,2)) == Interval(0.5,1)
@test (Interval(0,1) ∩ Interval(-0.5,2)) == Interval(0,1)
@test (Interval(1,1) ∩ Interval(-0.5,2)) == Interval(1,1)
@test (Interval(1.1,1) ∩ Interval(1.05,2)) == Interval(1,-1)
@test (Interval(1.1,1) ∩ Interval(1.05,2)) == TDD.∅
@test (Interval(1,2) ∩ Interval(1.05,1.5) ∩ Interval(0.7,1.3)) == Interval(1.05,1.3)


text:
using MinFEM

mesh = import_mesh("../meshes/poisson.msh")

L = asmElasticity(mesh,1.0,1.0)
Mb = asmBoundaryMassMatrix(mesh, Edges=mesh.Boundaries[1001].Edges, qdim=2)

f(x) = [0.1,0.0]
s = evaluateMeshFunction(mesh, f, region=mesh.Boundaries[1001].Nodes, qdim=2)

boundary = mesh.Boundaries[1003].Nodes

pde = PDESystem(A=L, b=Mb*s, bc=zeros(2*mesh.nnodes), DI=boundary, qdim=2)
solve(pde)

vtkfile = open_vtk_file(mesh, "elasticity.vtu")
write_point_data(vtkfile, [reshape(pde.state, 2, mesh.nnodes); zeros(1,mesh.nnodes)], "Y")
write_point_data(vtkfile, s, "S")
save_vtk_file(vtkfile)


text:
ArrayAbstract{N,T} = AbstractArray{T,N}

"""
Assumed PoU is Boolean.

"""
makecoherent

"""
    makecoherent!(x::CartesianDDMVector)

Make coherent.

"""
function makecoherent!(x::CartesianDDMVector)
    (; context, parent) = x
    nocontext = removeoverlap(context)

    dec, _ = ranges(context)
    part, _ = ranges(nocontext)

    makecoherent!(dec, part, parent)

    x
end

"""
Make coherent in 1d

"""
#=
makecoherent!(dec::ArrayAbstract{1},
              part::ArrayAbstract{1},
              parent::ArrayAbstract{1}) =
    _makecoherent!(dec, part, parent)
=#

"""
Make coherent across faces

"""
function makecoherent!(dec::ArrayAbstract{1},
                       part::ArrayAbstract{1},
                       parent::ArrayAbstract{1})
    for i in drop(axes(parent, 1), 1)
        _makecoherent!((parent[i-1], dec[i-1]),
                       (parent[i], dec[i], part[i]))
        _makecoherent!((parent[i], dec[i]),
                       (parent[i-1], dec[i-1], part[i-1]))
    end
end

"""
Make coherent in 2d

"""
function makecoherent!(dec::ArrayAbstract{2},
                       part::ArrayAbstract{2},
                       parent::ArrayAbstract{2})
    for j in axes(parent, 2)
        makecoherent!(@view(dec[:, j]),
                      @view(part[:, j]),
                      @view(parent[:, j]))
    end
    for i in axes(parent, 1)
        makecoherent!(@view(dec[i, :]),
                      @view(part[i, :]),
                      @view(parent[i, :]))
    end

    _makecoherent!(dec, part, parent)
end

"""
Make coherent across edges

"""
function _makecoherent!(dec::ArrayAbstract{2},
                        part::ArrayAbstract{2},
                        parent::ArrayAbstract{2})
    for j in drop(axes(parent, 2), 1)
        for i in drop(axes(parent, 1), 1)
            _makecoherent!((parent[i-1, j-1], dec[i-1, j-1]),
                           (parent[i, j], dec[i, j], part[i, j]))
            _makecoherent!((parent[i, j], dec[i, j]),
                           (parent[i-1, j-1], dec[i-1, j-1], part[i-1, j-1]))
            _makecoherent!((parent[i, j-1], dec[i, j-1]),
                           (parent[i-1, j], dec[i-1, j], part[i-1, j]))
            _makecoherent!((parent[i-1, j], dec[i-1, j]),
                           (parent[i, j-1], dec[i, j-1], part[i, j-1]))
        end
    end
end

"""
Make coherent in 3d

"""
function makecoherent!(dec::ArrayAbstract{3},
                       part::ArrayAbstract{3},
                       parent::ArrayAbstract{3})
    for k in axes(parent, 3)
        for j in axes(parent, 2)
            makecoherent!(@view(dec[:, j, k]),
                          @view(part[:, j, k]),
                          @view(parent[:, j, k]))
        end
    end
    for k in axes(parent, 3)
        for i in axes(parent, 1)
            makecoherent!(@view(dec[i, :, k]),
                          @view(part[i, :, k]),
                          @view(parent[i, :, k]))
        end
    end
    for j in axes(parent, 2)
        for i in axes(parent, 1)
            makecoherent!(@view(dec[i, j, :]),
                          @view(part[i, j, :]),
                          @view(parent[i, j, :]))
        end
    end

    for k in axes(parent, 3)
        makecoherent!(@view(dec[:, :, k]),
                      @view(part[:, :, k]),
                      @view(parent[:, :, k]))
    end
    for j in axes(parent, 2)
        makecoherent!(@view(dec[:, j, :]),
                      @view(part[:, j, :]),
                      @view(parent[:, j, :]))
    end
    for i in axes(parent, 1)
        makecoherent!(@view(dec[i, :, :]),
                      @view(part[i, :, :]),
                      @view(parent[i, :, :]))
    end

    _makecoherent!(dec, part, parent)
end

"""
Make coherent across vertices

"""
function _makecoherent!(dec::ArrayAbstract{3},
                        part::ArrayAbstract{3},
                        parent::ArrayAbstract{3})
    for k in drop(axes(parent, 3), 1)
        for j in drop(axes(parent, 2), 1)
            for i in drop(axes(parent, 1), 1)
                _makecoherent!((parent[i-1, j-1, k-1], dec[i-1, j-1, k-1]),
                               (parent[i, j, k], dec[i, j, k], part[i, j, k]))
                _makecoherent!((parent[i, j, k], dec[i, j, k]),
                               (parent[i-1, j-1, k-1], dec[i-1, j-1, k-1], part[i-1, j-1, k-1]))
                _makecoherent!((parent[i, j, k-1], dec[i, j, k-1]),
                               (parent[i-1, j-1, k], dec[i-1, j-1, k], part[i-1, j-1, k]))
                _makecoherent!((parent[i-1, j-1, k], dec[i-1, j-1, k]),
                               (parent[i, j, k-1], dec[i, j, k-1], part[i, j, k-1]))
                _makecoherent!((parent[i-1, j, k-1], dec[i-1, j, k-1]),
                               (parent[i, j-1, k], dec[i, j-1, k], part[i, j-1, k]))
                _makecoherent!((parent[i, j-1, k], dec[i, j-1, k]),
                               (parent[i-1, j, k-1], dec[i-1, j, k-1], part[i-1, j, k-1]))
                _makecoherent!((parent[i, j-1, k-1], dec[i, j-1, k-1]),
                               (parent[i-1, j, k], dec[i-1, j, k], part[i-1, j, k]))
                _makecoherent!((parent[i-1, j, k], dec[i-1, j, k]),
                               (parent[i, j-1, k-1], dec[i, j-1, k-1], part[i, j-1, k-1]))
            end
        end
    end
end

"""

"""
function _makecoherent!((out, gblout), (in, gblin, nogblin))
    rout = reshape(out, gblout...)
    rin = reshape(in, gblin...)

    iter = intersect(CartesianIndices.((gblout, nogblin))...)

    rout[iter] .= rin[iter]
end


text:
function maximum_weight_maximal_matching_lp(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
    return maximum_weight_maximal_matching_lp(g, solver, cutoff_weights(w, cutoff))
end

function maximum_weight_maximal_matching_lp(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}) where {T<:Real}
# TODO support for graphs with zero degree nodes
# TODO apply separately on each connected component
    bpmap = bipartite_map(g)
    length(bpmap) != nv(g) && error("Graph is not bipartite")
    v1 = findall(isequal(1), bpmap)
    v2 = findall(isequal(2), bpmap)
    if length(v1) > length(v2)
        v1, v2 = v2, v1
    end

    nedg = 0
    edgemap = Dict{Edge,Int}()

    for j in 1:size(w,2)
        for i in 1:size(w,1)
            if w[i,j] > 0.0
                nedg += 1
                edgemap[Edge(i,j)] = nedg
                edgemap[Edge(j,i)] = nedg
            end
        end
    end

    model = Model(solver=solver)
    @variable(model, x[1:length(w)] >= 0)

    for i in v1
        idx = Vector{Int}()
        for j in neighbors(g, i)
            if haskey(edgemap, Edge(i,j))
                push!(idx, edgemap[Edge(i,j)])
            end
        end
        if length(idx) > 0
            @constraint(model, sum(x[id] for id = idx) == 1)
        end
    end

    for j in v2
        idx = Vector{Int}()
        for i in neighbors(g, j)
            if haskey(edgemap, Edge(i,j))
                push!(idx, edgemap[Edge(i,j)])
            end
        end

        if length(idx) > 0
            @constraint(model, sum(x[id] for id = idx) <= 1)
        end
    end

    @objective(model, Max, sum(w[src(e),dst(e)] * x[edgemap[e]] for e in keys(edgemap)))

    status = solve(model)
    status != :Optimal && error("JuMP solver failed to find optimal solution.")
    sol = getvalue(x)

    all(Bool[s == 1 || s == 0 for s in sol]) || error("Found non-integer solution.")

    cost = getobjectivevalue(model)

    mate = fill(-1, nv(g))
    for e in edges(g)
        if w[src(e),dst(e)] > zero(T)
            inmatch = convert(Bool, sol[edgemap[e]])
            if inmatch
                mate[src(e)] = dst(e)
                mate[dst(e)] = src(e)
            end
        end
    end

    return MatchingResult(cost, mate)
end


text:
@testset "string" begin

  @test string(Double64(pi)) == "3.1415926535897932384626433832795058"
  @test string(Double32(pi)) == "3.141592653589796"
  @test string(Double16(pi)) == "3.1415925"

  @test string(ComplexDF64(1.0,1.0)) == "1.0 + 1.0im"
  @test string(ComplexDF32(1.0,1.0)) == "1.0 + 1.0im"
  @test string(ComplexDF16(1.0,1.0)) == "1.0 + 1.0im"

  @test stringtyped(Double64(pi)) == "Double64(3.141592653589793, 1.2246467991473532e-16)"
#  @test stringtyped(Double32(pi)) == "Double32(3.1415927, -8.742278e-8)"  julia version matters
  @test stringtyped(Double16(pi)) == "Double16(3.14, 0.0009675)"

  @test stringtyped(ComplexDF64(1.0,1.0)) == "ComplexDF64(Double64(1.0, 0.0), Double64(1.0, 0.0))"
  @test stringtyped(ComplexDF32(1.0,1.0)) == "ComplexDF32(Double32(1.0, 0.0), Double32(1.0, 0.0))"
  @test stringtyped(ComplexDF16(1.0,1.0)) == "ComplexDF16(Double16(1.0, 0.0), Double16(1.0, 0.0))"

end

@testset "sprintf" begin

  @test @sprintf("%f", Double64(3)) == "3.000000"

end


text:
export PermutationGraph, RandomPermutationGraph
export PermutationRepresentation

"""
`perm_adj(a,b)` is used to test for adajency in permutation graphs.
This function is not exported.
"""
function perm_adj(a::Tuple{S,T}, b::Tuple{S,T}) where {S<:Real,T<:Real}
    if a[1] < b[1]
        return a[2] > b[2]
    end
    return a[2] < b[2]
end

"""
`PermutationGraph(p::Permutation,q::Permutation)` builds a permutation
graph based on the permutations `p` and `q`. The vertex set is `1:n`
where `n` is the length of the permutations.
"""
function PermutationGraph(p::Permutation, q::Permutation)
    n = length(p)
    if length(q) != n
        error("Two permutations must be the same length")
    end
    G = IntGraph(n)
    for u = 1:n-1
        for v = u+1:n
            if perm_adj((p[u], q[u]), (p[v], q[v]))
                add!(G, u, v)
            end
        end
    end
    return G
end

"""
`PermutationGraph(p::Permutation)` is equivalent to
`PermutationGraph(p,id)` where `id` is the identity permutation.
"""
function PermutationGraph(p::Permutation)
    n = length(p)
    q = Permutation(n)
    return PermutationGraph(p, q)
end

"""
`PermutationGraph(d::Dict)` creates a permutation graph form a
dictionary that maps vertex names to pairs of real numbers.
The vertex type of the resulting graph is the key type of `d`.
The values in `d` must be all `Tuple{S,T}` where `S` and `T` are
subtypes of `Real`. For example, declare `d` like this:
`d = Dict{String, Tuple{Int,Int}}()`.
"""
function PermutationGraph(d::Dict{VV,Tuple{S,T}}) where {VV,S<:Real,T<:Real}
    vtcs = collect(keys(d))
    G = SimpleGraph{VV}()
    for v in vtcs
        add!(G, v)
    end

    n = length(vtcs)
    for i = 1:n-1
        u = vtcs[i]
        for j = i+1:n
            v = vtcs[j]
            if perm_adj(d[u], d[v])
                add!(G, u, v)
            end
        end
    end
    return G
end

"""
`PermutationGraph(f::Dict,g::Dict)` constructs a permutation graph
from a pair of mappings from a vertex set to real values.
"""
function PermutationGraph(f::Dict{T,R}, g::Dict{T,S}) where {T,R<:Real,S<:Real}
    # mush the two dictionaries into one
    h = Dict{T,Tuple{R,S}}()
    for k in keys(f)
        h[k] = f[k], g[k]
    end

    # invoke previous method
    G = PermutationGraph(h)
    cache_save(G, :PermutationRepresentation, (f, g))
    cache_save(G, :name, "Permutation graph")
    return G
end

"""
`RandomPermutationGraph(n)` creates a random permutation graph with
`n` vertices.
"""
function RandomPermutationGraph(n::Int)
    p = RandomPermutation(n)
    q = RandomPermutation(n)
    return PermutationGraph(p, q)
end

"""
`PermutationRepresentation(G)` returns a pair of dictionaries mapping
the vertices of `G` to integers. This pair of dictionaries form a
permutation representation of `G`. (If `G` is not a permutation graph,
an error is thrown.
"""
function PermutationRepresentation(G::SimpleGraph)
    if cache_check(G, :PermutationRepresentation)
        return cache_recall(G, :PermutationRepresentation)
    end
    A = SimpleDigraph()
    try
        A = transitive_orientation(complement(G))
    catch
        error("This graph does not have a permutation representation")
    end
    G1 = deepcopy(A)
    G2 = transitive_orientation(G)
    for e in elist(G2)
        add!(G1, e[1], e[2])
    end
    T = eltype(G)
    sigma = Dict{T,Int}()
    tau = Dict{T,Int}()
    vs1 = sort(vlist(G1), by = v -> out_deg(G1, v), rev = true)
    i = 1
    for v in vs1
        sigma[v] = i
        i = i + 1
    end
    G3 = deepcopy(A)
    for e in elist(G2)
        add!(G3, e[2], e[1])
    end
    vs2 = sort(vlist(G3), by = v -> out_deg(G3, v), rev = true)
    j = 1
    for v in vs2
        tau[v] = j
        j = j + 1
    end
    cache_save(G, :PermutationRepresentation, (sigma, tau))
    return sigma, tau
end


text:
"Splits the given string into arguments"
function split_compound_args(args::String)
    elems = Vector{String}()
    bracket_stack = Vector{Char}()
    last_index = 1
    ignore = false # for skipping inner parsing of string constants

    for ind = 1:length(args)
        if args[ind] == '\'' || args[ind] == '"'
            # if string constant is found, ignore any inner parsing until  you find closing '/"
            ignore = true
        elseif ignore && (args[ind] == '\'' || args[ind] == '"')
            # found closing '/", start splitting again
            ignore == false
        elseif args[ind] == '[' && !ignore
            # starting a list
            push!(bracket_stack, '[')
        elseif args[ind] == ']' && !ignore && ind != length(args)
            # found end of the list, not as at last character in the string
            elem = pop!(bracket_stack)
            if elem != '['
                error("wrongly formatted string $args")
            end
        elseif args[ind] == '(' && !ignore
            # start processing args of a structure
            push!(bracket_stack, '(')
        elseif args[ind] == ')' && !ignore && ind != length(args)
            # found end of the structure, not at the end of everything
            elem = pop!(bracket_stack)
            if elem != '('
                error("Wrongly formatted string $args")
            end
        elseif args[ind] == ',' && isempty(bracket_stack) && !ignore
            # found end of the argument
            push!(elems, args[last_index:ind-1])
            last_index = ind + 1
        elseif ind == length(args)
            # encountered end
            push!(elems, args[last_index:end])
        end
    end

    elems
end

"""
    converts a string to HerbLanguage

    [A-Z].* || _.*              -> Variable
    [a-z].* || '.*' || ".*"     -> Constant
    functor(args)               -> Structure if start_with_literals=false, otherwise Literal
    [...]                       -> List
    [X|Y]                       -> LPair
    \\+ ...                     -> Negation
"""
function parse_prolog(elem::String; start_with_literals=false)
    elem = String(strip(elem))
    if tryparse(Int,elem) != nothing
        tryparse(Int, elem)
    elseif tryparse(Float32, elem) != nothing
        tryparse(Float32, elem)
    elseif startswith(elem, "[")
        if occursin("|", elem)
            elem = elem[begin+1:end-1]
            separator_index = findfirst(isequal('|'), elem)
            head = elem[begin:separator_index-1]
            tail = elem[separator_index+1:end]
            LPair(parse_prolog(head), parse_prolog(tail))
        else
            args_to_parse = elem[begin+1:end-1]
            parsed_args = split_compound_args(args_to_parse)
            List([parse_prolog(a) for a in parsed_args])
        end
    elseif startswith(elem, "\\+")
        elem = strip(elem[begin+2:end])
        Negation(parse_prolog(elem, start_with_literals=true))
    elseif occursin("(", elem)
        open_bracket = findfirst(isequal('('), elem)
        closed_bracket = findlast(isequal(')'), elem)
        functor_name = elem[begin:open_bracket-1]
        args_to_parse = elem[open_bracket+1:closed_bracket-1]
        parsed_args = split_compound_args(args_to_parse)

        if start_with_literals
            Literal(c_pred(functor_name, length(parsed_args)), [parse_prolog(a) for a in parsed_args])
        else
            Structure(c_functor(functor_name, length(parsed_args)), [parse_prolog(a) for a in parsed_args])
        end
    elseif isuppercase(elem[1]) || startswith(elem, "_")
        c_var(elem)
    elseif (islowercase(elem[1]) || startswith(elem, "'") || startswith(elem, "\"")) && !start_with_literals
        c_const(elem)
    elseif islowercase(elem[1]) && start_with_literals
        c_prop(elem)
    else
        throw(DomainError(elem, "don't know how to convert it"))
    end
end

"""
    Parses a clause from a string
"""
function parse_clause(elem::String)
    sep = findfirst(":-", elem)
    head = elem[begin:sep[begin]-1]
    body = split_compound_args(elem[sep[begin+1]+1:end])
    Clause(parse_prolog(head, start_with_literals=true), Conj([parse_prolog(l, start_with_literals=true) for l in body]))
end

"""
    Parses a Prolog expression from a string
    asssumes a single clause

    base=true the first encounter of a compound structure is deemed Literal; otherwise Structure
              0-arity structures are interpreted as propositions; otherwise constants
"""
function from_string(elem::String; base=true)
    if occursin(":-", elem)
        parse_clause(elem)
    else
        parse_prolog(elem, start_with_literals=base)
    end
end

function macro_parse_term(expr)
    if isa(expr, Number)
        # numbers remain numbers
        return :($expr)
    elseif isa(expr, Symbol)
        if isuppercase(string(expr)[begin]) || string(expr)[1] == "_"
            # uppercased symbols are variables
            v = c_var(string(expr))
            return :($v)
        else
            # lowercase symbols are constants
            c = c_const(string(expr))
            return :($c)
        end
    elseif isa(expr, Expr)
        if expr.head == :call && expr.args[1] != :|
            # compounds term goes to structure
            functor_name = string(expr.args[begin])
            args = [macro_parse_term(e) for e in expr.args[2:end]]
            s = Structure(c_functor(functor_name, length(args)), args)
            return :($s)
        elseif expr.head == :call && expr.args[1] == :|
            # turns into a LPair
            head = macro_parse_term(expr.args[2])
            tail = macro_parse_term(expr.args[3])
            p = LPair(head, tail)
            return :($p)
        elseif expr.head == :vect
            # turn into list
            elems = [macro_parse_term(e) for e in expr.args]
            l = List(elems)
            return :($l)
        end
    end
end

function macro_parse_body(expr::Union{Symbol,Expr})
    if expr.args[begin] == :&
        # parse a conjunction
        left = macro_parse_body(expr.args[2])
        right = macro_parse_top(expr.args[3])
        push!(left, right)
        return :($left)
    else
        # parse single literal
        container = Vector{Union{Literal, Negation, Proposition}}()
        push!(container, macro_parse_top(expr))
        return container
    end
end

function macro_parse_top(expr::Union{Symbol,Expr})
    if isa(expr, Symbol)
        # a symbol on a top level turns into a propositions
        prop = c_prop(string(expr))
        return :($prop)
    elseif expr.head == :vect
        # list of clauses, turn into a Program
        clauses = [macro_parse_top(e) for e in expr.args]
        prog = Program(clauses)
        return :($prog)
    elseif expr.head == :call
        if expr.args[begin] == :<=
            #parse clause
            head = macro_parse_top(expr.args[2])
            body = macro_parse_body(expr.args[3])
            cl = Clause(head, Conj(body))
            return :($cl)
        elseif expr.args[begin] == :! && length(expr.args) == 2
            # parse negation
            lit = macro_parse_top(expr.args[2])
            neg = Negation(lit)
            return :($neg)
        else
            # parse literal
            pred_name = string(expr.args[begin])
            args = [macro_parse_term(e) for e in expr.args[2:end]]
            lit = Literal(c_pred(pred_name, length(args)), args)
            return :($lit)
        end
    end
end

macro prolog(expr)
    return macro_parse_top(expr)
end

function parse_disjunction(expr)
    if expr.head == :(::)
        #parse single literal
        container = Vector{ProbabilisticLiteral}()
        push!(container, macro_probabilistic_parse_top(expr))
        return container
    else
        #parse disjunction
        left = parse_disjunction(expr.args[2])
        right = macro_probabilistic_parse_top(expr.args[3])
        push!(left, right)
        return :($left)
    end
end

function isdisjunction(expr)
    if expr.args[begin] == :|
        true
    else
        false
    end
end

function macro_probabilistic_parse_top(expr::Union{Expr,Symbol})
    if expr.head == :(::)
        #parse probabilistic literal
        prob = expr.args[begin]
        literal = macro_parse_top(expr.args[2])
        plit = ProbabilisticLiteral(prob,literal)
        return :($plit)
    elseif expr.head == :call
        if isdisjunction(expr)
            heads = parse_disjunction(expr)
            ad = AnnotatedDisjunction(heads, true)
            return :($ad)
        elseif expr.args[begin] == :<= && isdisjunction(expr.args[2])
            heads = parse_disjunction(expr.args[2])

            if isa(expr.args[3], Union{Expr, Symbol})
                # if body is a conjunction
                body = macro_parse_top(expr.args[3])
                ad = AnnotatedDisjunction(heads, body)
                return :($ad)
            else
                # if body is true
                ad = AnnotatedDisjunction(heads, expr.args[3])
                return :($ad)
            end
        else
            #otherwise proceeds as normal
            return macro_parse_top(expr)
        end
    else
        # parse as normal logic program
        return macro_parse_top(expr)
    end
end

macro problog(expr)
    macro_probabilistic_parse_top(expr)
end

expr_to_prolog(ex) = macro_parse_top(ex)

expr_to_problog(ex) = macro_probabilistic_parse_top(ex)


text:
function dual_cost(rng::Random.AbstractRNG, c, v, μ::DiscreteMeasure, ν, ε; kwargs...)
    return dual_cost(rng, c, v, ν, μ, ε; kwargs...)
end
function dual_cost(
    rng::Random.AbstractRNG,
    c,
    v,
    μ,
    ν::DiscreteMeasure,
    ε;
    montecarlo_samples=10_000,
    kwargs...,
)
    # compute MC estimate of the expected c-transform with respect to `μ`
    mean_ctransform = Statistics.mean(
        ctransform(c, v, rand(rng, μ), ν, ε) for _ in 1:montecarlo_samples
    )

    return LinearAlgebra.dot(v, ν.ps) + mean_ctransform
end

function dual_v(rng::Random.AbstractRNG, c, μ::DiscreteMeasure, ν, ε; kwargs...)
    return dual_v(rng, c, ν, μ, ε; kwargs...)
end
function dual_v(
    rng::Random.AbstractRNG,
    c,
    μ,
    ν::DiscreteMeasure,
    ε;
    maxiters::Int=10_000,
    stepsize=1,
    warmup_phase=1,
    atol=0,
    rtol=iszero(atol) ? typeof(float(atol))(1//10_000) : 0,
)
    # initial iterates
    k = 1
    x = rand(rng, μ)
    τ = stepsize / (1 + √(0 / warmup_phase))
    ṽ = gradient_step(c, τ, ν, x, ε)

    # initial dual solution
    v = copy(ṽ)

    # error estimates
    Δv = similar(ṽ)
    errors = @. Δv / (atol + rtol * abs(v))

    converged = false
    while !converged && k < maxiters
        # update iterates
        k += 1
        x = rand(rng, μ)
        τ = stepsize / (1 + √((k - 1) / warmup_phase))
        gradient_step!(c, ṽ, τ, ṽ, ν, x, Δv, ε)

        # estimate error
        @. Δv = (ṽ - v) / k
        @. errors = Δv / (atol + rtol * abs(v))
        converged = maximum(abs, errors) < 1

        # update dual solution
        v .+= Δv
    end

    if k == maxiters && !converged
        @warn "Optimization algorithm did not converge. Try to adjust the parameters."
    end

    return v
end


text:
abstract type AbstractLProblem{Approx} end

struct LProblem{Approx, SCB, PP} <: AbstractLProblem{Approx}
    A::Approx
    cb::SCB
    prob::PP
    rules::ConvRules
    cf_vals::CFVals
    Tmax::Float64
    tstep::Float64
    delimiter::Int64
end

struct LProblem0{Approx, SCB, PP} <: AbstractLProblem{Approx}
    A::Approx
    cb::SCB
    prob::PP
    rules::ConvRules
    cf_vals::CFVals0
    Tmax::Float64
    tstep::Float64
    delimiter::Int64
end

function build_problem(L::Lattice, M::AbstractModel, args...; kwargs...)
    A = build_Approximation(M, L, args...)
    return build_problem(A; kwargs...)
end

function build_problem(A::abstractApproximation; Tmax::Float64 = 10.0, tstep::Float64 = 2.0^-7, delimiter::Int64 = 10, links = [(:all,:all)], axes = [:x, :y, :z])
    ARHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, links, axes)
    l = length(OSET.Observables)
    tspan = (0.0,Tmax*delimiter)
    t_means = collect(0.0:tstep:(Tmax*delimiter))
    t_cors = collect(0.0:tstep:Tmax)
    saved_values = SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    u0 = randomState(A, OSET)
    prob = ODEProblem(ARHS, u0, tspan)
    cb = SavingCallback(OSET, saved_values, u0; saveat = t_means)
    #cb = SavingCallback(OSET, saved_values; save_everystep=true)
    cf_vals = CFVals(rules, t_cors)
    return LProblem(A, cb, prob, rules, cf_vals, Tmax, tstep, delimiter)
end

function build_problem0(L::Lattice, M::AbstractModel, args...; kwargs...)
    A = build_Approximation(M, L, args...)
    return build_problem0(A; kwargs...)
end

function build_problem0(A::abstractApproximation; Tmax::Float64 = 10.0, tstep::Float64 = 2.0^-7, delimiter::Int64 = 10, links = [(:all,:all)], axes = [:x, :y, :z])
    ARHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, links, axes)
    l = length(OSET.Observables)
    tspan = (0.0,Tmax*delimiter)
    t_means = collect(0.0:tstep:(Tmax*delimiter))
    t_cors = collect(0.0:tstep:Tmax)
    saved_values = SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    u0 = randomState(A, OSET)
    prob = ODEProblem(ARHS, u0, tspan)
    cb = SavingCallback(OSET, saved_values, u0; saveat = t_means)
    #cb = SavingCallback(OSET, saved_values; save_everystep=true)
    cf_vals = CFVals0(rules, t_cors)
    return LProblem0(A, cb, prob, rules, cf_vals, Tmax, tstep, delimiter)
end

struct Logger{T1,T2} end
struct Saver{T}
    filename::String
end
Saver() = Saver{:no}("")
Saver(filename::String) = Saver{:file}(filename)
function Saver(filename::String, nTrials::Int64, cf_vals0::CFVals0)
    if !isfile(filename)
        init_file(filename, nTrials, cf_vals0)
    end
    return Saver{:fullfile}(filename)
end

function init_file(filename::AbstractString, nTrials::Int64, cf_vals0::CFVals0)
    jldopen(filename, "w") do file
        gdata = g_create(file, "data")
        file["ts"] = cf_vals0.ts
        file["str_vec"] = cf_vals0.str_vec
        file["iter"] = [0]
        for i in Base.OneTo(length(cf_vals0.str_vec))
            gdata[cf_vals0.str_vec[i]] = zeros(nTrials, size(cf_vals0.data, 1))
        end
    end
end

@inline save_data(s::Saver{:no}, cf_vals::CFVals, iter = 0) = nothing
@inline save_data(s::Saver{:file}, cf_vals::CFVals, iter) = save_data(s.filename, cf_vals, iter)
@inline save_data(s::Saver{:no}, cf_vals::CFVals0, iter) = nothing
@inline save_data(s::Saver{:fullfile}, cf_vals0::CFVals0, iter) = save_data!(s.filename, cf_vals0, iter)
function save_data(filename::AbstractString, cf_vals::CFVals, iter)
    jldopen(filename, "w") do file
        gdata = g_create(file, "data")
        gerrors = g_create(file, "errors")
        for i in Base.OneTo(length(cf_vals.str_vec))
            write(gdata, cf_vals.str_vec[i], cf_vals.data[i])
            write(gerrors, cf_vals.str_vec[i], cf_vals.errors[i])
        end
        write(file, "ts", cf_vals.ts)
        write(file, "str_vec", cf_vals.str_vec)
        write(file, "nTrials", iter)
    end
end
function save_data!(filename::AbstractString, cf_vals0::CFVals0, iter)
    jldopen(filename, "r+") do file
        gdata = file["data"]
        for i in Base.OneTo(length(cf_vals0.str_vec))
            gdata[cf_vals0.str_vec[i]][iter,:] = cf_vals0.data[i]
        end
        file["iter"][1] = iter
    end
end

function load_data(filename::AbstractString)
    jdata = load(filename)
    data = VectorOfArray([jdata["data"][key] for key in jdata["str_vec"]])
    errors = VectorOfArray([jdata["errors"][key] for key in jdata["str_vec"]])
    return CFVals(jdata["str_vec"], jdata["ts"], data, errors), jdata["nTrials"]
end
function load_data!(cf_vals::CFVals, filename::AbstractString)
    jdata = load(filename)
    for i in eachindex(jdata["str_vec"])
        key = jdata["str_vec"][i]
        cf_vals.data[:,i] = jdata["data"][key]
        cf_vals.errors[:,i] = jdata["errors"][key]
    end
    return jdata["nTrials"]
end
@inline load_data!(cf_vals::CFVals, saver::Saver{:file}) = load_data!(cf_vals, saver.filename)
@inline load_data!(cf_vals::CFVals0, saver::Saver{:fullfile}) = nothing

function load_data0(filename::AbstractString)
    file = jldopen(filename, "r")
    ts = read(file, "ts")
    str_vec = read(file, "str_vec")
    iter = read(file, "iter")[1]
    gdata = file["data"]
    datas = Dict{String, Matrix{Float64}}()
    for str in str_vec
        datas[str] = gdata[str][1:iter,:]
    end
    close(file)
    return CFArray(str_vec, ts, datas)
end

abstract type AbstractLog end
struct CMDLog <: AbstractLog end

@inline output(log_o::CMDLog, args...) = (print(args...);print("\n"))
@inline Base.close(log_o::CMDLog) = nothing

struct FileLog
    status::Int64
    offset::Float64
    filename::String
end

function FileLog(filename::AbstractString)
    if isfile(filename)
        io = open(filename, "r+")
        lc = countlines(io)
        if lc>0
            skip(io, -64)
            line = split(readline(io))
            status = parse(Int64, line[1])
            offset = status*parse(Float64, line[3])
            seekend(io)
        else
            status = 0
            offset = 0.0
        end
    else
        io = open(filename, "w")
        status = 0
        offset = 0.0
    end
    close(io)
    FileLog(status, offset, filename)
end

@inline output(log_o::FileLog, args...) = (io = open(log_o.filename, "a"); print(io, args...);print(io, "\n"); close(io))

@inline get_status(log_o::CMDLog) = 0
@inline get_status(log_o::FileLog) = log_o.status
@inline get_offset(log_o::CMDLog) = 0.0
@inline get_offset(log_o::FileLog) = log_o.offset

@inline get_string(nTrials::Int64, thread::Int64, alg::OrdinaryDiffEqAlgorithm) = @sprintf("nTr%d Alg(%s) thr%d", nTrials, alg, thread)

function get_path_string(A::abstractApproximation)
    return get_string(A.L) * "/" * get_string(A) * "/"
end
function get_path_string(A::CompositeApproximation)
    s1 = get_string(A.A1.L) * @sprintf(" + UnSp:%s", A.A2.L.cell_name)
    return s1 * "/" * get_string(A.A1) * "/"
end

function set_Logging(logger::Logger, A::abstractApproximation, path_part::AbstractString, filename, args...)
    path = set_path(logger, A, path_part)
    save_o = set_Saver(path, filename, args...)
    log_o = set_Logger(logger, path, filename)
    return log_o, save_o
end

@inline set_path(logger::Logger{:no, :cmd}, args...) = nothing
function set_path(logger::Logger, A::abstractApproximation, path_part::AbstractString)
    path = get_prefix(logger) * get_path_string(A) * path_part
    mkpath(path)
    return path
end

get_prefix(logger::Logger{:local}) = "./Data/"
get_prefix(logger::Logger{:cluster}) = ENV["WORK_DIR"] * "/Data/"

@inline set_Saver(path::Void, args...) = Saver()
set_Saver(path::AbstractString, filename::AbstractString, args...) = Saver(path * filename * ".jld", args...)
set_Saver(path::AbstractString, filename::Vector{T}, args...) where {T<:AbstractString} = [Saver(path * fname * ".jld", args...) for fname in filename]

@inline set_Logger(logger::Logger{T, :cmd}, args...) where {T} = CMDLog()
@inline set_Logger(logger::Logger{T, :cmd}, path::AbstractString, filename::AbstractString) where {T} = CMDLog()
@inline set_Logger(logger::Logger{T, :cmd}, path::AbstractString, filename::Vector{T1}) where {T, T1<:AbstractString} = CMDLog()
set_Logger(logger::Logger, path::AbstractString, filename::AbstractString) = FileLog(path * filename * ".log")
set_Logger(logger::Logger, path::AbstractString, filename::Vector{T}) where {T<:AbstractString} = FileLog(path * filename[1] * ".log")

set_Logging(LP::LProblem, nTrials::Int64, args...) = set_Logging(args...)
set_Logging(LP::LProblem0, nTrials::Int64, args...) = set_Logging(args..., nTrials, LP.cf_vals)

function simulate(LP::AbstractLProblem, nTrials::Int64, alg::OrdinaryDiffEqAlgorithm, logger::Logger{T1,T2} = Logger{:no, :cmd}(), thread::Int64 = 1, interval::Int64 = 1, RNG::AbstractRNG=Base.GLOBAL_RNG; abstol_i = 1e-14, reltol_i = 1e-6, kwargs...) where {T1,T2}
    assert(nTrials > 0)
    saved_values = LP.cb.affect!.saved_values
    rules = LP.rules
    cf_vals = LP.cf_vals
    A = LP.A
    OSET = LP.cb.affect!.save_func
    filename = (@sprintf("Tm%.2f dt%.3f del%d Obs%s/", LP.Tmax, LP.tstep, LP.delimiter, get_string(LP.rules.str_vec)), get_string(nTrials, thread, alg))
    log_o, save_o = set_Logging(LP, nTrials, logger, A, filename...)
    iter = get_status(log_o)
    offset = get_offset(log_o)
    if iter > 0
        load_data!(cf_vals, save_o)
    end
    if iter == nTrials && (T1 == :local || T1 == :cluster)
        nTrials = 2*nTrials
        #close(log_o)
        filename = (@sprintf("Tm%.2f dt%.3f del%d Obs%s/", LP.Tmax, LP.tstep, LP.delimiter, get_string(LP.rules.str_vec)), get_string(nTrials, thread, alg))
        log_o, save_o = set_Logging(LP, nTrials, logger, A, filename...)
    end
    iter += 1
    t1 = time()
    t0 = t1-offset
    integrator = init(LP.prob, alg; save_everystep = false, callback = LP.cb, abstol = abstol_i, reltol = reltol_i, kwargs...)
    reinit!(integrator, randomState(A,OSET,RNG))
    solve!(integrator)
    calculateCorrFunc!(A, saved_values, rules, cf_vals, iter)
    save_data(save_o, cf_vals, iter)
    t2 = time()
    output(log_o, @sprintf("%10d %20.5f %20.5f %10d", iter, t2-t1, (t2-t0)/iter, nTrials - iter))
    iter += 1
    while iter <= nTrials
        t1 = time()
        reinit!(integrator, randomState(A,OSET,RNG))
        solve!(integrator)
        calculateCorrFunc!(A, saved_values, rules, cf_vals, iter)
        t2 = time()
        if iter%interval == 0
            save_data(save_o, cf_vals, iter)
            output(log_o, @sprintf("%10d %20.5f %20.5f %10d", iter, t2-t1, (t2-t0)/iter, nTrials - iter))
        end
        iter += 1
    end
    save_data(save_o, cf_vals, iter)
    output(log_o, @sprintf("%10d %20.5f %20.5f %10d", nTrials, t2-t0, (t2-t0)/(iter-1), 0))
    #close(log_o)
    return cf_vals
end

export Logger

function parallel_simulate(LP::AbstractLProblem, nTrials::Int64, alg::OrdinaryDiffEqAlgorithm, logger::Logger{T1,T2} = Logger{:no, :cmd}(), interval::Int64 = 1, offset::Int64 = 0; kwargs...) where {T1,T2}
    RNGarray = initRNG(nprocs())
    @sync begin
        for i in workers()
            @spawnat i simulate(LP,nTrials, alg, logger, i+offset, interval, RNGarray[i]; kwargs...)
        end
        @spawnat 1 simulate(LP,nTrials, alg, logger, 1+offset, interval, RNGarray[1]; kwargs...)
    end
end

export parallel_simulate

function parse_dir(dirname::AbstractString = pwd())
    blacklist = ["aggregate.jld", "aggregate1.jld", "aggregate2.jld"]
    corfs = Vector{CFVals}()
    Ns = Vector{Int64}()
    for filename in readdir(dirname)
        file = joinpath(dirname, filename)
        if !(filename in blacklist) && isfile(file) && file[end-3:end] == ".jld"
            corf, Ntr = load_data(file)
            push!(corfs, corf)
            push!(Ns, Ntr)
        end
    end
    return corfs, Ns
end

function parse_dir(dirlist::Tuple{AbstractString,Vararg{AbstractString}})
    blacklist = ["aggregate.jld", "aggregate1.jld", "aggregate2.jld"]
    corfs = Vector{CFVals}()
    Ns = Vector{Int64}()
    for dirname in dirlist
        for filename in readdir(dirname)
            file = joinpath(dirname, filename)
            if !(filename in blacklist) && isfile(file) && file[end-3:end] == ".jld"
                corf, Ntr = load_data(file)
                push!(corfs, corf)
                push!(Ns, Ntr)
            end
        end
    end
    return corfs, Ns
end

function aggregate(dirname::AbstractString = pwd(), normalize::Val{T} = Val{true}()) where {T}
    corfs, Ns = parse_dir(dirname)
    return aggregate(corfs, Ns, dirname, "aggregate.jld", normalize)
end

function aggregate(dirlist::Tuple{AbstractString,Vararg{AbstractString}}, destdir::AbstractString = pwd(), normalize::Val{T} = Val{true}()) where {T}
    corfs, Ns = parse_dir(dirlist)
    mkpath(destdir)
    return aggregate(corfs, Ns, destdir, "aggregate.jld", normalize)
end

get_norm_coef(coef::Float64, ::Val{true}) = coef
get_norm_coef(coef::Float64, ::Val{false}) = 1.0

function aggregate(corfs::Vector{CFVals}, Ns::Vector{Int64}, dirname::AbstractString = pwd(), agrname = "aggregate.jld", normalize::Val{T} = Val{true}()) where {T}
    aggregate = zeros(first(corfs))
    N_total = sum(Ns)
    for i in eachindex(corfs)
        aggregate.data   .= aggregate.data  .+ corfs[i].data .* Ns[i]
        aggregate.errors .= aggregate.errors .+ corfs[i].errors
    end
    aggregate.data .= aggregate.data .* (1.0/N_total)
    for i = 1:length(corfs)
        aggregate.errors .= aggregate.errors .+ (corfs[i].data - aggregate.data).^2 .* Ns[i]
    end
    divisor = 1.0/(N_total * (N_total - 1))
    for i in 1:length(aggregate.data)
        coef = 1.0/aggregate.data[1,i]
        aggregate.data[:,i] .= aggregate.data[:,i] .* get_norm_coef(coef, normalize)
        aggregate.errors[:, i] .= sqrt.(aggregate.errors[:,i] .* divisor) .* get_norm_coef(coef, normalize)
    end
    save_data(joinpath(dirname, agrname), aggregate, N_total)
    return aggregate, N_total
end

function separate(Ns::Vector{Int64})
    i = length(Ns)
    lm = div(i,2)
    l1 = Vector{Int64}()
    l2 = Vector{Int64}()
    sum1 = 0
    sum2 = 0
    while (length(l1)<lm) && (length(l2)<lm)
        el = Ns[i]
        if sum1>sum2
            sum2 += el
            push!(l2, i)
        else
            sum1 += el
            push!(l1, i)
        end
        i -=1
    end
    if sum1 > sum2
        append!(l2, 1:i)
    else
        append!(l1, 1:i)
    end
    return l1, l2
end

function aggregate2(dirname::AbstractString = pwd(), normalize::Val{T} = Val{true}()) where {T}
    corfs, Ns = parse_dir(dirname)
    l1, l2 = separate(Ns)
    agr1, Ntot1 = aggregate(corfs[l1], Ns[l1], dirname, "aggregate1.jld", normalize)
    agr2, Ntot2 = aggregate(corfs[l2], Ns[l2], dirname, "aggregate2.jld", normalize)
    return agr1, Ntot1, agr2, Ntot2
end

function aggregate2(dirlist::Tuple{AbstractString,Vararg{AbstractString}}, destdir::AbstractString = pwd(), normalize::Val{T} = Val{true}()) where {T}
    corfs, Ns = parse_dir(dirlist)
    l1, l2 = separate(Ns)
    mkpath(destdir)
    agr1, Ntot1 = aggregate(corfs[l1], Ns[l1], destdir, "aggregate1.jld", normalize)
    agr2, Ntot2 = aggregate(corfs[l2], Ns[l2], destdir, "aggregate2.jld", normalize)
    return agr1, Ntot1, agr2, Ntot2
end


text:
using MatterPower
using Test

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


text:
using MPIMeasurements

# define Positions
positions = loadTDesign(8,36,42Unitful.mm,[16.0,0.0,0.0]Unitful.mm)

# create Scanner
scanner = MPIScanner("BrukerScanner.toml")
#scanner = MPIScanner("DummyScanner.toml")
robot = getRobot(scanner)
safety = getSafety(scanner)
gauss = getGaussMeter(scanner)
gausMeterRange = 1

mfMeasObj = MagneticFieldMeas(gauss, positions, gausMeterRange)

# perform measurement
res = performTour!(robot, safety, positions, mfMeasObj)

movePark(robot)

#filenameField = joinpath(homedir(),"measurementtmp","MagneticField_G2.4_FF0mm0mm0mm.h5")
#saveMagneticFieldAsHDF5(mfMeasObj, filenameField)

filenameField = joinpath(homedir(),"measurementtmp","MagneticField_G0.6_FF10mm10mm10mm.h5")
saveMagneticFieldAsHDF5(mfMeasObj, filenameField)


text:
using Revise
using Plots
using Basins
using DynamicalSystems

function newton_map(dz,z, p, n)
    f(x) = x^p[1]-1
    df(x)= p[1]*x^(p[1]-1)
    z1 = z[1] + im*z[2]
    dz1 = f(z1)/df(z1)
    z1 = z1 - dz1
    dz[1]=real(z1)
    dz[2]=imag(z1)
    return
end

# dummy function to keep the initializator happy
function newton_map_J(J,z0, p, n)
   return
end
ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3] , newton_map_J)
integ  = integrator(ds)

a=zeros(1,8)
b=zeros(1,8)
c=zeros(1,8)
d=zeros(1,8)

# for (k,res) in enumerate(range(100,800,step=100))
#
#     xg=range(-1.,1.,length=res)
#     yg=range(-1.,1.,length=res)
#
#     @time bsn=Basins.basins_map2D(xg, yg, integ)
#
#     a[k],De,e,f=uncertainty_exponent(bsn,integ)
#     b[k],e,f=Basins.static_estimate(xg,yg,bsn)
#     c[k],e,f=Basins.ball_estimate(xg,yg,bsn)
#     d[k]=2-box_counting_dim(xg,yg,bsn)
# end

xg=range(-1.,1.,length=300)
yg=range(-1.,1.,length=300)

@time bsn=Basins.basins_map2D(xg, yg, integ)

b=zeros(1,40)
c=zeros(1,40)

precision = 1e-5
for k in 1:40
        #a[k],De,e,f=uncertainty_exponent(bsn,integ)
        b[k],e,f=Basins.static_estimate(xg,yg,bsn; precision=precision)
        c[k],e,f=Basins.ball_estimate(xg,yg,bsn; precision=precision)
        #d[k]=2-box_counting_dim(xg,yg,bsn)
end

using Statistics
println(mean(b), "  ", mean(c))
#plot(xg,yg,bsn.basin',seriestype=:heatmap)
plot(b')
plot!(c')

#sa,sb = compute_saddle(integ, bsn, [1], [2,3]; N=100)

#s=Dataset(sa)
#plot!(s[:,1],s[:,2],seriestype=:scatter)


text:
# This file was generated, do not modify it. # hide
#hideall
println("From code:  **hello**")

max_stars_repo_name: abhishalya/julialang-judoc

max_issues_repo_name: abhishalya/julialang-judoc

max_forks_repo_name: abhishalya/julialang-judoc


text:
module Diagram

abstract type halfEdge end # this hack is currently necessary for definition of circular types in Julia, see https://github.com/JuliaLang/julia/issues/269

mutable struct Region
   generator::Tuple{Real, Real}
   borderHead::Union{halfEdge, Nothing} # first half edge of the border; we can traverse the border with he = he.next;
end

mutable struct HalfEdge <: halfEdge
   origin::Union{Tuple{Real, Real}, Nothing}

   isFixed::Bool

   twin::Union{HalfEdge, Nothing}
   next::Union{HalfEdge, Nothing}
   prev::Union{HalfEdge, Nothing}

   generator::Tuple{Real, Real}
end

function makeTwins(a::HalfEdge, b::HalfEdge)
   a.twin = b
   b.twin = a
end

function concat(a::HalfEdge, b::HalfEdge)
   a.next = b
   b.prev = a
end

mutable struct DCEL # Doubly-connected edge list
   regions::Array{Region}
   halfEdges::Array{HalfEdge}
end

function DCEL(points::Array{Tuple{Real, Real}, 1})
	n = size(points)[1]

   regions = Array{Region}(undef, 1, n)
   for i in 1:n
      regions[i] = Region(points[i], nothing)
   end

   halfEdges = HalfEdge[]

	return DCEL(regions, halfEdges)
end

function borderCoordinates(f::Region)
   x = []
   y = []

   if f.borderHead == nothing
      return x, y
   end

   he = f.borderHead

   if he.origin != nothing
      append!(x, he.origin[1])
      append!(y, he.origin[2])
   end
   he = he.next

   i = 0
   while he != nothing && he != f.borderHead
      if he.origin != nothing
         append!(x, he.origin[1])
         append!(y, he.origin[2])
      end
      he = he.next
      i += 1
   end

   return x, y
end

function regionBorders(V::DCEL)
   borders = []

   for region in V.regions
      append!(borders, [borderCoordinates(region)])
   end

   return borders
end

function centroidAndArea(region::Region)
   A = Cx = Cy = 0

   he = region.borderHead
   while true
      shoelace = he.origin[1]*he.next.origin[2] - he.origin[2]*he.next.origin[1]
      A += shoelace
      Cx += (he.origin[1] + he.next.origin[1]) * shoelace
      Cy += (he.origin[2] + he.next.origin[2]) * shoelace

      he = he.next
      if he == region.borderHead
         break
      end
   end

   A /= 2
   Cx /= 6*A
   Cy /= 6*A

   return ((Cx, Cy), A)
end

function centroidsAndAreas(V::DCEL)
   centroids = Array{Tuple{Real, Real}, 1}([])
   areas = Array{Real, 1}([])
   for region in V.regions
      res = centroidAndArea(region)
      push!(centroids, (res[1]))
      push!(areas, (res[2]))
   end

   return centroids, areas
end

end # module


text:
"""
    resample2(v, rate)

Performs time-domain resampling of `v` by the ratio `rate` using cubic interpolation.
When downsampling the signal is lowpassed before resampling. When up-sampling it is
lowpassed afterwards to remove interpolation noise.

Similarities to `DSP.resample`

- The signals are aligned such that the first sample of `v` corresponds to the first
  sample of `resample2(v, rate)`.

Differences from `DSP.resample`:

- this guarantees that `length(resample2(v, rate)) == round(Int, length(v)*rate)`
- `resample2` uses :reflect_at boundary conditions to avoid ringing when the signal
  starts far from 0.
"""
function resample2(v, rate, lpf=resample_filter2(rate))
    L = length(v)
    Lf = length(lpf)
    @assert isodd(Lf)
    fcomp = Lf÷2
    # reflecting adds Lf-1 samples
    # convolution adds Lf-1 samples
    if rate <= 1
        vf = conv(lpf, pad(v, :reflect_at, fcomp))
        map(0:round(Int, L*rate)-1) do i
            cubic(vf, i/rate+Lf, :reflect_at)
        end
    else
        out = map(0:round(Int, L*rate)-1) do i
            cubic(v, i/rate+1, :reflect_at)
        end
        conv(lpf, pad(out, :reflect_at, fcomp))[Lf:end-Lf+1]
    end
end

function pad(v, bc, pre, post=pre)
    # TODO: support padding longer than the signal, which is useful for large resampling
    # factors
    prepad, postpad = if bc == :reflect_at
        reverse(v[2:pre+1]), reverse(v[end-post:end-1])
    else
        throw(ArgumentError("Unsupported boundary condition type"))
    end

    [prepad; v; postpad]
end

function resample_filter2(rate)
    # DSP.jl doesn't force FIRWindow to be odd-length, so we do that here
    # see https://github.com/JuliaDSP/DSP.jl/issues/345
    cutoff = rate <= 1 ? rate : 1/rate
    ord, α = kaiserord(0.1*cutoff)
    iseven(ord) && (ord += 1)
    digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(ord, α)))
end

# generate data without information too close to nyquist
# x = resample2(rand(100).+10, 2)
# r = 0.5
# plot(0:length(x)-1, x)
# x2 = resample2(x, r)
# plot!((0:length(x2)-1)/r, x2)
#
# x = randn(100000)
# @benchmark resample(x, 0.05)
# @benchmark resample2(x, 0.05)
# currently `resample2` is about 4-5x faster than `resample` for this case

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
module SIMDTest

using Base.Test
using ForwardDiff: Dual, valtype

const DUALS = (Dual(1., 2., 3., 4.),
               Dual(1., 2., 3., 4., 5.),
               Dual(Dual(1., 2.), Dual(3., 4.)))

function simd_sum(x::Vector{T}) where T
    s = zero(T)
    @simd for i in eachindex(x)
        @inbounds s = s + x[i]
    end
    return s
end

for D in map(typeof, DUALS)
    plus_bitcode = sprint(io -> code_llvm(io, +, (D, D)))
    @test contains(plus_bitcode, "fadd <4 x double>")

    minus_bitcode = sprint(io -> code_llvm(io, -, (D, D)))
    @test contains(minus_bitcode, "fsub <4 x double>")

    times_bitcode = sprint(io -> code_llvm(io, *, (D, D)))
    @test ismatch(r"fadd \<.*?x double\>", times_bitcode)
    @test ismatch(r"fmul \<.*?x double\>", times_bitcode)

    div_bitcode = sprint(io -> code_llvm(io, /, (D, D)))
    @test ismatch(r"fadd \<.*?x double\>", div_bitcode)
    @test ismatch(r"fmul \<.*?x double\>", div_bitcode)

    exp_bitcode = sprint(io -> code_llvm(io, ^, (D, D)))
    @test ismatch(r"fadd \<.*?x double\>", exp_bitcode)
    if !(valtype(D) <: Dual)
        # see https://github.com/JuliaDiff/ForwardDiff.jl/issues/167
        @test ismatch(r"fmul \<.*?x double\>", exp_bitcode)

        # see https://github.com/JuliaDiff/ForwardDiff.jl/pull/201
        sum_bitcode = sprint(io -> code_llvm(io, simd_sum, (Vector{D},)))
        @test ismatch(r"fadd \<.*?x double\>", sum_bitcode)
    end
end

end # module