text:
module Html

import Markdown, Logging, Gumbo, Reexport, OrderedCollections, Millboard, HTTP, YAML

Reexport.@reexport using Genie

import Genie.Renderer
import Genie.Renderer: vars
Reexport.@reexport using HttpCommon

const DEFAULT_LAYOUT_FILE = :app
const LAYOUTS_FOLDER = "layouts"

const HTML_FILE_EXT = ".jl"
const TEMPLATE_EXT  = ".jl.html"
const MARKDOWN_FILE_EXT = [".md", ".jl.md"]

const SUPPORTED_HTML_OUTPUT_FILE_FORMATS = [TEMPLATE_EXT]

const HTMLString = String
const HTMLParser = Gumbo

const NBSP_REPLACEMENT = (" "=>"!!nbsp;;")

const NORMAL_ELEMENTS = [ :html, :head, :body, :title, :style, :address, :article, :aside, :footer,
                          :header, :h1, :h2, :h3, :h4, :h5, :h6, :hgroup, :nav, :section,
                          :dd, :div, :d, :dl, :dt, :figcaption, :figure, :li, :main, :ol, :p, :pre, :ul, :span,
                          :a, :abbr, :b, :bdi, :bdo, :cite, :code, :data, :dfn, :em, :i, :kbd, :mark,
                          :q, :rp, :rt, :rtc, :ruby, :s, :samp, :small, :spam, :strong, :sub, :sup, :time,
                          :u, :var, :wrb, :audio, :void, :embed, :object, :canvas, :noscript, :script,
                          :del, :ins, :caption, :col, :colgroup, :table, :tbody, :td, :tfoot, :th, :thead, :tr,
                          :button, :datalist, :fieldset, :label, :legend, :meter,
                          :output, :progress, :select, :option, :textarea, :details, :dialog, :menu, :menuitem, :summary,
                          :slot, :template, :blockquote, :center, :iframe] #TODO: Gumbo has a problem and strips away <template>
const VOID_ELEMENTS   = [:base, :link, :meta, :hr, :br, :area, :img, :track, :param, :source, :input]
const CUSTOM_ELEMENTS = [:form, :select]
const NON_EXPORTED = [:main, :map, :filter]
const SVG_ELEMENTS = [:animate, :circle, :animateMotion, :animateTransform, :clipPath, :defs, :desc, :discard,
                      :ellipse, :feComponentTransfer, :feComposite, :feDiffuseLighting, :feBlend, :feColorMatrix,
                      :feConvolveMatrix, :feDisplacementMap, :feDistantLight, :feDropShadow, :feFlood, :feFuncA, :feFuncB, :feFuncG,
                      :feFuncR, :feGaussianBlur, :feImage, :feMerge, :feMergeNode, :feMorphology, :feOffset, :fePointLight, :feSpecularLighting,
                      :feSpotLight, :feTile, :feTurbulence, :foreignObject, :g, :hatch, :hatchpath, :image, :line, :linearGradient,
                      :marker, :mask, :metadata, :mpath, :path, :pattern, :polygon, :polyline, :radialGradient, :rect, :set, :stop, :svg,
                      :switch, :symbol, :text, :textPath, :tspan, :use, :view]

const EMBEDDED_JULIA_PLACEHOLDER = "~~~~~|~~~~~"

export HTMLString, html, doc, doctype
export @yield, collection, view!, for_each
export partial, template

task_local_storage(:__yield, "")

include("MdHtml.jl")

"""
    normal_element(f::Function, elem::String, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString

Generates a HTML element in the form <...></...>
"""
function normal_element(f::Function, elem::Any, args::Vector = [], attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  normal_element(f(), string(elem), args, attrs...)
end
function normal_element(children::Union{String,Vector{String}}, elem::Any, args::Vector, attrs::Pair{Symbol,Any}) :: HTMLString
  normal_element(children, string(elem), args, Pair{Symbol,Any}[attrs])
end
function normal_element(children::Tuple, elem::Any, args::Vector, attrs::Pair{Symbol,Any}) :: HTMLString
  normal_element([children...], string(elem), args, Pair{Symbol,Any}[attrs])
end
function normal_element(children::Union{String,Vector{String}}, elem::Any, args::Vector, attrs...) :: HTMLString
  normal_element(children, string(elem), args, Pair{Symbol,Any}[attrs...])
end
function normal_element(children::Union{String,Vector{String}}, elem::Any, args::Vector = [], attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  children = join(children)

  idx = findfirst(x -> x == :inner, getindex.(attrs, 1))
  if !isnothing(idx)
    children *= join(attrs[idx][2])
    deleteat!(attrs, idx)
  end

  ii = (x -> x isa Symbol && occursin(Genie.config.html_parser_char_dash, "$x")).(args)
  args[ii] .= Symbol.(replace.(string.(args[ii]), Ref(Genie.config.html_parser_char_dash=>"-")))

  elem = normalize_element(string(elem))
  attribs = rstrip(attributes(attrs))
  string("<", elem, (isempty(attribs) ? "" : " $attribs"), (isempty(args) ? "" : " $(join(args, " "))"), ">", prepare_template(children), "</", elem, ">")
end
function normal_element(elem::Any, attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  normal_element("", string(elem), attrs...)
end
function normal_element(elems::Vector, elem::Any, args = [], attrs...) :: HTMLString
  io = IOBuffer()

  for e in elems
    e === nothing && continue

    if isa(e, Vector)
      print(io, join(e))
    elseif isa(e, Function)
      print(io, e(), "\n")
    else
      print(io, e, "\n")
    end
  end

  normal_element(String(take!(io)), string(elem), args, attrs...)
end
function normal_element(_::Nothing, __::Any) :: HTMLString
  ""
end
function normal_element(children::Vector{T} where T, elem::Any, args::Vector{T} where T) :: HTMLString
  normal_element(join([(isa(f, Function) ? f() : string(f)) for f in children]), elem, args)
end

"""
    prepare_template(s::String)
    prepare_template{T}(v::Vector{T})

Cleans up the template before rendering (ex by removing empty nodes).
"""
function prepare_template(s::String) :: String
  s
end
function prepare_template(v::Vector{T})::String where {T}
  filter!(v) do (x)
    ! isa(x, Nothing)
  end

  join(v)
end

"""
    attributes(attrs::Vector{Pair{Symbol,String}} = Vector{Pair{Symbol,String}}()) :: Vector{String}

Parses HTML attributes.
"""
function attributes(attrs::Vector{Pair{Symbol,Any}} = Vector{Pair{Symbol,Any}}()) :: String
  a = IOBuffer()

  for (k,v) in attrs
    v === nothing && continue # skip nothing values
    v == "false" && (v = false)

    print(a, attrparser(k, v))
  end

  String(take!(a))
end

function attrparser(k::Symbol, v::Bool) :: String
  v ? "$(k |> parseattr) " : ""
end

function attrparser(k::Symbol, v::Union{AbstractString,Symbol}) :: String
  isempty(string(v)) && return attrparser(k, true)
  "$(k |> parseattr)=\"$(v)\" "
end

function attrparser(k::Symbol, v::Any) :: String
  default_attrparser(k, v)
end

function default_attrparser(k::Symbol, v::Any) :: String
  "$(k |> parseattr)=$(v) "
end

"""
    parseattr(attr) :: String

Converts Julia keyword arguments to HTML attributes with illegal Julia chars.
"""
function parseattr(attr) :: String
  attr = string(attr)
  endswith(attr, Genie.config.html_parser_char_at) && (attr = string("@", attr[1:end-(length(Genie.config.html_parser_char_at))]))
  endswith(attr, Genie.config.html_parser_char_column) && (attr = string(":", attr[1:end-(length(Genie.config.html_parser_char_column))]))

  attr = replace(attr, Regex(Genie.config.html_parser_char_dash)=>"-")
  attr = replace(attr, Regex(Genie.config.html_parser_char_column * Genie.config.html_parser_char_column)=>":")
  replace(attr, Regex(Genie.config.html_parser_char_dot)=>".")
end

"""
    normalize_element(elem::String)

Cleans up problematic characters or DOM elements.
"""
function normalize_element(elem::String) :: String
  replace(elem, Genie.config.html_parser_char_dash=>"-")
end

"""
    denormalize_element(elem::String)

Replaces `-` with the char defined to replace dashes, as Julia does not support them in names.
"""
function denormalize_element(elem::String) :: String
  uppercase(elem) == elem && (elem = lowercase(elem))
  elem = replace(elem, "-"=>Genie.config.html_parser_char_dash)
  endswith(elem, "_") ? elem[1:end-1] : elem
end

"""
    void_element(elem::String, attrs::Vector{Pair{Symbol,String}} = Vector{Pair{Symbol,String}}()) :: HTMLString

Generates a void HTML element in the form <...>
"""
function void_element(elem::String, args = [], attrs::Vector{Pair{Symbol,Any}} = Pair{Symbol,Any}[]) :: HTMLString
  attribs = rstrip(attributes(attrs))
  string("<", normalize_element(elem), (isempty(attribs) ? "" : " $attribs"), (isempty(args) ? "" : " $(join(args, " "))"), Genie.config.html_parser_close_tag, ">")
end

"""
    get_template(path::String; partial::Bool = true, context::Module = @__MODULE__, vars...) :: Function

Resolves the inclusion and rendering of a template file
"""
function get_template(path::String; partial::Bool = true, context::Module = @__MODULE__, vars...) :: Function
  path, extension = Genie.Renderer.view_file_info(path, SUPPORTED_HTML_OUTPUT_FILE_FORMATS)

  f_name = Genie.Renderer.function_name(string(path, partial)) |> Symbol
  mod_name = Genie.Renderer.m_name(string(path, partial)) * ".jl"
  f_path = joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name)
  f_stale = Genie.Renderer.build_is_stale(path, f_path)

  if f_stale || ! isdefined(context, f_name)
    if extension in MARKDOWN_FILE_EXT
      path = MdHtml.md_to_html(path, context = context)
    end

    f_stale && Genie.Renderer.build_module(html_to_julia(path, partial = partial, extension = extension), path, mod_name)

    return Base.include(context, joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name))
  end

  getfield(context, f_name)
end

"""
Outputs document's doctype.
"""
function doctype(doctype::Symbol = :html) :: HTMLString
  "<!DOCTYPE $doctype>"
end

"""
Outputs document's doctype.
"""
function doc(html::String) :: HTMLString
  doctype() * "\n" * html
end
function doc(doctype::Symbol, html::String) :: HTMLString
  doctype(doctype) * "\n" * html
end

"""
    parseview(data::String; partial = false, context::Module = @__MODULE__) :: Function

Parses a view file, returning a rendering function. If necessary, the function is JIT-compiled, persisted and loaded into memory.
"""
function parseview(data::String; partial = false, context::Module = @__MODULE__) :: Function
  data_hash = hash(data)
  path = "Genie_" * string(data_hash)

  func_name = Genie.Renderer.function_name(string(data_hash, partial)) |> Symbol
  mod_name = Genie.Renderer.m_name(string(path, partial)) * ".jl"
  f_path = joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name)
  f_stale = Genie.Renderer.build_is_stale(f_path, f_path)

  if f_stale || ! isdefined(context, func_name)
    if f_stale
      Genie.Renderer.build_module(string_to_julia(data, partial = partial), path, mod_name)
    end

    return Base.include(context, joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name))
  end

  getfield(context, func_name)
end

"""
    render(data::String; context::Module = @__MODULE__, layout::Union{String,Nothing} = nothing, vars...) :: Function

Renders the string as an HTML view.
"""
function render(data::String; context::Module = @__MODULE__, layout::Union{String,Nothing} = nothing, vars...) :: Function
  Genie.Renderer.registervars(; context = context, vars...)

  if layout !== nothing
    task_local_storage(:__yield, parseview(data, partial = true, context = context))
    parseview(layout, partial = false, context = context)
  else
    parseview(data, partial = false, context = context)
  end
end

"""
    render(viewfile::Genie.Renderer.FilePath; layout::Union{Nothing,Genie.Renderer.FilePath} = nothing, context::Module = @__MODULE__, vars...) :: Function

Renders the template file as an HTML view.
"""
function render(viewfile::Genie.Renderer.FilePath; layout::Union{Nothing,Genie.Renderer.FilePath} = nothing, context::Module = @__MODULE__, vars...) :: Function
  Genie.Renderer.registervars(; context = context, vars...)

  if layout !== nothing
    task_local_storage(:__yield, get_template(string(viewfile); partial = true, context = context, vars...))
    get_template(string(layout); partial = false, context = context, vars...)
  else
    get_template(string(viewfile); partial = false, context = context, vars...)
  end
end

"""
    parsehtml(input::String; partial::Bool = true) :: String

"""
function parsehtml(input::String; partial::Bool = true) :: String
  parsehtml(
    HTMLParser.parsehtml(replace(input, NBSP_REPLACEMENT), preserve_whitespace = false).root,
    0, partial = partial)
end

function Genie.Renderer.render(::Type{MIME"text/html"}, data::String; context::Module = @__MODULE__, layout::Union{String,Nothing} = nothing, vars...) :: Genie.Renderer.WebRenderable
  try
    render(data; context = context, layout = layout, vars...) |> Genie.Renderer.WebRenderable
  catch ex
    isa(ex, KeyError) && Genie.Renderer.changebuilds() # it's a view error so don't reuse them
    rethrow(ex)
  end
end

function Genie.Renderer.render(::Type{MIME"text/html"}, viewfile::Genie.Renderer.FilePath; layout::Union{Nothing,Genie.Renderer.FilePath} = nothing, context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  try
    render(viewfile; layout = layout, context = context, vars...) |> Genie.Renderer.WebRenderable
  catch ex
    isa(ex, KeyError) && Genie.Renderer.changebuilds() # it's a view error so don't reuse them
    rethrow(ex)
  end
end

function html(resource::Genie.Renderer.ResourcePath, action::Genie.Renderer.ResourcePath; layout::Genie.Renderer.ResourcePath = DEFAULT_LAYOUT_FILE,
                context::Module = @__MODULE__, status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), vars...) :: Genie.Renderer.HTTP.Response
  html(Genie.Renderer.Path(joinpath(Genie.config.path_resources, string(resource), Renderer.VIEWS_FOLDER, string(action)));
        layout = Genie.Renderer.Path(joinpath(Genie.config.path_app, LAYOUTS_FOLDER, string(layout))),
        context = context, status = status, headers = headers, vars...)
end

"""
    html(data::String; context::Module = @__MODULE__, status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), layout::Union{String,Nothing} = nothing, vars...) :: HTTP.Response

Parses the `data` input as HTML, returning a HTML HTTP Response.

# Arguments
- `data::String`: the HTML string to be rendered
- `context::Module`: the module in which the variables are evaluated (in order to provide the scope for vars). Usually the controller.
- `status::Int`: status code of the response
- `headers::HTTPHeaders`: HTTP response headers
- `layout::Union{String,Nothing}`: layout file for rendering `data`

# Example
```jldoctest
julia> html("<h1>Welcome \$(vars(:name))</h1>", layout = "<div><% @yield %></div>", name = "Adrian")
HTTP.Messages.Response:
"
HTTP/1.1 200 OK
Content-Type: text/html; charset=utf-8

<html><head></head><body><div><h1>Welcome Adrian</h1>
</div></body></html>"
```
"""
function html(data::String; context::Module = @__MODULE__, status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(),
                            layout::Union{String,Nothing,Genie.Renderer.FilePath} = nothing, forceparse::Bool = false, vars...) :: Genie.Renderer.HTTP.Response
  isa(layout, Genie.Renderer.FilePath) && (layout = read(layout, String))

  if occursin(raw"$", data) || occursin("<%", data) || layout !== nothing || forceparse
    Genie.Renderer.WebRenderable(Genie.Renderer.render(MIME"text/html", data; context = context, layout = layout, vars...), status, headers) |> Genie.Renderer.respond
  else
    Genie.Renderer.WebRenderable(body = data, status = status, headers = headers) |> Genie.Renderer.respond
  end
end

"""
    html(md::Markdown.MD; context::Module = @__MODULE__, status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), layout::Union{String,Nothing} = nothing, forceparse::Bool = false, vars...) :: Genie.Renderer.HTTP.Response

Markdown view rendering
"""
function html(md::Markdown.MD; context::Module = @__MODULE__, status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), layout::Union{String,Nothing} = nothing, forceparse::Bool = false, vars...) :: Genie.Renderer.HTTP.Response
  data = MdHtml.eval_markdown(string(md)) |> Markdown.parse |> Markdown.html
  for kv in vars
    data = replace(data, ":" * string(kv[1]) => "\$" * string(kv[1]))
  end

  html(data; context = context, status = status, headers = headers, layout = layout, forceparse = forceparse, vars...)
end

"""
    html(viewfile::FilePath; layout::Union{Nothing,FilePath} = nothing,
          context::Module = @__MODULE__, status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), vars...) :: HTTP.Response

Parses and renders the HTML `viewfile`, optionally rendering it within the `layout` file. Valid file format is `.html.jl`.

# Arguments
- `viewfile::FilePath`: filesystem path to the view file as a `Renderer.FilePath`, ie `Renderer.FilePath("/path/to/file.html.jl")`
- `layout::FilePath`: filesystem path to the layout file as a `Renderer.FilePath`, ie `Renderer.FilePath("/path/to/file.html.jl")`
- `context::Module`: the module in which the variables are evaluated (in order to provide the scope for vars). Usually the controller.
- `status::Int`: status code of the response
- `headers::HTTPHeaders`: HTTP response headers
"""
function html(viewfile::Genie.Renderer.FilePath; layout::Union{Nothing,Genie.Renderer.FilePath} = nothing,
                context::Module = @__MODULE__, status::Int = 200, headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders(), vars...) :: Genie.Renderer.HTTP.Response
  Genie.Renderer.WebRenderable(Genie.Renderer.render(MIME"text/html", viewfile; layout = layout, context = context, vars...), status, headers) |> Genie.Renderer.respond
end

"""
    safe_attr(attr) :: String

Replaces illegal Julia characters from HTML attributes with safe ones, to be used as keyword arguments.
"""
function safe_attr(attr) :: String
  attr = string(attr)

  occursin("-", attr) && replace!(attr, "-"=>Genie.config.html_parser_char_dash)
  occursin(":", attr) && replace!(attr, ":"=>Genie.config.html_parser_char_column)
  occursin("@", attr) && replace!(attr, "@"=>Genie.config.html_parser_char_at)
  occursin(".", attr) && replace!(attr, "."=>Genie.config.html_parser_char_dot)

  attr
end

"""
    parsehtml(elem, output, depth; partial = true) :: String

Parses a HTML tree structure into a `string` of Julia code.
"""
function parsehtml(elem::HTMLParser.HTMLElement, depth::Int = 0; partial::Bool = true) :: String
  io = IOBuffer()

  tag_name = denormalize_element(string(HTMLParser.tag(elem)))

  invalid_tag = partial && (tag_name == "html" || tag_name == "head" || tag_name == "body")

  if tag_name == "script" && in("type", collect(keys(HTMLParser.attrs(elem)))) && HTMLParser.attrs(elem)["type"] == "julia/eval"
    isempty(HTMLParser.children(elem)) || print(io, repeat("\t", depth), string(HTMLParser.children(elem)[1].text), "\n")
  else
    mdl = isdefined(@__MODULE__, Symbol(tag_name)) ? string(@__MODULE__, ".") : ""
    print(io, repeat("\t", depth), ( invalid_tag ? "" : "$mdl$(tag_name)(" ) )

    attributes = IOBuffer()
    attributes_keys = String[]
    attributes_values = String[]

    for (k,v) in HTMLParser.attrs(elem)
      # x = v
      k = string(k) |> lowercase
      if occursin("\"", v)
        # we need to escape double quotes but only if it's not embedded Julia code
        mx = collect(eachmatch(r"\$\((.*?)\)|<%(.*?)%>", string(v)))
        if ! isempty(mx)
          index = 1
          for value in mx
            value === nothing && continue

            label = "$(EMBEDDED_JULIA_PLACEHOLDER)_$(index)"
            v = replace(string(v), value.match => label)
            index += 1
          end

          if occursin("\"", v) # do we still need to escape double quotes?
            v = replace(string(v), "\"" => "\\\"") # we need to escape " as this will be inside julia strings
          end

          # now we need to put back the embedded Julia
          index = 1
          for value in mx
            value === nothing && continue

            label = "$(EMBEDDED_JULIA_PLACEHOLDER)_$(index)"
            v = replace(string(v), label => value.match)
            index += 1
          end
        elseif occursin("\"", v) # do we still need to escape double quotes?
          v = replace(string(v), "\"" => "'") #
        end
      end

      if startswith(k, raw"$") # do not process embedded julia code
        print(attributes, k[2:end], ", ") # strip the $, this is rendered directly in Julia code
        continue
      end

      if occursin("-", k) ||
          occursin(":", k) ||
          occursin("@", k) ||
          occursin(".", k) ||
          occursin("for", k)

        push!(attributes_keys, Symbol(k) |> repr)

        v = string(v) |> repr
        occursin(raw"\$", v) && (v = replace(v, raw"\$"=>raw"$"))
        push!(attributes_values, v)
      else
        print(attributes, """$k="$v" """, ", ")
      end
    end

    attributes_string = String(take!(attributes))
    endswith(attributes_string, ", ") && (attributes_string = attributes_string[1:end-2])

    print(io, attributes_string)
    ! isempty(attributes_string) && ! isempty(attributes_keys) && print(io, ", ")
    ! isempty(attributes_keys) &&
      print(io, "; NamedTuple{($(join(attributes_keys, ", "))$(length(attributes_keys) == 1 ? ", " : ""))}(($(join(attributes_values, ", "))$(length(attributes_keys) == 1 ? ", " : "")))...")

    invalid_tag || print(io, ")")

    inner = ""
    if ! isempty(HTMLParser.children(elem))
      children_count = size(HTMLParser.children(elem))[1]

      invalid_tag ? print(io, " \n") : print(io, " do;[\n")

      idx = 0
      for child in HTMLParser.children(elem)
        idx += 1
        inner *= isa(child, HTMLParser.HTMLText) ? parsehtml(child, depth + 1) : parsehtml(child, depth + 1, partial = partial)
        if idx < children_count
          if  ( isa(child, HTMLParser.HTMLText) ) ||
              ( isa(child, HTMLParser.HTMLElement) &&
              ( ! in("type", collect(keys(HTMLParser.attrs(child)))) ||
                ( in("type", collect(keys(HTMLParser.attrs(child)))) && (HTMLParser.attrs(child)["type"] != "julia/eval") ) ) )
              isempty(inner) || (inner = string(repeat("\t", depth), inner, "\n"))
          end
        end
      end

      if ! isempty(inner)
        endswith(inner, "\n\n") && (inner = inner[1:end-2])
        print(io, inner, repeat("\t", depth))
      end

      invalid_tag ? print(io, "\n") : print(io, ";]end\n")
    end

  end

  String(take!(io))
end

function parsehtml(elem::HTMLParser.HTMLText, depth::Int = 0; partial::Bool = true) :: String
  content = elem.text
  endswith(content, "\"") && (content *= Char(0x0))
  content = replace(content, NBSP_REPLACEMENT[2]=>NBSP_REPLACEMENT[1])
  string(repeat("\t", depth), "\"\"\"$(content)\"\"\"")
end

"""
    html_to_julia(file_path::String; partial = true) :: String

Converts a HTML document to Julia code.
"""
function html_to_julia(file_path::String; partial = true, extension = TEMPLATE_EXT) :: String
  to_julia(file_path, parse_template; partial = partial, extension = extension)
end

"""
    string_to_julia(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts string view data to Julia code
"""
function string_to_julia(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend::String = "\n") :: String
  to_julia(content, parse_string, partial = partial, f_name = f_name, prepend = prepend)
end

"""
    to_julia(input::String, f::Function; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts an input file to Julia code
"""
  function to_julia(input::String, f::Union{Function,Nothing};
                  partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend::String = "\n", extension = TEMPLATE_EXT) :: String
  f_name = (f_name === nothing) ? Genie.Renderer.function_name(string(input, partial)) : f_name

  string("function $(f_name)(; $(Genie.Renderer.injectkwvars())) \n",
          "
          [
          ",

          prepend,

          (partial ? "" : "\nGenie.Renderer.Html.doctype() \n"),

          f !== nothing ? f(input; partial = partial, extension = extension) : input,

          "
          ]
          end
          ")
end

"""
    partial(path::String; context::Module = @__MODULE__, vars...) :: String

Renders (includes) a view partial within a larger view or layout file.
"""
function partial(path::String; context::Module = @__MODULE__, kwvars...) :: String
  for (k,v) in kwvars
    vars()[k] = v
  end

  template(path, partial = true, context = context)
end

"""
    template(path::String; partial::Bool = true, context::Module = @__MODULE__, vars...) :: String

Renders a template file.
"""
function template(path::String; partial::Bool = true, context::Module = @__MODULE__, vars...) :: String
  try
    get_template(path; partial = partial, context = context, vars...)() |> join
  catch ex
    if isa(ex, MethodError) && (string(ex.f) == "get_template" || startswith(string(ex.f), "func_"))
      Base.invokelatest(get_template(path; partial = partial, context = context, vars...)) |> join
    else
      rethrow(ex)
    end
  end
end

"""
    read_template_file(file_path::String) :: String

Reads `file_path` template from disk.
"""
function read_template_file(file_path::String; extension = TEMPLATE_EXT) :: String
  io = IOBuffer()

  extension == HTML_FILE_EXT && print(io, """
  \"\"\"
  """)

  open(file_path) do f
    for line in enumerate(eachline(f))
      print(io, parsetags(line), "\n")
    end
  end

  extension == HTML_FILE_EXT && print(io, """
  \"\"\"
  """)

  String(take!(io))
end

"""
    parse_template(file_path::String; partial = true) :: String

Parses a HTML file into Julia code.
"""
function parse_template(file_path::String; partial::Bool = true, extension = TEMPLATE_EXT) :: String
  parse(read_template_file(file_path; extension = extension)::String; partial = partial)
end

"""
    parse_string(data::String; partial = true) :: String

Parses a HTML string into Julia code.
"""
function parse_string(data::String; partial::Bool = true, extension = TEMPLATE_EXT) :: String
  parse(parsetags(data), partial = partial)
end

function parse(input::String; partial::Bool = true) :: String
  parsehtml(input, partial = partial)
end

"""
    parsetags(line::Tuple{Int,String}, strip_close_tag = false) :: String

Parses special HTML+Julia tags.
"""
function parsetags(line::Tuple{Int,String}) :: String
  parsetags(line[2])
end

function parsetags(code::String) :: String
  replace(
    replace(code, "<%"=>"""<script type="julia/eval">"""),
    "%>"=>"""</script>""")
end

"""
    register_elements() :: Nothing

Generated functions that represent Julia functions definitions corresponding to HTML elements.
"""
function register_elements(; context = @__MODULE__) :: Nothing
  for elem in NORMAL_ELEMENTS
    register_normal_element(elem)
  end

  for elem in VOID_ELEMENTS
    register_void_element(elem)
  end

  for elem in CUSTOM_ELEMENTS
    Core.eval(@__MODULE__, """include("html/$elem.jl")""" |> Meta.parse)
    elem in NON_EXPORTED || Core.eval(context, "export $elem" |> Meta.parse)
  end

  for elem in SVG_ELEMENTS
    register_normal_element(elem)
  end

  nothing
end

"""
    register_element(elem::Union{Symbol,String}, elem_type::Union{Symbol,String} = :normal; context = @__MODULE__) :: Nothing

Generates a Julia function representing an HTML element.
"""
function register_element(elem::Union{Symbol,String}, elem_type::Union{Symbol,String} = :normal; context = @__MODULE__) :: Nothing
  elem = string(elem)
  occursin("-", elem) && (elem = denormalize_element(elem))

  elem_type == :normal ? register_normal_element(elem) : register_void_element(elem)
end

"""
    register_normal_element(elem::Union{Symbol,String}; context = @__MODULE__) :: Nothing

Generates a Julia function representing a "normal" HTML element: that is an element with a closing tag, <tag>...</tag>
"""
function register_normal_element(elem::Union{Symbol,String}; context = @__MODULE__) :: Nothing
  Core.eval(context, """
    function $elem(f::Function, args...; attrs...) :: HTMLString
      \"\"\"\$(normal_element(f, "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  Core.eval(context, """
    function $elem(children::Union{String,Vector{String}} = "", args...; attrs...) :: HTMLString
      \"\"\"\$(normal_element(children, "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  Core.eval(context, """
    function $elem(children::Any, args...; attrs...) :: HTMLString
      \"\"\"\$(normal_element(string(children), "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  Core.eval(context, """
    function $elem(children::Vector{Any}, args...; attrs...) :: HTMLString
      \"\"\"\$(normal_element([string(c) for c in children], "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  elem in NON_EXPORTED || Core.eval(context, "export $elem" |> Meta.parse)

  nothing
end

"""
    register_void_element(elem::Union{Symbol,String}; context::Module = @__MODULE__) :: Nothing

Generates a Julia function representing a "void" HTML element: that is an element without a closing tag, <tag />
"""
function register_void_element(elem::Union{Symbol,String}; context::Module = @__MODULE__) :: Nothing
  Core.eval(context, """
    function $elem(args...; attrs...) :: HTMLString
      \"\"\"\$(void_element("$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  elem in NON_EXPORTED || Core.eval(context, "export $elem" |> Meta.parse)

  nothing
end

"""
    for_each(f::Function, v)

Iterates over the `v` Vector and applies function `f` for each element.
The results of each iteration are concatenated and the final string is returned.
"""
function for_each(f::Function, v)
  [f(x) for x in v] |> join
end

"""
    collection(template::Function, collection::Vector{T})::String where {T}

Creates a view fragment by repeateadly applying a function to each element of the collection.
"""
function collection(template::Function, collection::Vector{T})::String where {T}
  for_each(collection) do item
    template(item) |> string
  end
end

### === ###
### EXCEPTIONS ###

"""
    Genie.Router.error(error_message::String, ::Type{MIME"text/html"}, ::Val{500}; error_info::String = "") :: HTTP.Response

Returns a 500 error response as an HTML doc.
"""
function Genie.Router.error(error_message::String, ::Type{MIME"text/html"}, ::Val{500}; error_info::String = "") :: HTTP.Response
  serve_error_file(500, error_message, error_info = error_info)
end

"""
    Genie.Router.error(error_message::String, ::Type{MIME"text/html"}, ::Val{404}; error_info::String = "") :: HTTP.Response

Returns a 404 error response as an HTML doc.
"""
function Genie.Router.error(error_message::String, ::Type{MIME"text/html"}, ::Val{404}; error_info::String = "") :: HTTP.Response
  serve_error_file(404, error_message, error_info = error_info)
end

"""
    Genie.Router.error(error_code::Int, error_message::String, ::Type{MIME"text/html"}; error_info::String = "") :: HTTP.Response

Returns an error response as an HTML doc.
"""
function Genie.Router.error(error_code::Int, error_message::String, ::Type{MIME"text/html"}; error_info::String = "") :: HTTP.Response
  serve_error_file(error_code, error_message, error_info = error_info)
end

"""
    serve_error_file(error_code::Int, error_message::String = "", params::Dict{Symbol,Any} = Dict{Symbol,Any}()) :: Response

Serves the error file correspoding to `error_code` and current environment.
"""
function serve_error_file(error_code::Int, error_message::String = ""; error_info::String = "") :: HTTP.Response
  page_code = error_code in [404, 500] ? "$error_code" : "xxx"

  try
    error_page_file = isfile(joinpath(Genie.config.server_document_root, "error-$page_code.html")) ?
                        joinpath(Genie.config.server_document_root, "error-$page_code.html") :
                          joinpath(@__DIR__, "..", "..", "files", "static", "error-$page_code.html")

    error_page =  open(error_page_file) do f
                    read(f, String)
                  end

    if error_code == 500
      error_page = replace(error_page, "<error_description/>"=>split(error_message, "\n")[1])

      error_message = if Genie.Configuration.isdev()
                      """$("#" ^ 25) ERROR STACKTRACE $("#" ^ 25)\n$error_message                                     $("\n" ^ 3)""" *
                      """$("#" ^ 25)  REQUEST PARAMS  $("#" ^ 25)\n$(Millboard.table(Genie.Router.params()))                        $("\n" ^ 3)""" *
                      """$("#" ^ 25)     ROUTES       $("#" ^ 25)\n$(Millboard.table(Genie.Router.named_routes() |> Dict))  $("\n" ^ 3)""" *
                      """$("#" ^ 25)    JULIA ENV     $("#" ^ 25)\n$ENV                                               $("\n" ^ 1)"""
      else
        ""
      end

      error_page = replace(error_page, "<error_message/>"=>escapeHTML(error_message))

    elseif error_code == 404
      error_page = replace(error_page, "<error_message/>"=>error_message)

    else
      error_page = replace(replace(error_page, "<error_message/>"=>error_message), "<error_info/>"=>error_info)
    end

    HTTP.Response(error_code, ["Content-Type"=>"text/html"], body = error_page)
  catch ex
    @error ex
    HTTP.Response(error_code, ["Content-Type"=>"text/html"], body = "Error $page_code: $error_message")
  end
end

### === ###

"""
    @yield

Outputs the rendering of the view within the template.
"""
macro yield()
  :(view!()() |> join)
end
macro yield(value)
  :(view!($value))
end

function view!()
  haskey(task_local_storage(), :__yield) ? task_local_storage(:__yield) : task_local_storage(:__yield, String[])
end

function view!(value)
  task_local_storage(:__yield, value)
end

function el(; vars...)
  OrderedCollections.OrderedDict(vars)
end

register_elements()

end


text:
@testset "Proposals.Uniform -> Proposals.Rejection deprecation" begin
    prop = @test_deprecated Proposals.Uniform()
    @test prop === Proposals.Rejection()
end

ext: jl


text:
__precompile__()

module OffsetArrays

Base.@deprecate_binding (..) Colon()

using Base: Indices, LinearSlow, LinearFast, tail

export OffsetArray, @unsafe

immutable OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
end
typealias OffsetVector{T,AA<:AbstractArray} OffsetArray{T,1,AA}

OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) =
    OffsetArray{T,N,typeof(A)}(A, offsets)
OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) =
    OffsetArray(A, offsets)

(::Type{OffsetArray{T,N}}){T,N}(inds::Indices{N}) =
    OffsetArray{T,N,Array{T,N}}(Array{T,N}(map(length, inds)), map(indexoffset, inds))
(::Type{OffsetArray{T}}){T,N}(inds::Indices{N}) = OffsetArray{T,N}(inds)
OffsetArray{T}(A::AbstractArray{T,0}) = OffsetArray{T,0,typeof(A)}(A, ())
OffsetArray{T,N}(::Type{T}, inds::Vararg{UnitRange{Int},N}) = OffsetArray{T,N}(inds)

# The next two are necessary for ambiguity resolution. Really, the
# second method should not be necessary.
OffsetArray{T}(A::AbstractArray{T,0}, inds::Tuple{}) = OffsetArray{T,0,typeof(A)}(A, ())
OffsetArray{T,N}(A::AbstractArray{T,N}, inds::Tuple{}) = error("this should never be called")
function OffsetArray{T,N}(A::AbstractArray{T,N}, inds::NTuple{N,AbstractUnitRange})
    lA = map(length, indices(A))
    lI = map(length, inds)
    lA == lI || throw(DimensionMismatch("supplied indices do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
    OffsetArray(A, map(indexoffset, inds))
end
OffsetArray{T,N}(A::AbstractArray{T,N}, inds::Vararg{AbstractUnitRange,N}) =
    OffsetArray(A, inds)

Base.linearindexing{T<:OffsetArray}(::Type{T}) = Base.linearindexing(parenttype(T))
parenttype{T,N,AA}(::Type{OffsetArray{T,N,AA}}) = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

errmsg(A) = error("size not supported for arrays with indices $(indices(A)); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/")
Base.size(A::OffsetArray) = errmsg(A)
Base.size(A::OffsetArray, d) = errmsg(A)
Base.eachindex(::LinearSlow, A::OffsetArray) = CartesianRange(indices(A))
Base.eachindex(::LinearFast, A::OffsetVector) = indices(A, 1)

# Implementations of indices and indices1. Since bounds-checking is
# performance-critical and relies on indices, these are usually worth
# optimizing thoroughly.
@inline Base.indices(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? indices(parent(A))[d] + A.offsets[d] : (1:1)
@inline Base.indices(A::OffsetArray) = _indices(indices(parent(A)), A.offsets)  # would rather use ntuple, but see #15276
@inline _indices(inds, offsets) = (inds[1]+offsets[1], _indices(tail(inds), tail(offsets))...)
_indices(::Tuple{}, ::Tuple{}) = ()
Base.indices1{T}(A::OffsetArray{T,0}) = 1:1  # we only need to specialize this one

function Base.similar{T}(A::OffsetArray, ::Type{T}, dims::Dims)
    B = similar(parent(A), T, dims)
end
function Base.similar{T}(A::AbstractArray, ::Type{T}, inds::Tuple{UnitRange,Vararg{UnitRange}})
    B = similar(A, T, map(length, inds))
    OffsetArray(B, map(indexoffset, inds))
end

Base.similar(f::Union{Function,DataType}, shape::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(f(map(length, shape)), map(indexoffset, shape))

Base.reshape(A::AbstractArray, inds::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(reshape(A, map(length, inds)), map(indexoffset, inds))

Base.reshape(A::OffsetArray, inds::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(reshape(parent(A), map(length, inds)), map(indexoffset, inds))

function Base.reshape(A::OffsetArray, inds::Tuple{UnitRange,Vararg{Union{UnitRange,Int,Base.OneTo}}})
    throw(ArgumentError("reshape must supply UnitRange indices, got $(typeof(inds)).\n       Note that reshape(A, Val{N}) is not supported for OffsetArrays."))
end

# Don't allow bounds-checks to be removed during Julia 0.5
@inline function Base.getindex{T,N}(A::OffsetArray{T,N}, I::Vararg{Int,N})
    checkbounds(A, I...)
    @inbounds ret = parent(A)[offset(A.offsets, I)...]
    ret
end
@inline function Base._getindex(::LinearFast, A::OffsetVector, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[offset(A.offsets, (i,))[1]]
    ret
end
@inline function Base._getindex(::LinearFast, A::OffsetArray, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[i]
    ret
end
@inline function Base.setindex!{T,N}(A::OffsetArray{T,N}, val, I::Vararg{Int,N})
    checkbounds(A, I...)
    @inbounds parent(A)[offset(A.offsets, I)...] = val
    val
end
@inline function Base._setindex!(::LinearFast, A::OffsetVector, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[offset(A.offsets, (i,))[1]] = val
    val
end
@inline function Base._setindex!(::LinearFast, A::OffsetArray, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[i] = val
    val
end

### Convenience functions ###

Base.fill(x, inds::Tuple{UnitRange,Vararg{UnitRange}}) = fill!(OffsetArray{typeof(x)}(inds), x)
@inline Base.fill(x, ind1::UnitRange, inds::UnitRange...) = fill(x, (ind1, inds...))

### Low-level utilities ###

# Computing a shifted index (subtracting the offset)
offset{N}(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) = _offset((), offsets, inds)
_offset(out, ::Tuple{}, ::Tuple{}) = out
@inline _offset(out, offsets, inds) = _offset((out..., inds[1]-offsets[1]), Base.tail(offsets), Base.tail(inds))

indexoffset(r::Range) = first(r) - 1
indexoffset(i::Integer) = 0

macro unsafe(ex)
    esc(unsafe(ex))
end
unsafe(ex) = ex
function unsafe(ex::Expr)
    if ex.head ∈ (:+=, :-=, :*=, :/=)
        ex = Expr(:(=), ex.args[1], Expr(:call, Symbol(string(ex.head)[1]), ex.args...))
    end
    if ex.head == :(=)
        a = ex.args[1]
        if isa(a, Expr) && (a::Expr).head == :ref
            # setindex!
            newargs = map(unsafe, ex.args[2:end])
            @assert length(newargs) == 1
            return Expr(:call, :(OffsetArrays.unsafe_setindex!), (a::Expr).args[1], newargs[1], (a::Expr).args[2:end]...)
        end
    end
    newargs = map(unsafe, ex.args)
    if ex.head == :ref
        # getindex
        return Expr(:call, :(OffsetArrays.unsafe_getindex), newargs...)
    end
    Expr(ex.head, newargs...)
end

@inline unsafe_getindex(a::AbstractArray, I...) = (@inbounds ret = a[I...]; ret)
@inline unsafe_setindex!(a::AbstractArray, val, I...) = (@inbounds a[I...] = val; val)

@inline unsafe_getindex(a::OffsetArray, I::Int...) = unsafe_getindex(parent(a), offset(a.offsets, I)...)
@inline unsafe_setindex!(a::OffsetArray, val, I::Int...) = unsafe_setindex!(parent(a), val, offset(a.offsets, I)...)
@inline unsafe_getindex(a::OffsetArray, I...) = unsafe_getindex(a, Base.IteratorsMD.flatten(I)...)
@inline unsafe_setindex!(a::OffsetArray, val, I...) = unsafe_setindex!(a, val, Base.IteratorsMD.flatten(I)...)

# Indexing a SubArray which has OffsetArray indices
typealias OffsetSubArray{T,N,P,I<:Tuple{OffsetArray,Vararg{OffsetArray}}} SubArray{T,N,P,I,false}
@inline function unsafe_getindex{T,N}(a::OffsetSubArray{T,N}, I::Vararg{Int,N})
    J = map(unsafe_getindex, a.indexes, I)
    unsafe_getindex(parent(a), J...)
end
@inline function unsafe_setindex!{T,N}(a::OffsetSubArray{T,N}, val, I::Vararg{Int,N})
    J = map(unsafe_getindex, a.indexes, I)
    unsafe_setindex!(parent(a), val, J...)
end
@inline unsafe_getindex(a::OffsetSubArray, I::Union{Integer,CartesianIndex}...) = unsafe_getindex(a, Base.IteratorsMD.flatten(I)...)
@inline unsafe_setindex!(a::OffsetSubArray, val, I::Union{Integer,CartesianIndex}...) = unsafe_setindex!(a, val, Base.IteratorsMD.flatten(I)...)

# Deprecations
import Base: zeros, ones
@deprecate zeros(T::Type, ind1::UnitRange, ind2::UnitRange, inds::UnitRange...) fill!(OffsetArray{T}((ind1, ind2, inds...)), zero(T))
@deprecate  ones(T::Type, ind1::UnitRange, inds2::UnitRange, inds::UnitRange...) fill!(OffsetArray{T}((ind1, ind2, inds...)), one(T))
@deprecate zeros(ind1::UnitRange, ind2::UnitRange, inds::UnitRange...) fill!(OffsetArray{Float64}((ind1, ind2, inds...)), 0)
@deprecate  ones(ind1::UnitRange, ind2::UnitRange, inds::UnitRange...) fill!(OffsetArray{Float64}((ind1, ind2, inds...)), 1)

end # module


text:
"Get the WKT of an Integer EPSG code"
function epsg2wkt(epsgcode::Int)
    srs = ArchGDAL.GDAL.osrnewspatialreference(C_NULL)
    ArchGDAL.GDAL.osrimportfromepsg(srs, epsgcode)
    wkt_ptr = Ref(Cstring(C_NULL))
    ArchGDAL.GDAL.osrexporttowkt(srs, wkt_ptr)
    return unsafe_string(wkt_ptr[])
end

"Get the WKT of an Proj string"
function proj2wkt(projstring::AbstractString)
    srs = ArchGDAL.GDAL.osrnewspatialreference(C_NULL)
    ArchGDAL.GDAL.osrimportfromproj4(srs, projstring)
    wkt_ptr = Ref(Cstring(C_NULL))
    ArchGDAL.GDAL.osrexporttowkt(srs, wkt_ptr)
    return unsafe_string(wkt_ptr[])
end

function wkt2wkt(wktstring::AbstractString)
    srs = ArchGDAL.GDAL.osrnewspatialreference(C_NULL)
    ArchGDAL.GDAL.osrimportfromwkt(srs, [wktstring])
    wkt_ptr = Ref(Cstring(C_NULL))
    ArchGDAL.GDAL.osrexporttowkt(srs, wkt_ptr)
    return unsafe_string(wkt_ptr[])
end

"""Parse CRS string into WKT."""
function str2wkt(crs_string::AbstractString)
    if startswith(crs_string, "+proj=")
        return proj2wkt(crs_string)
    elseif startswith(crs_string, "EPSG:")
        epsg_code = parse(Int, crs_string[findlast("EPSG:", crs_string).stop+1:end])
        return epsg2wkt(epsg_code)
    else
        # Fallback method to validate string
        wkt =  wkt2wkt(crs_string)
        return wkt
    end
end

"Set CRS on GeoArray by epsgcode"
epsg!(ga::GeoArray, epsgcode::Int) = crs!(ga, EPSG(epsgcode))
epsg!(ga::GeoArray, epsgstring::AbstractString) = crs!(ga, EPSG(epsgstring))

function crs!(ga::GeoArray, crs::GFT.CoordinateReferenceSystemFormat)
    ga.crs = convert(GFT.WellKnownText, GFT.CRS(), crs)
    ga
end


text:
include("lib_for_robot.jl")
function task19(r)
    steps = 1
    side = Nord
    while ismarker(r)==false
        for _ in 1:2
            special_move!(r,side,steps)
            side=counterclockwise_side(side)
        end
        steps += 1
    end
    println("Marker found!")
end


text:
"""
    closeness_centrality(g, distmx=weights(g); normalize=true)

Calculate the [closeness centrality](https://en.wikipedia.org/wiki/Centrality#Closeness_centrality)
of the graph `g`. Return a vector representing the centrality calculated for each node in `g`.

### Optional Arguments
- `normalize=true`: If true, normalize the centrality value of each
node `n` by ``\\frac{|δ_n|}{|V|-1}``, where ``δ_n`` is the set of vertices reachable
from node `n`.

# Examples
```jldoctest
julia> using LightGraphs

julia> closeness_centrality(StarGraph(5))
5-element Array{Float64,1}:
 1.0
 0.5714285714285714
 0.5714285714285714
 0.5714285714285714
 0.5714285714285714

julia> closeness_centrality(PathGraph(4))
4-element Array{Float64,1}:
 0.5
 0.75
 0.75
 0.5
```
"""
function closeness_centrality(g::AbstractGraph,
    distmx::AbstractMatrix=weights(g);
    normalize=true)

    n_v = nv(g)
    closeness = zeros(n_v)

    for u in vertices(g)
        if degree(g, u) == 0     # no need to do Dijkstra here
            closeness[u] = 0.0
        else
            d = dijkstra_shortest_paths(g, u, distmx).dists
            δ = filter(x -> x != typemax(x), d)
            σ = sum(δ)
            l = length(δ) - 1
            if σ > 0
                closeness[u] = l / σ
                if normalize
                    n = l / (n_v - 1)
                    closeness[u] *= n
                end
            end
        end
    end
    return closeness
end


text:
#===============================================================================

Let r be the remainder when (a−1)^n + (a+1)^n is divided by a^2.

For example, if a = 7 and n = 3, then r = 42: 6^3 + 8^3 = 728 ≡ 42 mod 49. And
as n varies, so too will r, but for a = 7 it turns out that r_max = 42.

For 3 ≤ a ≤ 1000, find ∑ r_max.

--------------------------------------------------------------------------------

Let r(a,n) := ((a−1)^n + (a+1)^n) % a^2
It comes: r(a, n) = 0 for a=1 (because of the modulo 1)
                  = 2 for n even
                  = 2an % a^2 = a(2n % a) for n odd

          max_n(2n % a) = a-1 if a is odd
                        = a-2 if a is even

This max_n(r(a, n)) = 0 for a = 1
                    = 2 for a = 2
                    = a(a-1) for 3 <= a and a odd
                    = a(a-2) for 3 <= a and a even

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

rmax(a) = isodd(a) ? a*(a-1) : a*(a-2) # It is assumed that a>=3

function compute(limit)
    result = 0
    for a in 3:limit
        result += rmax(a)
    end
    result;
end

println("Euler 120: $(compute(1000))")

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
using NamedArrays

export create_z

function create_z(model::EoS, z::AbstractArray)
    return NamedArray(z, model.components)
end

"""
    normalize_units(val)
On normal numbers, it is the identity, but on numbers or vectors of `Unitful.Quantity`,it converts the unit to an equivalent SI unit and strips the unit information.

# Examples
```julia-repl
julia> normalize_units(0.0u"°C")
273.15

```
```julia-repl
julia> normalize_units(273.15)
273.15
```
"""
normalize_units(x) =  Unitful.ustrip.(Unitful.upreferred.(x))


text:
#
# Denoise an S2-valued image Example with Cyclic Proximal Point applied to the
#
# L2-TV functional with anisotropic TV
#
# where the example is the same data as for the corresponding CP algorithm
#
using Manopt, Manifolds
using Images, CSV, DataFrames, LinearAlgebra, JLD2
#
# Settings
experiment_name = "S2_Whirl_CPPA"
results_folder = joinpath(@__DIR__, "S2_TV")
export_result = true
export_orig = true
export_table = true
export_function_value = true
asymptote_render_detail = 2
!isdir(results_folder) && mkdir(results_folder)
#
# Manifold & Data
f = artificial_S2_whirl_image(64)
pixelM = Sphere(2);

if export_orig
    orig_filename = joinpath(results_folder, experiment_name * "-orig.asy")
    asymptote_export_S2_data(orig_filename; data=f)
    render_asymptote(orig_filename; render=asymptote_render_detail)
end
#
# Parameters
α = 1.5
maxIterations = 4000
#
# Build Problem for L2-TV
M = PowerManifold(pixelM, NestedPowerRepresentation(), size(f)...)
d = length(size(f))
iRep = [Integer.(ones(d))..., d]
fidelity(M, x) = 1 / 2 * distance(M, x, f)^2
Λ(M, x) = forward_logs(M, x) # on T_xN
prior(M, x) = norm(norm.(Ref(pixelM), repeat(x, iRep...), Λ(M, x)), 1)
#
# Setup and Optimize
cost(M, x) = fidelity(M, x) + α * prior(M, x)
proxes = ((M, λ, x) -> prox_distance(M, λ, f, x, 2), (M, λ, x) -> prox_TV(M, α * λ, x, 1))
x0 = f
@time o = cyclic_proximal_point(
    M,
    cost,
    proxes,
    x0;
    debug=[
        :Iteration,
        " | ",
        DebugProximalParameter(),
        " | ",
        :Change,
        " | ",
        :Cost,
        "\n",
        100,
        :Stop,
    ],
    record=[:Iteration, :Iterate, :Cost],
    stopping_criterion=StopAfterIteration(maxIterations),
    λ=i -> π / (2 * i),
    return_options=true,
)
y = get_solver_result(o)
yRec = get_record(o)
#
# Results
if export_result
    result_filename = joinpath(
        results_folder,
        experiment_name * "-result-$(maxIterations)-α$(replace(string(α), "." => "-")).asy",
    )
    asymptote_export_S2_data(result_filename; data=y)
    render_asymptote(result_filename; render=asymptote_render_detail)
end
if export_table
    A = cat([y[1] for y in yRec], [y[3] for y in yRec]; dims=2)
    table_filename = joinpath(results_folder, experiment_name * "-recorded-cost.csv")
    CSV.write(table_filename, DataFrame(A); writeheader=false)
end
if export_function_value
    fctval_filename = joinpath(results_folder, experiment_name * "-cost.jld2")
    values = Dict("cost_function_value" => last(yRec)[3], "iterations" => length(yRec) - 1)
    save(fctval_filename, values)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using TOML, Test
using TOML: ParserError

@testset "TOML.(try)parse(file) entrypoints" begin
    dict = Dict{String,Any}("a" => 1)
    str = "a = 1"; invalid_str = "a"
    path, io = mktemp(); write(io, str); close(io)
    invalid_path, io = mktemp(); write(io, invalid_str); close(io)
    p = TOML.Parser()
    # TOML.parse
    @test TOML.parse(str) == TOML.parse(SubString(str)) ==
          TOML.parse(IOBuffer(str)) ==
          TOML.parse(p, str) == TOML.parse(p, SubString(str)) ==
          TOML.parse(p, IOBuffer(str)) == dict
    @test_throws ParserError TOML.parse(invalid_str)
    @test_throws ParserError TOML.parse(SubString(invalid_str))
    @test_throws ParserError TOML.parse(IOBuffer(invalid_str))
    @test_throws ParserError TOML.parse(p, invalid_str)
    @test_throws ParserError TOML.parse(p, SubString(invalid_str))
    @test_throws ParserError TOML.parse(p, IOBuffer(invalid_str))
    # TOML.tryparse
    @test TOML.tryparse(str) == TOML.tryparse(SubString(str)) ==
          TOML.tryparse(IOBuffer(str)) ==
          TOML.tryparse(p, str) == TOML.tryparse(p, SubString(str)) ==
          TOML.tryparse(p, IOBuffer(str)) == dict
    @test TOML.tryparse(invalid_str) isa ParserError
    @test TOML.tryparse(SubString(invalid_str)) isa ParserError
    @test TOML.tryparse(IOBuffer(invalid_str)) isa ParserError
    @test TOML.tryparse(p, invalid_str) isa ParserError
    @test TOML.tryparse(p, SubString(invalid_str)) isa ParserError
    @test TOML.tryparse(p, IOBuffer(invalid_str)) isa ParserError
    # TOML.parsefile
    @test TOML.parsefile(path) == TOML.parsefile(SubString(path)) ==
          TOML.parsefile(p, path) == TOML.parsefile(p, SubString(path)) == dict
    @test_throws ParserError TOML.parsefile(invalid_path)
    @test_throws ParserError TOML.parsefile(SubString(invalid_path))
    @test_throws ParserError TOML.parsefile(p, invalid_path)
    @test_throws ParserError TOML.parsefile(p, SubString(invalid_path))
    @test_throws ErrorException TOML.parsefile(homedir())
    @test_throws ErrorException TOML.parsefile(p, homedir())
    # TOML.tryparsefile
    @test TOML.tryparsefile(path) == TOML.tryparsefile(SubString(path)) ==
          TOML.tryparsefile(p, path) == TOML.tryparsefile(p, SubString(path)) == dict
    @test TOML.tryparsefile(invalid_path) isa ParserError
    @test TOML.tryparsefile(SubString(invalid_path)) isa ParserError
    @test TOML.tryparsefile(p, invalid_path) isa ParserError
    @test TOML.tryparsefile(p, SubString(invalid_path)) isa ParserError
    @test_throws ErrorException TOML.tryparsefile(homedir())
    @test_throws ErrorException TOML.tryparsefile(p, homedir())
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
module AllShapesTest

using  ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

filenameTiger = joinpath(Modia3D.path, "objects", "fish", "SiameseTiger0.3ds")
filenameBunny = joinpath(Modia3D.path, "objects", "bunny", "bunny.obj")

# Properties
mat1 = VisualMaterial(color="Green", transparency=0.5)
mat2 = "Red"

font2 = Font(fontFamily="FreeSans", bold=false, italic=false, charSize=0.1,
    color = "LightBlue", transparency=0.0)

# material for massProperties where all MassProperties are defined
m   = 2.0
rCM = [1,2,3]
massDefined = MassProperties(mass=m, centerOfMass=rCM, Ixx=3.0, Iyy=4.0, Izz=4.0, Ixy=6.0, Ixz=7.0, Iyz=9.0)

# mass properties are computed from shape and mass
massAndGeo = MassPropertiesFromShapeAndMass(mass=4.5)

# mass properties are computed from shape and solid material
geomAndSolidMat = MassPropertiesFromShape()

# Objects3D
AllShapes = Model(
    gravField = UniformGravityField(g=9.81, n=[0, 0,-1]),
    world = Object3D(feature=Scene(gravityField=:gravField)),
    worldFrame = Object3D(parent=:world,feature=Visual(shape=CoordinateSystem())),

    # visual elements
    sphere = Object3D(parent=:world, feature=Visual(shape=Sphere(diameter=0.5)),
        translation=[5.5, 0.0, 2.5]),
    ellipsoid = Object3D(parent=:world, feature=Visual(shape =
        Ellipsoid(lengthX=0.6,lengthY=0.9, lengthZ=0.3)), translation=[0.0, 0.0, 2.5]),
    box = Object3D(parent=:world, feature=Visual(shape=Box(lengthX=0.9, lengthY=0.5,
        lengthZ=0.3)), translation=[3.5, 0.0, 2.5]),
    cylinder = Object3D(parent=:world, feature=Visual(shape=Cylinder(diameter=0.5,
        length=0.8)), translation=[1.5, 0.0, 2.5]),
    pipe = Object3D(parent=:world, feature=Visual(shape=Cylinder(diameter=0.5,
        length=0.8, innerDiameter=0.45)), translation=[-4.5, 0.0, 2.5]),
    capsule = Object3D(parent=:world, feature=Visual(shape=Capsule(diameter=0.4,
        length=0.45)), translation=[-1.5, 0.0, 2.5]),
    beam = Object3D(parent=:world, feature=Visual(shape=Beam(axis=1, length=0.4, width=0.5, thickness=0.3)),
        translation=[ -2.6, 0.0, 2.5]),
    cone = Object3D(parent=:world, feature=Visual(shape=Cone(axis=3, diameter=0.3, length=0.7, topDiameter=0.15)),
        translation=[ -3.6, 0.0, 2.5], rotation=Modia3D.rot1(-pi/3)),
    bunny = Object3D(parent=:world, feature=Visual(shape =
        FileMesh(filename = filenameBunny, scale=[0.1, 0.1, 0.1])), translation=[-6.0, 0.0, 2.5]),
    bunnyPartition = Object3D(parent=:world, feature = Visual(shape =
        FileMesh(filename = filenameBunny, scale=[0.1, 0.1, 0.1], convexPartition=true)), translation=[-8.0, 0.0, 2.5]),

    # Solid objects
    sphereSolid = Object3D(parent=:world, feature=Solid(shape=Sphere(diameter=0.5),
        visualMaterial=mat1, solidMaterial="Steel"), translation=[5.5, 0.0, 0.0]),

    ellipsoidSolid = Object3D(
        feature=Solid(shape=Ellipsoid(lengthX=0.6, lengthY=0.9, lengthZ=0.3), visualMaterial=mat1, solidMaterial="DryWood")),
    free = FreeMotion(obj1=:world, obj2=:ellipsoidSolid),

    boxSolid = Object3D(
        feature=Solid(shape=Box(lengthX=0.9, lengthY=0.5, lengthZ=0.3), visualMaterial=mat1, solidMaterial="Steel")),
    fixed = Fix(obj1=:world, obj2=:boxSolid, translation=[3.5, 0.0, 0.0]),

    cylinderSolid = Object3D(parent=:world, feature=Solid(shape=Cylinder(diameter=0.5,
        length=0.8), visualMaterial=mat1, solidMaterial="Steel"), translation=[1.5, 0.0, 0.0]),
    pipeSolid = Object3D(parent=:world, feature=Solid(shape=Cylinder(diameter=0.5,
        length=0.8, innerDiameter=0.45), visualMaterial=mat1, solidMaterial="Steel", massProperties=massDefined), translation=[-4.5, 0.0, 0.0]),

    capsuleSolid = Object3D(parent=:world, feature=Solid(shape=Capsule(diameter=0.4,
        length=0.45), visualMaterial=mat1, massProperties=massAndGeo), translation=[-1.5, 0.0, 0.0]),

    beamSolid = Object3D(parent=:world, feature=Solid(shape=Beam(axis=1, length=0.4, width=0.5, thickness=0.3),visualMaterial=mat1, solidMaterial="Steel", massProperties=geomAndSolidMat), translation =[ -2.6, 0.0, 0.0]),

    coneSolid = Object3D(parent=:world, feature=Solid(shape=Cone(axis=3, diameter=0.3, length=0.7, topDiameter=0.15),
        visualMaterial=mat1, solidMaterial="Steel", collision=true, massProperties = geomAndSolidMat), translation=[ -3.6, 0.0, 0.0], rotation=Modia3D.rot1(-pi/3)),
    bunnySolid = Object3D(parent=:world, feature=Solid(shape =
        FileMesh(filename=filenameBunny, scale=[0.1, 0.1, 0.1]), solidMaterial="Steel", massProperties = geomAndSolidMat ), translation=[-6.0, 0.0, 0.0]),
    bunnyPartitionSolid = Object3D(parent=:world, feature=Solid(shape =
        FileMesh(filename = filenameBunny, scale=[0.1, 0.1, 0.1], convexPartition=true), solidMaterial="Steel", collision=true), translation=[-8.0, 0.0, 0.0]),

    # pure visual objects
    spring    = Object3D(parent=:world, feature=Visual(shape =
        SpringShape(axis=3, length=0.8, diameter=0.3, wireDiameter=0.04, windings=5), visualMaterial=mat2), translation=[3.0, 0.0, -2.5]),
    gearWheel = Object3D(parent=:world, feature=Visual(shape =
        GearWheel(axis=3, diameter=0.5, length=0.8, innerDiameter=0.25, angle=15/180*pi, teeth=16), visualMaterial=mat2), translation=[4.5, 0.0, -2.5]),
    grid      = Object3D(parent=:world, feature=Visual(shape =
        Grid(axis=3, length=1.0, width=0.6, distance=0.05, lineWidth=1.0)), translation=[0.0, 0.0, -2.5]),
    coSystem  = Object3D(parent=:world, feature = Visual(shape =
        CoordinateSystem(length=0.7)), translation=[-1.5, 0.0, -2.5]),
    fileMeshTiger = Object3D(parent=:world, feature=Visual(shape =
        FileMesh(filename=filenameTiger, scale=[4.0, 4.0, 2.5])), translation=[-3.0, 0.0, -2.5]),
    # Text object
    text = Object3D(parent=:world, feature = Visual(shape =
    TextShape(text="Hello world!")), translation=[1.0, 0.0, -2.5]),
)

allShapes = @instantiateModel(buildModia3D(AllShapes), unitless=true, log=false, logStateSelection=false, logCode=false)

simulate!(allShapes, stopTime=0.0, log=false, logStates=false, logEvents=false)

println("... Basics/AllShapes.jl completed.")

end


text:
@testset "SimpleVCE" begin
    @test sprint(show, MIME("text/plain"), SimpleVCE()) ==
        "Simple variance-covariance estimator"
end

@testset "HRVCE" begin
    @test criticalvalue(HRVCE(), 0.95, 100) ≈ 1.9599639845400576
    @test pvalue(HRVCE(), 1.96, 100) ≈ 0.04999579029644087
    @test sprint(show, MIME("text/plain"), HRVCE()) ==
        "Heteroskedasticity-robust variance-covariance estimator"
end

@testset "EWC" begin
    @test EWC() === EqualWeightedCosine()
    @test sprint(show, MIME("text/plain"), EWC()) == "Equal-weighted cosine transform"
    v = HARVCE(EWC())
    @test criticalvalue(v, 0.95, 100) ≈ 2.262157162798205
    @test pvalue(v, 1.96, 100) ≈ 0.08164440546041651
end

@testset "HARVCE" begin
    v = HARVCE(EWC())
    @test sprint(show, MIME("text/plain"), v) == """
        Heteroskedasticity-autocorrelation-robust variance-covariance estimator:
          Long-run variance estimator: Equal-weighted cosine transform"""
end


text:
verbose = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using NeidSolarScripts
 using NeidSolarScripts.SolarRotation
 using NeidSolarScripts.DifferentialExtinction
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, Dates
 #using StatsBase, Statistics, Dates

fits_target_str = "Sun"
 if !@isdefined(target_subdir)
      target_subdir = ""
 end
 if !@isdefined(create_missing_continuum_files)
     create_missing_continuum_files = false
 end
 paths_to_search_for_param = [pwd(),pkgdir(NeidSolarScripts)]

 if verbose println("# Finding what data files are avaliable.")  end
 eval(read_data_paths(paths_to_search=paths_to_search_for_param))
 @assert isdefined(Main,:neid_data_path)
 @assert isdefined(Main,:output_dir)
 output_path = joinpath(neid_data_path, target_subdir,output_dir)
 manifest_filename = joinpath(output_path,"manifest.csv")
 manifest_calib_filename = joinpath(output_path,"manifest_calib.csv")

can_skip_generating_manifest = false
if isfile(manifest_filename) || islink(manifest_filename)
    df_files  = CSV.read(manifest_filename, DataFrame)
    warn_link =  islink(manifest_filename) ? " (a link)" : ""
    println("# Read ", manifest_filename, warn_link, ".")
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
    @assert hasproperty(df_files,:alt_sun)
    @assert hasproperty(df_files,:Δv_diff_ext)
    @assert hasproperty(df_files,:Δfwhm²)
    @assert hasproperty(df_files,:order_snrs)
    if eltype(df_files[!,:order_snrs]) == String
        df_files[!,:order_snrs] = map(i->parse.(Float64,split(df_files[i,:order_snrs][2:end-1],',')),1:size(df_files,1))
    end
    @assert eltype(df_files[!,:order_snrs]) == Vector{Float64}
    println("# Required fields present ", size(df_files), ".  No need to regenerate manifest.")
    can_skip_generating_manifest = true
    if any(isnan.(df_files[!,:order_snrs])) ||
       any(isnan.(df_files[!,:Δfwhm²]))    ||
       any(isnan.(df_files[!,:Δv_diff_ext]))
          can_skip_generating_manifest = false
    end
else
    println("# Will need to generate manifest at ", manifest_filename, ".")
end

#can_skip_generating_manifest = false
if can_skip_generating_manifest && !create_missing_continuum_files
    exit()
end

df_files = NEID.make_manifest(joinpath(neid_data_path, target_subdir))

calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]
df_files_calib = df_files |>
        @filter( any(map(ss->contains(_.target,ss),calibration_target_substrings)) ) |>
        DataFrame

df_files_calib[!,:airmass] = fill(NaN,size(df_files_calib,1))
CSV.write(manifest_calib_filename, df_files_calib)

df_files_obs = df_files |>
    @filter( !any(map(ss->contains(_.target,ss),calibration_target_substrings)) ) |>
    DataFrame

df_files_use = df_files_obs |>
      @filter( _.target == fits_target_str ) |>
      #@filter(bjd_first_good <= _.bjd < bjd_last_good) |>
      #@filter( is_good_day(_.Filename) ) |>
      #@take(max_spectra_to_use) |>
      DataFrame

if fits_target_str == "Sun" || fits_target_str == "Solar"
    df_files_use[!,:alt_sun] = calc_solar_alt.(df_files_use.bjd,obs=:WIYN)
    df_files_use[!,:airmass] = DifferentialExtinction.f_airmass.(DifferentialExtinction.f_z.(deg2rad.(df_files_use.alt_sun)))
    println("# Computing differential extinction")
    try
        df_files_use[!,:Δv_diff_ext] = calc_Δv_diff_extinction.(df_files_use.bjd, obs=:WIYN)
    catch
       try
          println("# Encountered error with call to JplHorizons, sleeping and will try again.")
          sleep(60)
          df_files_use[!,:Δv_diff_ext] = calc_Δv_diff_extinction.(df_files_use.bjd, obs=:WIYN)
       catch
          println("# Encountered error with call to JplHorizons, using NaNs to proceed with rest of pipeline.")
          df_files_use[!,:Δv_diff_ext] = fill(NaN,length(df_files_use.bjd))
       end
    end
    println("# FWHM effect")
    df_files_use[!,:Δfwhm²] = CalculateFWHMDifference_SolarRotation_from_obs.(df_files_use.bjd,obs=:WIYN)
end

if create_missing_continuum_files
   df_files_use.continuum_filename = map(fn->joinpath(output_path,"continuum",match(r"(neidL1_\d+[T_]\d+)\.fits$", fn)[1] * "_continuum=afs.jld2"),df_files_use.Filename)
end

#=
if verbose println("# Reading in customized parameters from param.jl.")  end
   if !@isdefined(idx_day_to_use)
       idx_day_to_use = 1
   end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))
   if match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] ==  match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1]
      date_str = match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]
    else
      date_str = string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]) * "-" * string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1])
   end
   outputs["df_files_use"] = df_files_use

   outputs_filename = joinpath(output_dir,"solar_" * date_str * "_new.jld2")
   if isfile(outputs_filename) && false
     times_already_processed = load(outputs_filename, "times")
     files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1)
      if files_in_day_to_process == length(times_already_processed)
         println("# Already processed all ", length(times_already_processed), " files for ", date_str)
         exit()
      end
   end
=#

#using NaNMath
println("# Computing order SNRs")
df_files_use[!,:order_snrs] = fill(zeros(0),size(df_files_use,1))
   for (i,row) in enumerate(eachrow(df_files_use))
       #if row.target != "Sun" continue   end
       if any(map(ss->contains(row.target,ss),calibration_target_substrings))
           println("# Warning: target = ", row.target)
           continue
       end
       spec = NEID.read_data(row)
       all_orders = min_order(NEID2D()):max_order(NEID2D())
       pixels = min_pixel_in_order(NEID2D()):max_pixel_in_order(NEID2D())
       order_snr = map(ord->RvSpectMLBase.calc_snr(spec, pixels, ord), all_orders)
       df_files_use[i,:order_snrs] = order_snr
   end

println("# Writing manifest file.")
CSV.write(manifest_filename, df_files_use)
      # Can extract vector of Reals turned into a string via:
      # df2 = CSV.read("test.csv", DataFrame)
      # map(i->parse.(Float64,split(df2[i,:order_snrs][2:end-1],',')),1:size(df2,1))

#snr_matrix = hcat(df_files_use[!,:order_snrs]...)

#=
using JLD2, FileIO
if create_missing_continuum_files
   println("# Creating missing continuum files.")
   for row in eachrow(df_files_use)
     if isfile(row.continuum_filename) continue end
     println("# Need to make ", row.continuum_filename )

     output_filename = row.continuum_filename
     tmpdir = basename(output_filename)
     #=
     if !isdir(tmpdir)
        println("# Making ", tmpdir)
        mkdir(tmpdir)
     end
     =#
     #continue
     spec = NEID.read_data(row.Filename)
     try
        continuum = Continuum.calc_continuum_model(spec)
        @save output_filename continuum
     catch
        println("# ERROR computing continuum for ", output_filename)
     end
   end
end

=#


text:
#=
   The following function wraps a call that ascertains the floating point value
   deemed to be most representative of the valuespan an interval signifies.

   If you prefer, you may provide a functional substitute by changing
   the function called inside of both mostRepresentativeValue()s.
=#

function mostRepresentativeValue{S<:Sculpt, Q<:Qualia, C<:Clay}(x::Flex{S,Q,C})
    lo,hi = (Q==CLCL) ? value(x) : value(opened(x))
    domainExtendedLogarithmicMean(lo,hi)
end
function mostRepresentativeValue{C<:Clay}(lo::C, hi::C)
    domainExtendedLogarithmicMean(lo,hi)
end

function domainExtendedLogarithmicMean{R<:Real}(x::R,y::R)
       lo,hi = minmax(x,y)
       if lo > 0.0
           #=
               using the logarithmic mean to condense the spanned interval into a representative floating point value
           =#
           (hi - lo) / (log(hi) - log(lo))
        else # determine the logarithmic mean indirectly
           if signbit(hi) == signbit(lo)
                abs(hi - lo) / (log(abs(hi)) - log(abs(lo)))
           else
              lo = abs(hi) / log(abs(lo))
              hi = abs(lo) / log(abs(hi))
              if abs(lo) <= abs(hi)
                 copysign( abs(lo-hi), lo)
              else
                 copysign( abs(lo-hi), hi)
              end
            end
        end
end


text:
export AbBasis, LinearUnit, Relu, Requ

abstract type AbBasis{F<:AbstractFloat} end
abstract type AbUniDirBasis{F}<:AbBasis{F} end
calc{F<:AbUniDirBasis}(::Type{F}, X, u) = calc(F, X*u)

################################################################# Linear
struct LinearUnit{F} <: AbUniDirBasis{F}
  u::Vector{F}
end
calc{F}(::Type{LinearUnit{F}}, x) = x
derv{F}(::Type{LinearUnit{F}}, x) = ones(x)
calc{F}(b::LinearUnit{F}, X) = X*b.u

################################################################# Signed
abstract type AbSignedUniDirBasis{F}<:AbUniDirBasis{F} end
calc(b::B, X) where B<:AbSignedUniDirBasis = b.sgn * calc(B, X*b.u)
derv(b::B, X) where B<:AbSignedUniDirBasis = b.sgn * derv.(B, X*b.u)

################################################################# Relu
mutable struct Relu{F} <: AbSignedUniDirBasis{F}
  u::Vector{F}
  sgn::F
end
Relu{F}(u::Vector{F}) = Relu(u, one(F))
calc{F}(::Type{Relu{F}}, x) = max.(x, 0.)
derv{F}(::Type{Relu{F}}, x) = (x .> 0.) * 1.

################################################################# ReQU
struct Requ{F} <: AbSignedUniDirBasis{F}
  u::Vector{F}
  sgn::F
end
Requ{F}(u::Vector{F}) = Requ(u, one(F))
calc{F}(::Type{Requ{F}}, x) = max.(x, 0) .^ 2
derv{F}(::Type{Requ{F}}, x) = 2. * (x.>0) .* x

################################################################# Huber
struct Huber{t, F} <: AbSignedUniDirBasis{F}
    u::Vector{F}
    sgn::F
end
Huber{t}(u::Vector{F}, s::F) where {t, F} = Huber{t, F}(u, s)
Huber{t}(u::Vector{F}) where {t, F} = Huber{t, F}(u, one(F))
Huber(u::Vector{F}, s::F) where {F} = Huber{1., F}(u, s)
Huber(u::Vector{F}) where {F} = Huber{1., F}(u, one(F))
calc(::Type{Huber{t, F}}, x) where {F, t} = (x.>0) .* ((x*t .< 1) .* (t*x.*x/2) + (x*t .≥ 1) .* (x-.5/t))
derv(::Type{Huber{t, F}}, x) where {F, t} = (x.>0) .* ((x*t .< 1) .* (t*x)      + (x*t .≥ 1)            )

################################################################# Softplus
struct Softplus{t, F} <: AbSignedUniDirBasis{F}
    u::Vector{F}
    sgn::F
end
Softplus{t}(u::Vector{F}, s::F) where {t, F} = Softplus{t, F}(u, s)
Softplus{t}(u::Vector{F}) where {t, F} = Softplus{t, F}(u, one(F))
Softplus(u::Vector{F}, s::F) where {F} = Softplus{1., F}(u, s)
Softplus(u::Vector{F}) where {F} = Softplus{1., F}(u, one(F))
calc(::Type{Softplus{t, F}}, x) where {F, t} = log1plusexp(t*x)/t
derv(::Type{Softplus{t, F}}, x) where {F, t} = sigmoid(t*x)
###############################################################################


text:
using ChainRulesCore
using ZygoteRules: @adjoint

using ..ITensorAutoHOOT
using ..ITensorAutoHOOT: AbstractNetwork

struct SimpleMPSTensor <: AbstractNetwork
  mps::MPS
end

function Base.:+(A::SimpleMPSTensor, B::SimpleMPSTensor; kwargs...)
  if length(A.mps) == 1 && order(A.mps[1]) == 0
    return SimpleMPSTensor(MPS([ITensor(A[] + B[])]))
  else
    return SimpleMPSTensor(+(A.mps, B.mps; kwargs...))
  end
end

function Base.getindex(t::SimpleMPSTensor)
  @assert length(t.mps) == 1
  @assert order(t.mps[1]) == 0
  return t.mps[1][]
end

function ChainRulesCore.rrule(::typeof(getindex), x::SimpleMPSTensor)
  y = x[]
  function getindex_pullback(ȳ)
    x̄ = SimpleMPSTensor(MPS([ITensor(ȳ)]))
    return (NoTangent(), x̄)
  end
  return y, getindex_pullback
end

"""
initialization methods
"""
function SimpleMPSTensor(tensor::ITensor; cutoff, maxdim)
  mps_out = if size(tensor) == ()
    MPS([tensor])
  else
    MPS(tensor, inds(tensor); cutoff=cutoff, maxdim=maxdim)
  end
  return SimpleMPSTensor(mps_out)
end

@adjoint function SimpleMPSTensor(tensor::ITensor; cutoff, maxdim)
  MPSTensor_pullback(dtensor) = (ITensor(dtensor),)
  return SimpleMPSTensor(tensor; cutoff=cutoff, maxdim=maxdim), MPSTensor_pullback
end

"""
contract
"""
function simple_mps_contract(mps1::MPS, mps2::MPS; cutoff, maxdim)
  ## TODO: modify this function based on https://arxiv.org/pdf/1912.03014.pdf
  tensor = contract(vcat(collect(mps1), collect(mps2))...)
  return if size(tensor) == ()
    MPS([tensor])
  else
    MPS(tensor, inds(tensor); cutoff=cutoff, maxdim=maxdim)
  end
end

function ITensors.contract(t1::SimpleMPSTensor, t2::SimpleMPSTensor; cutoff, maxdim)
  mps_out = simple_mps_contract(t1.mps, t2.mps; cutoff=cutoff, maxdim=maxdim)
  return SimpleMPSTensor(mps_out)
end


text:
module ExportData

export AnsembleData, SingleData, MineData, VanDerWaalsData

cd(dirname(@__FILE__))
include("Q1-MD.jl")
datapath = "../Data/"

using ProgressMeter, JLD, Statistics

function AnsembleData(sn::Integer, rn::Integer, Parameters::Dict{Symbol,Real})
    stepnum = sn
    runnum = rn
    Prog = Progress(stepnum * runnum)

    TotUColl = zeros(stepnum, runnum)
    TotKColl = zeros(stepnum, runnum)
    TotTColl = zeros(stepnum, runnum)
    TotPColl = zeros(stepnum, runnum)
    TotlsColl = zeros(stepnum, runnum)

    for j ∈ 1:runnum
        sys = MDSim.init(; Parameters...)
        for i ∈ 1:stepnum
            MDSim.simulate!(sys)
            lsN = count(<=(sys.l / 2), sys.r[1, :]) / sys.N
            TotlsColl[i, j] = lsN
            TotUColl[i, j] = sys.U
            TotKColl[i, j] = sys.K
            TotTColl[i, j] = sys.T
            TotPColl[i, j] = sys.P
            next!(Prog)
            update!(Prog)
        end
    end
    save(datapath * "DataTD.jld", "TotUColl", TotUColl, "TotKColl", TotKColl,
        "TotTColl", TotTColl, "TotPColl", TotPColl, "TotlsColl", TotlsColl)
end

function SingleData(sn::Integer, Parameters::Dict)
    sys = MDSim.init(; Parameters...)

    stepnum = sn
    Prog = Progress(stepnum)
    PosColl = zeros(stepnum, 2, Parameters[:N])
    VelColl = zeros(stepnum, 2, Parameters[:N])

    UColl = zeros(stepnum)
    KColl = zeros(stepnum)
    TColl = zeros(stepnum)
    PColl = zeros(stepnum)

    for i ∈ 1:stepnum
        MDSim.simulate!(sys)
        PosColl[i, :, :] = sys.r
        VelColl[i, :, :] = sys.v
        UColl[i] = sys.U
        KColl[i] = sys.K
        TColl[i] = sys.T
        PColl[i] = sys.P
        next!(Prog)
        update!(Prog)
    end
    save(datapath * "DataSD.jld",
        "PosColl", PosColl, "VelColl", VelColl,
        "UColl", UColl, "KColl", KColl, "TColl", TColl, "PColl", PColl)
end

function MineData(TotUColl::Matrix{T}, TotKColl::Matrix{T}, TotTColl::Matrix{T}, TotPColl::Matrix{T}, TotlsColl::Matrix{T}) where {T<:AbstractFloat}
    TUC = [[TotUColl[i, :]...] for i ∈ 1:size(TotUColl)[1]]
    for i ∈ TUC
        deleteat!(i, findall(x -> isnan(x), i))
    end
    TKC = [[TotKColl[i, :]...] for i ∈ 1:size(TotKColl)[1]]
    for i ∈ TKC
        deleteat!(i, findall(x -> isnan(x), i))
    end
    TTC = [[TotTColl[i, :]...] for i ∈ 1:size(TotTColl)[1]]
    for i ∈ TTC
        deleteat!(i, findall(x -> isnan(x), i))
    end
    TPC = [[TotPColl[i, :]...] for i ∈ 1:size(TotPColl)[1]]
    for i ∈ TPC
        deleteat!(i, findall(x -> isnan(x), i))
    end
    TLSC = [[TotlsColl[i, :]...] for i ∈ 1:size(TotlsColl)[1]]
    for i ∈ TLSC
        deleteat!(i, findall(x -> isnan(x), i))
    end

    save(datapath * "MeanData.jld", "MEANUColl", mean.(TUC), "MEANKColl", mean.(TKC),
        "MEANTColl", mean.(TTC), "MEANPColl", mean.(TPC), "MEANlsColl", mean.(TLSC))
    save(datapath * "STDData.jld", "STDUColl", std.(TUC), "STDKColl", std.(TKC),
        "STDTColl", std.(TTC), "STDPColl", std.(TPC), "STDlsColl", std.(TLSC))
end

function VanDerWaalsData(Parameters::Dict{Symbol,Real})
    sys = MDSim.init(; Parameters...)

    TColl = zeros(round(Int, 1 / sys.h))
    PColl = zeros(round(Int, 1 / sys.h))
    sys.h *= 10
    for i ∈ 1:round(Int, 20 / sys.h)
        MDSim.update_phase!(sys)
        MDSim.check_boundary!(sys)
    end

    sys.h /= 10

    for i ∈ 1:round(Int, 1 / sys.h)
        MDSim.update_phase!(sys)
        MDSim.check_boundary!(sys)
        MDSim.update_temperature!(sys)
        MDSim.update_pressure!(sys)
        TColl[i] = sys.T
        PColl[i] = sys.P
    end

    return std(TColl), mean(TColl), std(PColl), mean(PColl)
end

end


text:
immutable Laplace <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Laplace(l::Real, s::Real)
        if s > 0.0
            new(float64(l), float64(s))
        else
            error("scale must be positive")
        end
    end
end

Laplace(location::Real) = Laplace(location, 1.0)
Laplace() = Laplace(0.0, 1.0)

const Biexponential = Laplace

function cdf(d::Laplace, q::Real)
    return 0.5 * (1.0 + sign(q - d.location) *
           (1.0 - exp(-abs(q - d.location) / d.scale)))
end

entropy(d::Laplace) = log(2.0 * e * d.scale)

insupport(d::Laplace, x::Number) = isreal(x) && isfinite(x)

kurtosis(d::Laplace) = 3.0

mean(d::Laplace) = d.location

median(d::Laplace) = d.location

function mgf(d::Laplace, t::Real)
    m, b = d.location, d.scale
    return exp(t * m) / (1.0 - b^2 * t^2)
end

function cf(d::Laplace, t::Real)
    m, b = d.location, d.scale
    return exp(im * t * m) / (1.0 + b^2 * t^2)
end

modes(d::Laplace) = [d.location]

function pdf(d::Laplace, x::Real)
    (1.0 / (2.0 * d.scale)) * exp(-abs(x - d.location) / d.scale)
end

function logpdf(d::Laplace, x::Real)
    -log(2.0 * d.scale) - abs(x - d.location) / d.scale
end

function quantile(d::Laplace, p::Real)
    d.location - d.scale * sign(p - 0.5) * log(1.0 - 2.0 * abs(p - 0.5))
end

# Need to see whether other RNG strategies are more efficient:
# (1) Difference of two Exponential(1/b) variables
# (2) Ratio of logarithm of two Uniform(0.0, 1.0) variables
function rand(d::Laplace)
    u = rand() - 0.5
    return d.location - d.scale * sign(u) * log(1.0 - 2.0 * abs(u))
end

skewness(d::Laplace) = 0.0

std(d::Laplace) = sqrt(2.0) * d.scale

var(d::Laplace) = 2.0 * d.scale^2


text:
function get_convergence_checker(obj::OptimizationProblem)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluation_counter(obj::OptimizationProblem)
    return jcall(obj, "getEvaluationCounter", Incrementor, ())
end

function get_iteration_counter(obj::OptimizationProblem)
    return jcall(obj, "getIterationCounter", Incrementor, ())
end


text:
########################################
## File Name: belief_transition_models.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Belief Transition Models for SACBP
########################################

import Plots
using LinearAlgebra

abstract type BeliefTransModel <: TransModel end
abstract type GaussianFilter <: BeliefTransModel end
# abstract type DiscreteBayesianFilter <: BeliefTransModel end

function covObservePos(p::Vector{<:Real},qi::Vector{<:Real}) # State-dependent observation noise model used in multi-target localization problem.
    R_scale = 0.01*Matrix(1.0I, 2, 2);
    R_nominal = 0.001*Matrix(1.0I, 2, 2);
    dist = norm(p-qi);
    return R_nominal + dist*R_scale;
end

struct UKFPosRange <: GaussianFilter
# Discrete-time EKF model for belief transition in position estimation with range-only measurements.
end
function trans(model::UKFPosRange,b::BelMvNormal{T},p_robot::PhysPos{U},dt::Real,y::Real,
               Q::Matrix{<:Real}) where T <: Real where U <: Real
    d = Distributions.MvNormal(b);
    μ,Σ = d.μ, Matrix(d.Σ);

    # Predict.
    μ = μ;
    Σ = Σ + Q*dt;

    # Generate Sigma Points.
    n = 2 + 2;
    λ = 2.0;
    x0 = [μ;zeros(2)];
    x_sigma = [];
    push!(x_sigma,x0);
    sqrtm_P = Matrix(cholesky([Σ          zeros(2,2);
                               zeros(2,2) covObservePos(p_robot.pos,μ)]).L);
    for ii = 1:n
        push!(x_sigma,x0 + sqrt(n+λ)*sqrtm_P[:,ii])
    end
    for ii = 1:n
        push!(x_sigma,x0 - sqrt(n+λ)*sqrtm_P[:,ii])
    end
    w_sigma = 1/(2*(n+λ))*ones(2*n+1);
    w_sigma[1] = λ/(n+λ);
    y_sigma = Vector{Real}(undef, 2*n+1)
    for ii = 1:2*n+1
        y_sigma[ii] = norm(p_robot.pos - x_sigma[ii][1:2] + x_sigma[ii][3:4]);
    end
    y_pred_mean = dot(w_sigma,y_sigma);

    # Update
    σ_YY = dot(w_sigma,(y_sigma .- y_pred_mean).^2);
    Σ_XY = w_sigma[1]*(x_sigma[1][1:2] - μ)*(y_sigma[1] - y_pred_mean);
    for ii = 1:2*n
        Σ_XY += w_sigma[ii+1]*(x_sigma[ii+1][1:2] - μ)*(y_sigma[ii+1] - y_pred_mean);
    end
    K = 1/σ_YY*Σ_XY;
    μ = μ + K*(y - y_pred_mean);
    Σ = Σ - K*Σ_XY';
    Σ = (Σ + Σ')/2;
    return BelMvNormal(p_robot.t,μ,Σ)
end
function trans_jacobi_auto_p(model::UKFPosRange,b::BelMvNormal{T},p_robot::PhysPos{U},
                             dt::Real,y::Real,Q::Matrix{<:Real}) where T <: Real where U <: Real
    function f(pVec)
        # if isa(p_robot,PhysPos)
            return trans(model,b,PhysPos(p_robot.t,pVec),dt,y,Q).params
        # else
        #     return trans(model,b,PhysPosHead2D(p_robot.t,pVec),dt,y,Q).params
        # end
    end
    cfg = ForwardDiff.JacobianConfig(f,vec(p_robot),ForwardDiff.Chunk{2}());
    return ForwardDiff.jacobian(f,vec(p_robot),cfg)
    #return ForwardDiff.jacobian(f,vec(p_robot))
end
function trans_jacobi_auto_b(model::UKFPosRange,b::BelMvNormal{T},p_robot::PhysPos{U},
                             dt::Real,y::Real,Q::Matrix{<:Real}) where T <: Real where U <: Real
    function f(bVec)
        b1 = BelMvNormal(b.t,bVec)
        return trans(model,b1,p_robot,dt,y,Q).params
    end
    cfg = ForwardDiff.JacobianConfig(f,b.params,ForwardDiff.Chunk{6}());
    return ForwardDiff.jacobian(f,b.params,cfg)
    #return ForwardDiff.jacobian(f,b.params)
end
function trans(model::UKFPosRange,s::AugState{PhysPos{U},BelMvNormal{T}},dt::Real,y::Real,
               Q::Matrix{<:Real})::AugState{PhysPos{U},BelMvNormal{T}} where T <: Real where U <: Real
    b_new = trans(model,s.b,s.p,dt,y,Q);
    return AugState(s.p,b_new)
end
function trans_jacobi_auto_p(model::UKFPosRange,s::AugState{PhysPos{U},BelMvNormal{T}},
                             dt::Real,y::Real,Q::Matrix{<:Real}) where T <: Real where U <: Real
    return trans_jacobi_auto_p(model,s.b,s.p,dt,y,Q);
end
function trans_jacobi_auto_b(model::UKFPosRange,s::AugState{PhysPos{U},BelMvNormal{T}},
                             dt::Real,y::Real,Q::Matrix{<:Real}) where T <: Real where U <: Real
    return trans_jacobi_auto_b(model,s.b,s.p,dt,y,Q);
end
function trans(model::UKFPosRange,bvec::VecBelMvNormal{T},p_robot::PhysPos{U},dt::Real,yvec::Vector{<:Real},
               Q::Matrix{<:Real})::VecBelMvNormal{T} where T <: Real where U <: Real
    if length(bvec) != length(yvec)
        error(ArgumentError("Invalid number of observations."))
    else
        bvecNew = similar(bvec);
        for ii = 1:length(bvec)
            bvecNew[ii] = trans(model,bvec[ii],p_robot,dt,yvec[ii],Q)
        end
        return bvecNew
    end
end
function trans_jacobi_auto_p(model::UKFPosRange,bvec::VecBelMvNormal{T},
                             p_robot::PhysPos{U},dt::Real,yvec::Vector{<:Real},Q::Matrix{<:Real}) where T <: Real where U <: Real
    JArray = Vector{Matrix{T}}(undef, length(bvec));
    for ii = 1:length(bvec)
        JArray[ii] = trans_jacobi_auto_p(model,bvec[ii],p_robot,dt,yvec[ii],Q);
    end
    return JArray
end
function trans_jacobi_auto_b(model::UKFPosRange,bvec::VecBelMvNormal{T},
                             p_robot::PhysPos{U},dt::Real,yvec::Vector{<:Real},Q::Matrix{<:Real}) where T <: Real where U <: Real
    JArray = Vector{Matrix{T}}(undef, length(bvec));
    for ii = 1:length(bvec)
        JArray[ii] = trans_jacobi_auto_b(model,bvec[ii],p_robot,dt,yvec[ii],Q)
    end
    return JArray
end
function trans(model::UKFPosRange,s::AugState{PhysPos{U},VecBelMvNormal{T}},dt::Real,yvec::Vector{<:Real},
               Q::Matrix{<:Real})::AugState{PhysPos{U},VecBelMvNormal{T}} where T <: Real where U <: Real
    if length(s.b) != length(yvec)
        error(ArgumentError("Invalid number of observations."))
    else
        bvecNew = similar(s.b);
        for ii = 1:length(s.b)
            bvecNew[ii] = trans(model,s.b[ii],s.p,dt,yvec[ii],Q);
        end
        return AugState(s.p,bvecNew)
    end
end
function trans_jacobi_auto_p(model::UKFPosRange,s::AugState{PhysPos{U},VecBelMvNormal{T}},
                             dt::Real,yvec::Vector{<:Real},Q::Matrix{<:Real}) where T <: Real where U <: Real
    JArray = Vector{Matrix{T}}(undef, length(s.b));
    for ii = 1:length(s.b)
        JArray[ii] = trans_jacobi_auto_p(model,s.b[ii],s.p,dt,yvec[ii],Q);
    end
    return JArray
end
function trans_jacobi_auto_b(model::UKFPosRange,s::AugState{PhysPos{U},VecBelMvNormal{T}},
                             dt::Real,yvec::Vector{<:Real},Q::Matrix{<:Real}) where T <: Real where U <: Real
    JArray = Vector{Matrix{T}}(undef, length(s.b));
    for ii = 1:length(s.b)
        JArray[ii] = trans_jacobi_auto_b(model,s.b[ii],s.p,dt,yvec[ii],Q);
    end
    return JArray
end
function trans(model::UKFPosRange,b::BelMvNormal{T},p_robot::PhysPos{U},dt::Real,
               Q::Matrix{<:Real},rng::MersenneTwister)::Tuple{Real,BelMvNormal{T}} where T <: Real where U <: Real
    d = Distributions.MvNormal(b);
    μ,Σ = d.μ,Matrix(d.Σ);

    # Predict.
    μ = μ;
    Σ = Σ + Q*dt;

    # Generate Sigma Points.
    n = 2 + 2;
    λ = 2.0;
    x0 = [μ;zeros(2)];
    x_sigma = [];
    push!(x_sigma,x0);
    sqrtm_P = Matrix(cholesky([Σ          zeros(2,2);
                               zeros(2,2) covObservePos(p_robot.pos,μ)]).L);
    for ii = 1:n
        push!(x_sigma,x0 + sqrt(n+λ)*sqrtm_P[:,ii])
    end
    for ii = 1:n
        push!(x_sigma,x0 - sqrt(n+λ)*sqrtm_P[:,ii])
    end
    w_sigma = 1/(2*(n+λ))*ones(2*n+1);
    w_sigma[1] = λ/(n+λ);
    y_sigma = Vector{Real}(undef, 2*n+1)
    for ii = 1:2*n+1
        y_sigma[ii] = norm(p_robot.pos - x_sigma[ii][1:2] + x_sigma[ii][3:4]);
    end
    y_pred_mean = dot(w_sigma,y_sigma);

    # Update
    σ_YY = dot(w_sigma,(y_sigma .- y_pred_mean).^2);
    Σ_XY = w_sigma[1]*(x_sigma[1][1:2] - μ)*(y_sigma[1] - y_pred_mean);
    for ii = 1:2*n
        Σ_XY += w_sigma[ii+1]*(x_sigma[ii+1][1:2] - μ)*(y_sigma[ii+1] - y_pred_mean);
    end
    K = 1/σ_YY*Σ_XY;
    observeModel = ObserveModel_Range();
    y = observe(observeModel,p_robot,μ,Σ+covObservePos(p_robot.pos,μ),rng);
    μ = μ + K*(y - y_pred_mean);
    Σ = Σ - K*Σ_XY';
    Σ = (Σ + Σ')/2;
    return y,BelMvNormal(p_robot.t,μ,Σ)
end
function trans(model::UKFPosRange,s::AugState{PhysPos{U},BelMvNormal{T}},dt::Real,
               Q::Matrix{<:Real},rng::MersenneTwister)::Tuple{Real,AugState{PhysPos{U},BelMvNormal{T}}} where T <: Real where U <: Real
    y,b_new = trans(model,s.b,s.p,dt,Q,rng);
    return y,AugState(s.p,b_new)
end
function trans(model::UKFPosRange,bvec::VecBelMvNormal{T},p_robot::PhysPos{U},dt::Real,
               Q::Matrix{<:Real},rng::MersenneTwister)::Tuple{Vector{<:Real},VecBelMvNormal{T}} where T <: Real where U <: Real
    temp = map(b -> trans(model,b,p_robot,dt,Q,rng),bvec)
    yvec = Vector{Float64}(undef, length(bvec));
    bvecNew = similar(bvec);
    for ii = 1:length(bvec)
        yvec[ii] = temp[ii][1];
        bvecNew[ii] = temp[ii][2];
    end
    return yvec,bvecNew
end
function trans(model::UKFPosRange,s::AugState{PhysPos{U},VecBelMvNormal{T}},dt::Real,
               Q::Matrix{<:Real},rng::MersenneTwister)::Tuple{Vector{<:Real},AugState{PhysPos{U},VecBelMvNormal{T}}} where T <: Real where U <: Real
    temp = map(b -> trans(model,b,s.p,dt,Q,rng),s.b)
    yvec = Vector{Float64}(undef, length(s.b))
    bvecNew = similar(s.b);
    for ii = 1:length(s.b)
        yvec[ii] = temp[ii][1];
        bvecNew[ii] = temp[ii][2];
    end
    sNew = AugState(s.p,bvecNew);
    return yvec,sNew
end

struct CPredictManipulate2D <: GaussianFilter
# Linearized Continuous Gaussian Prediction model for belief transition in planar manipulation problem.
end
function trans(model::CPredictManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},dt::Real,
               Q::Matrix{Float64},mode=1)::BelMvNormal{T} where T <: Real
# Continuous Time Transition. Euler Approximation.
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."))
    else
        if mode!=1  # mode n != 1 is with naive euler approximation (for automatic jacobian computation with dt == 1.)
            transModel = TransModel_Manipulate2D();
            d = Distributions.MvNormal(b);
            μ,Σ = d.μ,Matrix(d.Σ);
            μMP = PhysManipulate2D(b.t,μ)
            A = trans_jacobi(transModel,μMP,u);
            Σ = Σ + (A*Σ + Σ*A' + Q)*dt; #This Euler Approximation is numerically unstable for large dt (e.g. dt == 0.01). Using Square-Root form.
            μ = vec(trans(transModel,PhysManipulate2D(b.t,μ),u,dt));
            t = round(b.t + dt,digits=5);
            return BelMvNormal(t,μ,Σ)
        else #mode 1 is with square-root form update.
            #println("In trans(): t = $(b.t), $(b.params[b.dim+1:end])")
            transModel = TransModel_Manipulate2D();
            try d = Distributions.MvNormal(b);
            catch
                # This block is solely for debugging and needs to be deleted before release.
                println("In trans(): t = $(b.t), μ = $(b.params[1:b.dim]),Σ = $(b.params[b.dim+1:end])")
                d = Distributions.MvNormal(b.params[1:b.dim],reshape(b.params[b.dim+1:end],11,11)+1e-5*Matrix(1.0I, 11, 11))
            end
            #d = Distributions.MvNormal(b);
            μ,Σ = d.μ,Matrix(d.Σ);
            μMP = PhysManipulate2D(b.t,μ)
            A = trans_jacobi(transModel,μMP,u);
            L = cholesky(Σ).L;
            L_inv = inv(L);
            L = L + (A*L + 1/2*Q*L_inv)*dt;
            Σ = Matrix(L*L');
            #A_discrete_approx = expm(A*dt);
            #Σ = A_discrete_approx*Σ*A_discrete_approx' + Q*dt; # Use Matrix Exponential and treat like a discrete-time system.
            #Σ = (Σ + Σ')/2;
            μ = vec(trans(transModel,PhysManipulate2D(b.t,μ),u,dt));
            #if μ[7] <= 0. || μ[8] <= 0. || μ[11] <= 0.
            #    println(μ)
            #end
            t = round(b.t + dt,digits=5);
            return BelMvNormal(t,μ,Σ)
        end
    end
end
function trans_jacobi(model::CPredictManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},Q::Matrix{Float64}) where T <: Real
    function f(μ)
        d = Distributions.MvNormal(b);
        Σ = Matrix(d.Σ);
        A = trans_jacobi(TransModel_Manipulate2D(),PhysManipulate2D(b.t,μ),u);
        return vcat(A*Σ + Σ*A' + Q...);
    end
    Jacobian = zeros(T,11^2+11,11^2+11);
    A = trans_jacobi(TransModel_Manipulate2D(),PhysManipulate2D(b.t,b.params[1:11]),u);
    Jacobian[1:11,1:11] = A;
    for ii = 1:11
        Jacobian[12+11*(ii-1):11+11*ii,12+11*(ii-1):11+11*ii] += A;
        for jj = 1:11
            Jacobian[12+11*(ii-1):11+11*ii,12+11*(jj-1):11+11*jj] += diagm(A[ii,jj]*ones(11));
        end
    end
    cfg = ForwardDiff.JacobianConfig(f,b.params[1:11],ForwardDiff.Chunk{11}())
    Jacobian[12:end,1:11] = ForwardDiff.jacobian(f,b.params[1:11])#,cfg);
    return Jacobian
end
function trans_jacobi_auto(model::CPredictManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},Q::Matrix{Float64}) where T <: Real
    function f(params)
        bel = BelMvNormal(b.t,params);
        return trans(model,bel,u,1.,Q,0).params - params;
    end
    cfg = ForwardDiff.JacobianConfig(f,b.params,ForwardDiff.Chunk{40}());
    return ForwardDiff.jacobian(f,b.params,cfg)
    #return ForwardDiff.jacobian(f,b.params)
end
function trans_u_coeff(model::CPredictManipulate2D,b::BelMvNormal{T},Q::Matrix{Float64}) where T <: Real
    uVecs = [[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]];
    u0 = MControl2D(b.t,0.,0.,0.);
    H = Matrix{T}(undef,length(b.params),3);
    transModel = CPredictManipulate2D();
    for ii = 1:3
        u = MControl2D(b.t,uVecs[ii]);
        H[:,ii] = trans(transModel,b,u,1.,Q,0).params - trans(transModel,b,u0,1.,Q,0).params
    end
    return H
end
struct DUpdateManipulate2D <: GaussianFilter
# Discrete-time EKF Update model for belief transition in planar manipulation problem.
end
function trans(model::DUpdateManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},y::Vector{<:Real},
               R::Matrix{Float64})::BelMvNormal{T} where T <: Real
# Instantaneous Transition (EKF Update Step)
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."))
    else
        observeModel = ObserveModel_Manipulate2D();
        d = Distributions.MvNormal(b);
        μ,Σ = d.μ,Matrix(d.Σ);
        μMP = PhysManipulate2D(b.t,μ);
        C = observe_jacobi(observeModel,μMP,u);
        K = Σ*C'/(C*Σ*C' + R);
        y_pred_mean = observe(observeModel,μMP,u);
        if y[3] - y_pred_mean[3] > pi   # Correct angular difference.
            y[3] -= 2*pi;
        elseif y[3] - y_pred_mean[3] < -pi
            y[3] += 2*pi;
        end
        μ += K*(y - y_pred_mean);
        μ[3] = mod2pi(μ[3]);
        Σ -= K*C*Σ;
        Σ = (Σ + Σ')/2;
        t = b.t;
        return BelMvNormal(t,μ,Σ)
    end
end
function trans(model::DUpdateManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},
               R::Matrix{Float64},rng::MersenneTwister)::Tuple{Vector{<:Real},BelMvNormal{T}} where T <: Real
    # Instantaneous Transition with sampled observation. (EKF Update Step)
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."))
    else
        observeModel = ObserveModel_Manipulate2D();
        d = Distributions.MvNormal(b);
        μ,Σ = d.μ,Matrix(d.Σ);
        μMP = PhysManipulate2D(b.t,μ);
        C = observe_jacobi(observeModel,μMP,u);
        H = C*Σ*C' + R;
        H = (H+H')/2;
        K = Σ*C'/H;
        y = observe(observeModel,μMP,u,H,rng);
        y_pred_mean = observe(observeModel,μMP,u);
        if y[3] - y_pred_mean[3] > pi   # Correct angular difference.
            y[3] -= 2*pi;
        elseif y[3] - y_pred_mean[3] < -pi
            y[3] += 2*pi;
        end
        μ += K*(y - y_pred_mean);
        μ[3] = mod2pi(μ[3]);
        U = cholesky((H+H')/2).L; # Square-root form (Andrews, 1968)
        V = cholesky(R).L;
        L = cholesky(Σ).L;
        L -= Σ*C'/((U+V)*U')*C*L;
        Σ = Matrix(L*L');
        #Σ -= K*C*Σ;
        #Σ = (Σ + Σ')/2;
        #try MvNormal(zeros(11),Σ)
        #catch
        #    Σ += eye(11)*2*minimum(eig(Σ)[1]);
        #end
        t = b.t;
        return y,BelMvNormal(t,μ,Σ)
    end
end
function trans_jacobi_auto(model::DUpdateManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},y::Vector{<:Real},R::Matrix{Float64}) where T <: Real
    function f(params)
        bel = BelMvNormal(b.t,params);
        return trans(model,bel,u,y,R).params
    end
    cfg = ForwardDiff.JacobianConfig(f,b.params,ForwardDiff.Chunk{40}());
    #return ForwardDiff.jacobian(f,b.params,cfg)
    return ForwardDiff.jacobian(f,b.params)
end

struct CDEKFManipulate2D <: GaussianFilter
# Continuous-Discrete EKF model for belief transition in planar manipulation problem.
end
function trans(model::CDEKFManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},dtc::Real,dto::Real,y::Vector{<:Real},
               Q::Matrix{Float64},R::Matrix{Float64})::BelMvNormal{T} where T <: Real
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."))
    else
        predictModel = CPredictManipulate2D();
        for ii = 1:Int64(dto/dtc)
            b = trans(predictModel,b,MControl2D(b.t,u.fx,u.fy,u.tr),dtc,Q);
        end
        updateModel = DUpdateManipulate2D();
        b = trans(updateModel,b,MControl2D(b.t,u.fx,u.fy,u.tr),y,R);
        return b
    end
end
function trans(model::CDEKFManipulate2D,b::BelMvNormal{T},u::MControl2D{<:Real},dtc::Real,dto::Real,
               Q::Matrix{Float64},R::Matrix{Float64},rng::MersenneTwister)::Tuple{Vector{<:Real},BelMvNormal{T}} where T <: Real
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."))
    else
        predictModel = CPredictManipulate2D();
        for ii = 1:Int64(dto/dtc)
            b = trans(predictModel,b,MControl2D(b.t,u.fx,u.fy,u.tr),dtc,Q);
        end
        updateModel = DUpdateManipulate2D();
        y,b = trans(updateModel,b,MControl2D(b.t,u.fx,u.fy,u.tr),R,rng);
        return y,b
    end
end