text:
"""
# DocGenerator

This is an "experimental" generator for the Markdown documentation.

"experimental" because some (a lot) of undocumented Julia methods are used.

## Why a generator?

Because I want to support (or I want to have)

* (automatically generated) cross references (without an additional syntax)
* inclusion of functions doc-string into the documentation (with an easy
  syntax)
* generation of an index

## How?

There are 3 main ingredients

* where something is documented (in the end), that's a `DocRef`
* what is documented, that's a `ObjSpecification`
* the relationship of "where-is-what" (i.e. what objects are documented where)
  thats `refs` inside `DocCreationEnvironment`

"""
module DocGenerator

using Base.Docs
using Printf
using Markdown
using Cairo
using ComplexVisual
@ComplexVisual.import_huge()

"""
all "registerted" doc-Modules producing some documentation file(s).

Registration at the bottom of this file.
"""
const doc_modules = Vector{Module}()

"""
Source of Documentation.
"""
struct DocSource
    filename   :: AbstractString  # without path and without suffix
    doc_module :: Module          # Module which created the documentation
end

"""
Reference to a point/anchor in the (generated) documentation.
"""
struct DocRef
    source   :: DocSource
    anchor   :: AbstractString
end

"""
Specifies an object (which will be documented). `sig` is a signature
in order to identify a special method if `obj_sym` is the symbol for
a function.
"""
struct ObjSpecification
    mod      :: Module
    obj_sym  :: Symbol
    sig      :: Type
end

"""
Regex for how objects can be specified: `[Module.]Name[:Signature]`.
"""
const re_parse_obj_specification = r"""
    ^\s*                      # maybe space(s)
    (                         # group 1: optional: the module
      ([^\s]+)                # group 2: module name
      \.                      # followed by dot
    )?
    ([^\s]+)                  # group 3: object/symbol name
    (                         # group 4: optional: the signature
      :\s*                    # colon (with maybe space(s))
      (                       # group 5: signature
        Union(.*)  |          # Starting with Union or
        Tuple(.*)             # Starting with Tuple
      )
    )?
    \s*$                      # maybe space(s) to the end
    """x

function parse_obj_specification(str::AbstractString;
        default_module::Module=ComplexVisual,
        default_sig::Type=Union{})
    match_obj = match(re_parse_obj_specification, str)
    if match_obj === nothing
        error("Cannot parse object specification: " * str)
    end
    return ObjSpecification(
        match_obj[2] === nothing ?
            default_module :
            getproperty(@__MODULE__, Symbol(match_obj[2]))::Module,
        Symbol(match_obj[3]),
        match_obj[5] === nothing ?
            default_sig :
            eval(Meta.parse(match_obj[5]))::Type)
end

"""
saves data for generating documents.

Gathers all the `DocRef` for object specifications.
"""
struct DocCreationEnvironment
    refs :: Dict{Module, Dict{Symbol, Dict{Type, Vector{DocRef}}}}
end
DocCreationEnvironment() = DocCreationEnvironment(Dict{Module,
    Dict{Symbol, Dict{Type, Vector{DocRef}}}}())

"""
context while generating the documentation.
"""
struct DocContext
    env      :: DocCreationEnvironment
    source   :: DocSource               # current source
end

"""
A produced (Markdown) document.
"""
struct Document
    source     :: DocSource
    markdown   :: Markdown.MD
end

# Helper functions {{{

eval_in_context(context::DocContext) = context.source.doc_module.eval

"""
```
extract_lines(filename, start_re, stop_re; include_start_stop=true)

filename             AbstractString
start_re             Regexp
stop_re              Regexp
include_start_stop   Bool
```

extract some lines of a file.
"""
function extract_lines(filename::AbstractString, start_re::Regex, stop_re::Regex;
        include_start_stop=true) # {{{
    inside_flag = false
    found_lines = Vector{AbstractString}()
    open(filename, "r") do fio
        while !eof(fio)
            line = readline(fio)

            if match(start_re, line) !== nothing
                if include_start_stop
                    push!(found_lines, line)
                end
                inside_flag = true
            elseif inside_flag && match(stop_re, line) !== nothing
                if include_start_stop
                    push!(found_lines, line)
                end
                inside_flag = false
            elseif inside_flag
                push!(found_lines, line)
            end
        end
    end
    return found_lines
end # }}}

"""
```
extract_lines(doc_context, start_re, stop_re; include_start_stop=true)

doc_context          DocContext
start_re             Regexp
stop_re              Regexp
include_start_stop   Bool
```

extract some lines of a file.
"""

function extract_lines(doc_context::DocContext, start_re::Regex, stop_re::Regex;
        include_start_stop=true)
    filename = doc_context.source.filename * ".jl"
    return extract_lines(filename, start_re, stop_re; include_start_stop)
end

"""
append to Markdown other Markdown-elements (and strip off `MD`s).
"""
function append_md(where_to_append, what_to_append)
    while what_to_append isa Markdown.MD  &&
            length(what_to_append.content) == 1
        what_to_append = what_to_append.content[1]
    end
    if what_to_append isa Markdown.MD
        push!(where_to_append.content, what_to_append.content...)
    else
        push!(where_to_append.content, what_to_append)
    end
    return nothing
end

# }}}

# Methods for refs in DocCreationEnvironment {{{
"""
returns `Vector{DocRef}` for the given object specification or `nothing`
if no `DocRef` can be found.

if `sig === UnionAll` then all the `DocRef` for every saved signature
are put in the result vector.
"""
function get_doc_refs(env::DocCreationEnvironment,
        mod::Module, sym::Symbol, sig::Type)  # {{{
    refs = env.refs
    sym_dict = get(refs, mod, nothing)
    if sym_dict === nothing
        return nothing
    end
    sig_dict = get(sym_dict, sym, nothing)
    if sig_dict === nothing
        return nothing
    end
    if sig === UnionAll
        docref_vec = Vector{DocRef}()
        append!(docref_vec, values(sig_dict)...)
        return docref_vec
    end
    docref_vec = get(sig_dict, sig, nothing)
    return docref_vec
end  # }}}

function get_doc_refs(env::DocCreationEnvironment, obj::ObjSpecification)
    return get_doc_refs(env, obj.mod, obj.obj_sym, obj.sig)
end

"""
appends the `doc_ref` to the refs of the specified object specification.
"""
function append_doc_ref(env::DocCreationEnvironment,
        mod::Module, sym::Symbol, sig::Type, doc_ref::DocRef)   # {{{
    refs = env.refs
    if !haskey(refs, mod)
        refs[mod] = Dict{Symbol, Dict{Type, Vector{DocRef} } }()
    end
    sym_dict = refs[mod]
    if !haskey(sym_dict, sym)
        sym_dict[sym] = Dict{Type, Vector{DocRef}}()
    end
    sig_dict = sym_dict[sym]
    if !haskey(sig_dict, sig)
        sig_dict[sig] = Vector{DocRef}()
    end
    push!(sig_dict[sig], doc_ref)
    return nothing
end # }}}
function append_doc_ref(env::DocCreationEnvironment,
        obj::ObjSpecification, doc_ref::DocRef)
    return append_doc_ref(env, obj.mod, obj.obj_sym, obj.sig, doc_ref)
end
# }}}

# Anchors {{{
"""
test if a `DocRef` is already saved in the `refs` in the DocCreationEnvironment.
"""
function has_docref(context::DocContext, search_ref::DocRef) # {{{
    refs = context.env.refs
    for doc_module in keys(refs)
        for sym in keys(refs[doc_module])
            for docref_vec in values(refs[doc_module][sym])
                if findfirst(
                    doc_ref -> doc_ref == search_ref, docref_vec) !== nothing
                    return true
                end
            end
        end
    end
    return false
end # }}}

"""
try to guess the anchor name which github's flavored markdown will use if
this content is used (again) in this context (as anchor).

Takes into account if there were already anchors for the same content and
use the "-1", "-2", ... suffixes used.
"""
function get_anchor_name(context::DocContext,
        content::AbstractString, prefix::AbstractString="user-content-") # {{{
    name = lowercase(content)                #  ⎫ only a rough approximation
    name = replace(name, r"[^\w\- ]" => "")  #  ⎬ what github/gitlab are
    name = replace(name, " " => "-")         #  ⎭ doing

    refs = context.env.refs
    result = prefix * name
    if result != ""
        base, next_index = result, 1
        while true
            try_ref = DocRef(context.source, result)
            !has_docref(context, try_ref) && break
            result = base * "-" * string(next_index)
            next_index += 1
        end
    end
    return result
end # }}}

doc_ref_to_linktarget(ref::DocRef) = @sprintf(
    "./%s.md#%s", ref.source.filename, ref.anchor)

doc_ref_to_markdownlink(ref::DocRef, name) =
    Markdown.Link(Markdown.Code(String(name)), doc_ref_to_linktarget(ref))
# }}}

# Functions to extract Markdown subparts {{{

const MD_has_content = Union{Markdown.Paragraph, Markdown.BlockQuote,
        Markdown.Admonition, Markdown.MD}
const MD_has_items = Union{Markdown.List}
const MD_has_rows = Union{Markdown.Table}
const MD_has_subparts = Union{MD_has_content, MD_has_items, MD_has_rows}

get_subparts(md) = []
get_subparts(md::MD_has_content) = md.content
get_subparts(md::MD_has_items) = md.items
get_subparts(md::MD_has_rows) = md.rows
# }}}

# Subst `{func: name_of_function}`  {{{
const re_substitute_code_func = r"^\s*\{func:\s*([^\s}]+)\s*\}"

"""
look for `Markdown.Code` with content "{func: name_of_function}"
and replace the content with the acutal code of the function.
"""
substitute_code_func(context::DocContext, md) = nothing
function substitute_code_func(context::DocContext, v::Vector)
    for part in v
        substitute_code_func(context, part)
    end
    return nothing
end
substitute_code_func(context::DocContext, md::MD_has_subparts) =
    substitute_code_func(context, get_subparts(md))
function substitute_code_func(context::DocContext, md::Markdown.Code)
    match_obj = match(re_substitute_code_func, md.code)
    if match_obj !== nothing
        start = Regex("^function "*match_obj.captures[1])
        stop = r"^end"
        func_lines = extract_lines(context, start, stop)
        md.code = join(func_lines, "\n")
    end
    return nothing
end
# }}}

# Subst `{image_from_canvas: expression}` {{{
const re_substitute_canvas_image = r"^\s*\{image_from_canvas:\s*([^\s}]+)\s*\}"

"""
look for `Markdown.Image` with `url` of the form

```
  {image_from_canvas: expression}
```

Then

* execute expression (expect a canvas)
* create a png (with the content of the canvas) and save it unter the
  alt-tag-name
* replace the url with the name of the created png.
"""
substitute_canvas_image(context::DocContext, md) = nothing
function substitute_canvas_image(context::DocContext, v::Vector)
    for part in v
        substitute_canvas_image(context, part)
    end
    return nothing
end
substitute_canvas_image(context::DocContext, md::MD_has_subparts) =
    substitute_canvas_image(context, get_subparts(md))
substitute_canvas_image(context::DocContext, md::Union{
    Markdown.Header, Markdown.Link}) =
    substitute_canvas_image(context, md.text)

function substitute_canvas_image(context::DocContext, md::Markdown.Image)
    match_obj = match(re_substitute_canvas_image, md.url)
    if match_obj !== nothing
        to_eval = Meta.parse(match_obj.captures[1])
        img_file = md.alt
        println("  producing image ", img_file, " for ", to_eval)
        canvas = eval_in_context(context)(to_eval) :: CV_2DCanvas
        write_to_png(canvas.surface, img_file)
        md.url = img_file
    end
    return nothing
end
# }}}

# Subst Paragraphs of the form `[inline](<name>)` {{{
"""
look for Paragraphs which have exactly one link of the form
`[inline](<object specification>)`.

Then replace the whole paragraph with the Markdown documentation of the
specified object.
"""
get_inline_replacement_for(context::DocContext, md) = nothing
function get_inline_replacement_for(context::DocContext,
        md::Markdown.Paragraph)
    if length(md.content) == 1 && md.content[1] isa Markdown.Link
        link = md.content[1]
        if link.text isa Vector && length(link.text) == 1 &&
                link.text[1] == "inline"
            obj_spec = parse_obj_specification(link.url; default_sig=UnionAll)

            if obj_spec.sig !== UnionAll
                # special signature is given
                binding = Base.Docs.Binding(obj_spec.mod, obj_spec.obj_sym)
                doc_md = doc(binding, obj_spec.sig)
            else
                # all docs that can be found for this symbol
                doc_md = doc(getproperty(obj_spec.mod, obj_spec.obj_sym))
            end
            if doc_md === nothing
                error("Cannot find doc for " * link.url)
            end
            return doc_md.content
        end
    end
    return nothing
end
substitute_inline_markdown(context::DocContext, md) = nothing
function substitute_inline_markdown(context::DocContext, v::Vector)
    index = 1
    while index <= length(v)
        part = v[index]
        replacement = get_inline_replacement_for(context, part)
        if replacement !== nothing
            deleteat!(v, index)
            for elem in reverse(replacement)
                insert!(v, index, elem)
            end
            continue
        end
        substitute_inline_markdown(context, part)
        index += 1
    end
    return nothing
end
substitute_inline_markdown(context::DocContext, md::MD_has_subparts) =
    substitute_inline_markdown(context, get_subparts(md))
# }}}

# Extend Headers of the form `doc: <obj specification>` with object documentation {{{
"""
rewrite header content (by appending the code-points in "U+..."-notation)
if needed. This is needed if no anchor-name could be deduced.

An example where this is needed if the content of the headline is `"⇒"`.
"""
function rewrite_header(context::DocContext, content)
    result = content
    if get_anchor_name(context, result, "") == ""
        # Need to rewrite it (a litte bit)
        result *= " ("
        for (index, character) in enumerate(content)
            if index != 1
                result *= ' '
            end
            result *= "U+" * uppercase(string(
                codepoint(character), base=16, pad=4))
        end
        result *= ')'
    end
    return result
end

const re_substitute_doc = r"^\s*doc:(.*)$"
substitute_obj_header(context::DocContext, md) = nothing
function substitute_obj_header(context::DocContext, v::Vector)
    index = 1
    while index <= length(v)
        part = v[index]
        if part isa Markdown.Header  &&  part.text isa Vector  &&
                length(part.text) == 1  && part.text[1] isa Markdown.Code
            code_md = part.text[1]
            match_obj = match(re_substitute_doc, code_md.code)
            if match_obj !== nothing
                obj_spec = parse_obj_specification(match_obj[1])
                content = String(obj_spec.obj_sym)

                # rewrite (if there is no nice anchor)
                content_in_header = rewrite_header(context, content)
                code_md.code = content_in_header

                doc_ref = DocRef(context.source,
                    get_anchor_name(context, content_in_header))
                append_doc_ref(context.env, obj_spec, doc_ref)

                sub_md = Markdown.parse("[inline](" * match_obj[1] * ")")
                substitute_marker_in_markdown(context, sub_md)

                for elem in reverse(sub_md.content)
                    insert!(v, index+1, elem)
                end
            end
        end
        substitute_obj_header(context, part)
        index += 1
    end
    return nothing
end
substitute_obj_header(context::DocContext, md::MD_has_subparts) =
    substitute_obj_header(context, get_subparts(md))
# }}}

# Subst Code-Cross-References {{{
substitute_code_ref(context::DocContext, md) = nothing
function substitute_code_ref(context::DocContext, v::Vector)
    index = 1
    while index <= length(v)
        part = v[index]
        if part isa Markdown.Code
            obj_spec = nothing
            try
                obj_spec = parse_obj_specification(part.code; default_sig=UnionAll)
            catch
            end
            if obj_spec !== nothing
                rendered = render_ref(context.env, obj_spec; error_if_not_found=false)
                if rendered !== nothing
                    deleteat!(v, index)
                    for elem in reverse(rendered.content)
                        insert!(v, index, elem)
                    end
                end
            end
        end
        substitute_code_ref(context, part)
        index += 1
    end
end
substitute_code_ref(context::DocContext, md::MD_has_subparts) =
    substitute_code_ref(context, get_subparts(md))
# }}}

"""
replace some "markers" inside markdown with actual code from
the file here. This is done to make sure the code in the documentation
is exactly the same code as the one used in the example functions.
"""
function substitute_marker_in_markdown(context::DocContext, md)
    substitute_code_func(context, md)
    substitute_canvas_image(context, md)
    substitute_inline_markdown(context, md)
    substitute_obj_header(context, md)
    return nothing
end

# Fadeout {{{
"""
```
create_fadeout_surface(w, h, dist=30)

w     Integer           width
h     Integer           height
dist  Integer           distance to boundary where to decrease alpha
```

creates a surface where the alpha component is decreased (to 0) near
the boundary.
"""
function create_fadeout_surface(w::Integer, h::Integer, dist::Integer=30)
    w, h, dist = Int32(w), Int32(h), Int32(dist)
    if 2*dist > w || 2*dist > h
        cv_error("d too large: 2*dist > w  or 2*dist > h")
    end
    canvas = CV_Std2DCanvas(w, h)
    Cairo.flush(canvas.surface)
    data = canvas.surface.data
    for x = 1:dist
        p1 = (dist-x)^2
        data[x    , dist+1:h-dist] .= p1
        data[w-x+1, dist+1:h-dist] .= p1
        for y = 1:dist
            p2 = (dist-y)^2
            data[x, y] = data[x, h-y+1] =
            data[w-x+1, y] = data[w-x+1, h-y+1] = UInt32(p1 + p2)
        end
    end
    for y = 1:dist
        p2 = (dist-y)^2
        data[dist+1:w-dist, y] .= p2
        data[dist+1:w-dist, h-y+1] .= p2
    end
    data[dist+1:w-dist, dist+1:h-dist] .= 0x0

    # normalize
    data[:] .= round.(UInt32,
        255*(
            1.0 .- max.(min.(Float64.(data[:])/((dist-3)^2), 1.0), 0.0)
            )) .<< 24

    Cairo.mark_dirty(canvas.surface)
    return canvas
end

"""
default fadeout canvas for `create_doc_icon`.
"""
const cv_default_fadeout_canvas = create_fadeout_surface(70, 70, 20)

"""
```
create_doc_icon(src_canvas, src_rect=src_canvas.bounding_box,
                fadeout_canvas=cv_default_fadeout_canvas)

src_canvas      CV_2DCanvas         canvas with the source image
src_rect        CV_Rectangle{Int32} rectangle to use for icon
fadeout_canvas  CV_2DCanvas         canvas used for fadeout effect
```

creates a canvas with an "icon" for the documentation. The size of the
icon is determined by the size of the `fadeout_canvas`. The `src_rect`
is scaled to the icon size.
"""
function create_doc_icon(
        src_canvas::CV_2DCanvas,
        src_rect::CV_Rectangle{Int32}=src_canvas.bounding_box,
        fadeout_canvas::CV_2DCanvas=cv_default_fadeout_canvas)  # {{{

    w, h = fadeout_canvas.pixel_width, fadeout_canvas.pixel_height
    icon = CV_Std2DCanvas(w, h)

    sw, sh = cv_width(src_rect), cv_height(src_rect)
    cv_create_context(icon) do con
        ctx = con.ctx

        set_operator(ctx, Cairo.OPERATOR_SOURCE)
        set_source_rgb(ctx, 1, 1, 1)
        rectangle(ctx, 0, 0, w, h)
        fill_preserve(ctx)

        pat = CairoPattern(src_canvas.surface)
        pattern_set_filter(pat, Cairo.FILTER_BEST)
        pat_mat = CairoMatrix(sw/w, 0, 0, sh/h, src_rect.left, src_rect.bottom)
        set_matrix(pat, pat_mat)
        set_source(ctx, pat)
        set_operator(ctx, Cairo.OPERATOR_OVER)
        fill_preserve(ctx)

        set_source(ctx, fadeout_canvas.surface)
        set_operator(ctx, Cairo.OPERATOR_DEST_IN)
        fill(ctx)
    end

    return icon
end # }}}

# }}}

# Write short descritions for arguments of methods {{{

mutable struct StringIO
    buffer     :: String
    indent     :: Vector{AbstractString}
    max_width  :: Int
    cur_width  :: Int
end

function StringIO(max_width::Integer, cur_indent::AbstractString="")
    indent = [cur_indent]
    return StringIO("", indent, Int(max_width), 0)
end

function append_newline(sio::StringIO)
    sio.buffer *= "\n" * sio.indent[end]
    sio.cur_width = textwidth(sio.indent[end])
    return sio
end

function push_indent(sio::StringIO, new_indent::AbstractString)
    push!(sio.indent, new_indent)
    return sio
end

function pop_indent(sio::StringIO)
    if length(sio.indent) > 1
        pop!(sio.indent)
    end
    return sio
end

function append_nc(sio::StringIO, to_append::AbstractString)
    sio.buffer *= to_append
    sio.cur_width += textwidth(to_append)
    return sio
end

function append(sio::StringIO, to_append::AbstractString,
        sep_if_fits=""::AbstractString, sep_for_newline=""::AbstractString)
    inline = sep_if_fits * to_append
    inline_width = textwidth(inline)
    if inline_width + sio.cur_width >= sio.max_width
        append_newline(append_nc(sio, sep_for_newline))
        inline = to_append
        inline_width = textwidth(inline)
    end
    append_nc(sio, inline)
    return sio
end

"""
Highly experimental

I use it to get a first version of a function description for the docstring.

Many undocumented methods from Base are used
"""
function get_arg_description(m::Method; indent="")
    args = Base.arg_decl_parts(m)[2]     # args[1]: function itself

    result = StringIO(79, indent)
    append_newline(append_nc(result, "```"))
    header = indent * args[1][2] * '('
    push_indent(result, repeat(' ', textwidth(header)))
    append(result, header)

    isfirst = true
    for arg in args[2:end]
        append(result, arg[1], isfirst ? "" : ", ", isfirst ? "" : ",")
        isfirst = false
    end

    kwargs = Base.kwarg_decl(m)
    if !isempty(kwargs)
        append_newline(append(result, ";"))
        isfirst = true
        for arg in kwargs
            append(result, String(arg), isfirst ? "" : ", ", isfirst ? "" : ",")
            isfirst = false
        end
    end
    pop_indent(result)
    append_newline(append_newline(append(result, ")")))

    max_arg_width = length(args) > 1 ?
        maximum(map(arg -> textwidth(arg[1]), args[2:end])) : 4
    max_type_width = length(args) > 1 ?
        maximum(map(arg -> textwidth(arg[2]), args[2:end])) : 4

    arg_fmt = @sprintf("%%-%is", max_arg_width)
    pad_arg = @eval str -> @sprintf($arg_fmt, str)

    type_fmt = @sprintf("%%-%is", max_type_width)
    pad_type = @eval str -> @sprintf($type_fmt, str)

    for arg in args[2:end]
        append_nc(result, Base.invokelatest(pad_arg ,arg[1]))
        append_nc(result, "   ")
        append_nc(result, Base.invokelatest(pad_type, arg[2]))
        append_newline(result)
    end
    append_newline(append_nc(result, "```"))

    return result.buffer
end

function get_arg_description(something)
    result = ""
    for method in methods(something)
        result *= get_arg_description(method)
        result *= "\n\n"
    end
    return result
end

function list_all_sigs(something)
    for m in methods(something)
        println(tuple(m.sig.types[2:end]...))
    end
end
# }}}

# write Markdown files {{{
"""
write Markdown-files for all Documents.
"""
function write_markdown_documentations(documents::Vector{Document})
    for document in documents
        filename = document.source.filename * ".md"
        open(filename, "w") do fio
            write(fio, string(document.markdown))
            write(fio, "\n\n")
        end
    end
end
# }}}

const re_find_cv_start = r"^(cv|CV)_"

function isless_ignore_cv(sym1::Symbol, sym2::Symbol)
    name1 = replace(lowercase(String(sym1)), re_find_cv_start => "")
    name2 = replace(lowercase(String(sym2)), re_find_cv_start => "")
    return isless(name1, name2)
end

function render_ref(doc_env::DocCreationEnvironment,
        obj_spec::ObjSpecification; error_if_not_found::Bool=true) # {{{
    ref_vec = get_doc_refs(doc_env, obj_spec)
    if ref_vec === nothing
        if error_if_not_found
            error("Cannot find ", obj_spec)
        else
            return nothing
        end
    end

    sym_name = String(obj_spec.obj_sym)
    parts = []
    if length(ref_vec) == 1
        push!(parts, doc_ref_to_markdownlink(ref_vec[1], sym_name))
    else
        push!(parts, Markdown.Code(sym_name))
        push!(parts, " (")
        first = true
        for ref in ref_vec
            if first
                first = false
            else
                push!(parts, ", ")
            end
            push!(parts, doc_ref_to_markdownlink(ref, ref.source.filename))
        end
        push!(parts, ")")
    end
    return Markdown.Paragraph(parts)
end # }}}

function create_letter_list(doc_env::DocCreationEnvironment)
    letter_list = Dict{AbstractChar, Markdown.List}()
    refs = doc_env.refs

    sym_vec = Vector{Symbol}()  # all the symbols (that are documented)
    for sym_dict in values(refs)
        push!(sym_vec, keys(sym_dict)...)
    end
    sort!(sym_vec; lt=isless_ignore_cv)

    for entry in sym_vec
        entry_name = String(entry)
        fchar = first(replace(uppercase(entry_name), re_find_cv_start => ""))
        if !isletter(fchar)
            fchar = '…'
        end
        if !haskey(letter_list, fchar)
            letter_list[fchar] = Markdown.List()
        end

        para = Markdown.Paragraph(Vector{Any}([Markdown.Code(entry_name)]))
        push!(letter_list[fchar].items, para)
    end
    return letter_list
end

function create_index(doc_env::DocCreationEnvironment)
    doc_source = DocSource("Index", @__MODULE__)
    context = DocContext(doc_env, doc_source)

    header = Markdown.Header{1}("Index")

    letter_list = create_letter_list(doc_env)

    md = Markdown.MD(header)
    for fchar in sort(collect(keys(letter_list)))
        hdl = "" * fchar
        md = Markdown.MD(md, Markdown.Header{2}(hdl), letter_list[fchar])
    end

    return md
end

function create_overview(doc_env::DocCreationEnvironment,
        docs::Vector{Document})

    header = Markdown.Header{1}("Overview")

    list = Markdown.List()
    for doc in docs
        content = doc.markdown.content
        if length(content) > 0 && content[1] isa Markdown.Header
            para = Markdown.Paragraph(content[1].text)
            push!(list.items, para)
        end
    end

    return Markdown.MD(header, list)
end

function create_readme(doc_env::DocCreationEnvironment,
        docs::Vector{Document})

    doc_source = DocSource("README", @__MODULE__)

    md = Markdown.MD([create_overview(doc_env, docs), create_index(doc_env)])
    return Document(doc_source, md)
end

function create_documents()
    doc_env = DocCreationEnvironment()

    documents = Vector{Document}()
    for mod in doc_modules
        println("generating markdown for ", mod)
        push!(documents, mod.create_document(doc_env))
    end
    println("generating markdown for Readme")
    push!(documents, create_readme(doc_env, documents))

    # Cross-Ref
    println("inserting cross references")
    context = DocContext(doc_env, DocSource("cross-refs", @__MODULE__))
    for doc in documents
        substitute_code_ref(context, doc.markdown)
    end

    write_markdown_documentations(documents)
    return nothing
end

for filename in (
    "./Canvas.jl",
    "./Context.jl",
    "./Style.jl",
    "./PixelCoordinates.jl",
    "./Painter.jl",
    "./Axis.jl")
    "./LayoutTutorial.jl",
    push!(doc_modules, include(filename))
end

# create_documents()

end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
"Computes the hamming distance between two strands."
function distance(a, b)
    if length(a) != length(b)
        throw(ArgumentError("Strands $a and $b are not of same length"))
    end
    diff = 0
    for (idx, nucleotide) in enumerate(a)
        if nucleotide != b[idx]
            diff += 1
        end
    end
    return diff
end

ext: jl

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
using Documenter
using Mutabilities

makedocs(
    sitename = "Mutabilities",
    format = Documenter.HTML(),
    modules = [Mutabilities],
)

deploydocs(repo = "github.com/tkf/Mutabilities.jl", push_preview = true)


text:
using IterativeSolvers
using LinearMaps

"""
Implementation based on the presentation in
Holtz, Sebastian, Thorsten Rohwedder, and Reinhold Schneider. "The alternating linear scheme for tensor optimization in the tensor train format." SIAM Journal on Scientific Computing 34.2 (2012): A683-A713.
"""

function updateH_mals!(x_vec::Array{T,3}, A_vec::Array{T,4}, Hi::AbstractArray{T,5}, Him::AbstractArray{T,5}) where T<:Number
	@tensor(Him[a,i,α,l,β] = conj.(x_vec)[j,α,x]*(Hi[z,j,x,k,y]*x_vec[k,β,y])*A_vec[i,l,a,z]) #size(rAim,ri,ni,ri,ni)
	nothing
end

function init_H_mals(x_tt::ttvector{T},A::ttoperator{T},rmax::Int) where T<:Number
	d = length(x_tt.ttv_dims)
	H = Array{Array{T,5}}(undef, d-1)
	H[d-1] = reshape(permutedims(A.tto_vec[d],[3,1,2,4]), :, x_tt.ttv_dims[d], 1, x_tt.ttv_dims[d], 1) #size(R^A_d, n_d, r_{d+1}, n_d, r_{d+1})
	for i = (d-1) : -1 : 2
		# Update H[i-1]
		rmax_i = min(rmax,prod(x_tt.ttv_dims[1:i]),prod(x_tt.ttv_dims[i+1:end]))
		H[i-1] = zeros(T,A.tto_rks[i],x_tt.ttv_dims[i],rmax_i,x_tt.ttv_dims[i],rmax_i)
		Hi = @view(H[i][:,:,1:x_tt.ttv_rks[i+2],:,1:x_tt.ttv_rks[i+2]])
		Him = @view(H[i-1][:,:,1:x_tt.ttv_rks[i+1],:,1:x_tt.ttv_rks[i+1]])
		updateH_mals!(x_tt.ttv_vec[i+1], A.tto_vec[i], Hi, Him)
	end
	return H
end

function updateHb_mals!(xtt_vec::Array{T,3}, btt_vec::Array{T,3}, Hbi::AbstractArray{T,3}, Hbim::AbstractArray{T,3}) where T<:Number
	@tensor Hbim[β,i,χ] = conj.(xtt_vec)[j,χ,a]*Hbi[γ,j,a]*btt_vec[i,β,γ]
	nothing
end

function init_Hb_mals(x_tt::ttvector{T},b::ttvector{T},rmax::Int) where T<:Number
	d = length(x_tt.ttv_dims)
	Hb = Array{Array{T,3}}(undef, d-1)
	Hb[d-1] = reshape(permutedims(b.ttv_vec[d],[2,1,3]),b.ttv_rks[d],b.ttv_dims[d],1)
	for i in d-1:-1:2
		rmax_i = min(rmax,prod(x_tt.ttv_dims[1:i]),prod(x_tt.ttv_dims[i+1:end]))
		Hb[i-1] = zeros(T,b.ttv_rks[i+1],b.ttv_dims[i+1],rmax_i)
		Hbi = @view(Hb[i][:,:,1:x_tt.ttv_rks[i+2]])
		Hbim = @view(Hb[i-1][:,:,1:x_tt.ttv_rks[i+1]])
		updateHb_mals!(x_tt.ttv_vec[i+1],b.ttv_vec[i],Hbi,Hbim)
	end
	return Hb
end

function left_core_move_mals(xtt::ttvector{T},i::Integer,V::Array{T,4},tol::Float64,rmax::Integer) where T<:Number
	# Perform the truncated svd
	u_V, s_V, v_V, = svd(reshape(V, prod(size(V)[1:2]), :))
	# Determine the truncated rank
	s_trunc = sv_trunc(s_V,tol)
	# Update the ranks to the truncated one
	xtt.ttv_rks[i+1] = min(length(s_trunc),rmax)
	println("Rank: $(xtt.ttv_rks[i+1]),	Max rank=$rmax")
	println("Discarded weight: $((norm(s_V)-norm(s_V[1:xtt.ttv_rks[i+1]]))/norm(s_V))")

	# xtt.ttv_vec[i+1] = truncated Transpose(v_V)
	xtt.ttv_vec[i+1] = permutedims(reshape(v_V'[1:xtt.ttv_rks[i+1],:],xtt.ttv_rks[i+1],size(V,3),size(V,4)),[2 1 3])
	# xtt.ttv_vec[i] = truncated u_V * Diagonal(s_V)
	xtt.ttv_vec[i] = reshape(u_V[:, 1:xtt.ttv_rks[i+1]] * Diagonal(s_trunc[1:xtt.ttv_rks[i+1]]),size(V,1),size(V,2),:)
	xtt.ttv_ot[i+1] = 1
	xtt.ttv_ot[i] = 0
	return xtt
end

function right_core_move_mals(xtt::ttvector{T},i::Integer,V::Array{T,4},tol::Float64,rmax::Integer) where T<:Number
	# Perform the truncated svd
	u_V, s_V, v_V, = svd(reshape(V, prod(size(V)[1:2]), :))
	# Determine the truncated rank
	s_trunc = sv_trunc(s_V,tol)
	# Update the ranks to the truncated one
	xtt.ttv_rks[i+1] = min(length(s_trunc),rmax)
	println("Rank: $(xtt.ttv_rks[i+1]),	Max rank=$rmax")
	println("Discarded weight: $((norm(s_V)-norm(s_V[1:xtt.ttv_rks[i+1]]))/norm(s_V))")

	# xtt.ttv_vec[i] = truncated u_V
	xtt.ttv_vec[i] = reshape(u_V[:, 1:xtt.ttv_rks[i+1]], size(V,1), size(V,2), xtt.ttv_rks[i+1])
	xtt.ttv_ot[i] = -1

	# xtt.ttv_vec[i+1] = truncated Diagonal(s_V) * Transpose(v_V)
	xtt.ttv_vec[i+1] = permutedims(reshape(Diagonal(s_trunc[1:xtt.ttv_rks[i+1]]) * v_V'[1:xtt.ttv_rks[i+1],:],xtt.ttv_rks[i+1],size(V,3),size(V,4)), [2 1 3])
	xtt.ttv_ot[i+1] = 0
	return xtt
end

function K_full_mals(Gi::AbstractArray{T,5},Hi::AbstractArray{T,5},K_dims::NTuple{4,Int}) where T<:Number
	K = zeros(T,prod(K_dims),prod(K_dims))
	Krshp = reshape(K,(K_dims...,K_dims...))
	@tensor Krshp[a,b,c,d,e,f,g,h] = Gi[a,b,e,f,z]*Hi[z,c,d,g,h] #size(K)=(ni,rim,nip,rip)
	return Hermitian(K)
end

function Ksolve_mals(Gi::AbstractArray{T,5}, Hi::AbstractArray{T,5}, G_bi::AbstractArray{T,3}, H_bi::AbstractArray{T,3}) where T<:Number
	K_dims = (size(Gi,1),size(Gi,2),size(Hi,2),size(Hi,3))
	K = K_full_mals(Gi,Hi,K_dims)
	Pb = zeros(T,K_dims)
	@tensor Pb[a,b,c,d] = G_bi[a,b,z]*H_bi[z,c,d] #size(ni,rim,nip,rip)
	V = reshape(K,prod(K_dims),:) \ Pb[:]
	return reshape(V,K_dims)
end

function K_eigmin_mals(Gi::Array{T,5},Hi::Array{T,5},ttv_vec_i::Array{T,3},ttv_vec_ip::Array{T,3};it_solver=false,itslv_thresh=256::Int64,maxiter=200::Int64,tol=1e-6::Float64) where T<:Number
	K_dims = (size(ttv_vec_i,1),size(ttv_vec_i,2),size(ttv_vec_ip,1),size(ttv_vec_ip,3))
	Gtemp = @view(Gi[:,1:K_dims[2],:,1:K_dims[2],:])
	Htemp = @view(Hi[:,:,1:K_dims[4],:,1:K_dims[4]])
	if it_solver || prod(K_dims) > itslv_thresh
		H = zeros(T,prod(K_dims))
		function K_matfree(V::AbstractArray{S,1};K_dims=K_dims::NTuple{4,Int},H=H::AbstractArray{S,1},Gtemp=Gtemp::AbstractArray{S,5},Htemp=Htemp::AbstractArray{S,5}) where S<:Number
			Hrshp = reshape(H,K_dims)
			@tensoropt((f,h), Hrshp[a,b,c,d] = Gtemp[a,b,e,f,z]*reshape(V,K_dims)[e,f,g,h]*Htemp[z,c,d,g,h])
			return H::AbstractArray{S,1}
		end
		X0 = zeros(T,prod(K_dims))
		X0_temp = reshape(X0,K_dims)
		@tensor X0_temp[a,b,c,d] = ttv_vec_i[a,b,z]*ttv_vec_ip[c,z,d]
		r = lobpcg(LinearMap(K_matfree,prod(K_dims);ishermitian = true),false,X0,1;maxiter=maxiter,tol=tol)
		return r.λ[1]::Float64, reshape(r.X[:,1],K_dims)::Array{T,4}
	else
		K = K_full_mals(Gtemp,Htemp,K_dims)
		F = eigen(K,1:1)
		return real(F.values[1])::Float64,reshape(F.vectors[:,1],K_dims)::Array{T,4}
	end
end

"""
Returns the solution `tt_opt :: ttvector` of Ax=b using the MALS algorithm where A is given as `ttoperator` and `b`, `tt_start` are `ttvector`.
The ranks are adapted at each microstep by keeping the singular values larger than `tol`.
"""
function mals_linsolv(A :: ttoperator{T}, b :: ttvector{T}, tt_start :: ttvector{T}; tol=1e-12::Float64,rmax=round(Int,sqrt(prod(tt_start.ttv_dims)))) where T<:Number
	# mals finds the minimum of the operator J(x)=1/2*<Ax,x> - <x,b>
	# input:
	# 	A: the tensor operator in its tensor train format
	#   b: the tensor in its tensor train format
	#	tt_start: start value in its tensor train format
	#	opt_rks: rank vector considered to be optimal enough
	#	tol: tolerated inaccuracy
	# output:
	#	tt_opt: stationary point of J up to tolerated inaccuracy
	# 			in its tensor train format

	# Initialize the to be returned tensor in its tensor train format
	tt_opt = orthogonalize(tt_start)
	dims = tt_start.ttv_dims
	d = length(dims)
	# Define the array of ranks of A [R_0=1,R_1,...,R_d]
	A_rks = A.tto_rks
	# Define the array of ranks of b [R^b_0=1,R^b_1,...,R^b_d]
	b_rks = b.ttv_rks

	# Initialize the arrays of G and H
	G = Array{Array{T,5}}(undef, d)
	# Initialize the arrays of G_b and H_b
	G_b = Array{Array{T,3}}(undef, d)
	# Initialize G[1], G_b[1], H[d] and H_b[d]
	for i in 1:d
		rmax_i = min(rmax,prod(dims[1:i-1]),prod(dims[i:end]))
		G[i] = zeros(dims[i],rmax_i,dims[i],rmax_i,A_rks[i+1])
		G_b[i] = zeros(dims[i],rmax_i,b_rks[i+1])
	end
	G[1][:,1:1,:,1:1,:] = reshape(A.tto_vec[1][:,:,1,:], dims[1],1,dims[1], 1, :)
	G_b[1] = reshape(b.ttv_vec[1], dims[1], 1, :)

	H = init_H_mals(tt_opt,A,rmax)
	H_b = init_Hb_mals(tt_opt,b,rmax)

	#while 1==1 #TODO make it work for real
		# First half sweep
		for i = 1:(d-1)
			Gi = @view(G[i][:,1:tt_opt.ttv_rks[i],:,1:tt_opt.ttv_rks[i],:])
			Hi = @view(H[i][:,:,1:tt_opt.ttv_rks[i+2],:,1:tt_opt.ttv_rks[i+2]])
			G_bi = @view(G_b[i][:,1:tt_opt.ttv_rks[i],:])
			H_bi = @view(H_b[i][:,:,1:tt_opt.ttv_rks[i+2]])
			# Define V as solution of K*x=P2b in x
			V = Ksolve_mals(Gi,Hi,G_bi,H_bi)
			tt_opt = right_core_move_mals(tt_opt,i,V,tol,rmax)
			# Update G[i+1],G_b[i+1]
			Gip = @view(G[i+1][:,1:tt_opt.ttv_rks[i+1],:,1:tt_opt.ttv_rks[i+1],:])
			G_bip = @view(G_b[i+1][:,1:tt_opt.ttv_rks[i+1],:])
			update_G!(tt_opt.ttv_vec[i],A.tto_vec[i+1],Gi,Gip)
			update_Gb!(tt_opt.ttv_vec[i],b.ttv_vec[i+1],G_bi,G_bip)
		end

		# Second half sweep
		for i = d-1:(-1):1
			Gi = @view(G[i][:,1:tt_opt.ttv_rks[i],:,1:tt_opt.ttv_rks[i],:])
			Hi = @view(H[i][:,:,1:tt_opt.ttv_rks[i+2],:,1:tt_opt.ttv_rks[i+2]])
			G_bi = @view(G_b[i][:,1:tt_opt.ttv_rks[i],:])
			H_bi = @view(H_b[i][:,:,1:tt_opt.ttv_rks[i+2]])
			# Define V as solution of K*x=P2b in x
			V = Ksolve_mals(Gi,Hi,G_bi,H_bi)
			tt_opt = left_core_move_mals(tt_opt,i,V,tol,rmax)
			# Update H[i-1], H_b[i-1]
			if i > 1
				Him = @view(H[i-1][:,:,1:tt_opt.ttv_rks[i+1],:,1:tt_opt.ttv_rks[i+1]])
				updateH_mals!(tt_opt.ttv_vec[i+1], A.tto_vec[i], Hi, Him)
				H_bim = @view(H_b[i-1][:,:,1:tt_opt.ttv_rks[i+1]])
				updateHb_mals!(tt_opt.ttv_vec[i+1], b.ttv_vec[i], H_bi,H_bim)
			end
		end
	#end
	return tt_opt
end

"""
Returns the list of the approximate smallest eigenvalue at each microstep, the corresponding eigenvector as a `ttvector` and the list of the maximum rank at each microstep.

`A` is given as `ttoperator` and `tt_start` is a `ttvector`.
The ranks are adapted at each microstep by keeping the singular values larger than `tol`.
The number of total sweeps is given by `sweep_schedule[end]`. The maximum rank is prescribed at each sweep `sweep_schedule[k]≤ i <sweep_schedule[k+1]` by `rmax_schedule[k]`.
"""
function mals_eigsolv(A :: ttoperator{T}, tt_start :: ttvector{T}; tol=1e-12::Float64,sweep_schedule=[2]::Array{Int64,1},rmax_schedule=[round(Int,sqrt(prod(tt_start.ttv_dims)))]::Array{Int64,1},it_solver=false::Bool,linsolv_maxiter=200::Int64,linsolv_tol=max(sqrt(tol),1e-8)::Float64,itslv_thresh=256::Int) where T<:Number

	@assert(length(rmax_schedule)==length(sweep_schedule),"Sweep schedule error")
	# Initialize the to be returned tensor in its tensor train format
	tt_opt = orthogonalize(tt_start)
	dims = tt_start.ttv_dims
	d = length(dims)
	# Initialize the output objects
	E = Float64[]
	r_hist = Int64[]
	# Initialize the arrays of G
	G = Array{Array{T}}(undef, d)
	rmax = maximum(rmax_schedule)
	# Initialize G[i]
	for i in 1:d
		rmax_i = min(rmax,prod(dims[1:i-1]),prod(dims[i:end]))
		G[i] = zeros(dims[i],rmax_i,dims[i],rmax_i,A.tto_rks[i+1])
	end
	G[1][:,1:1,:,1:1,:] = reshape(A.tto_vec[1][:,:,1,:], dims[1],1,dims[1], 1, :)

	H = init_H_mals(tt_opt,A,rmax)

	nsweeps = 0 #sweeps counter
	i_schedule = 1
	while i_schedule <= length(sweep_schedule)
		nsweeps+=1
		println("Macro-iteration $nsweeps; bond dimension $(rmax_schedule[i_schedule])")

		if nsweeps == sweep_schedule[i_schedule]
			i_schedule+=1
			if i_schedule > length(sweep_schedule)
				return E::Array{Float64,1}, tt_opt::ttvector{T}, r_hist::Array{Int,1}
			end
		end

		# First half sweep
		for i = 1:(d-1)
			println("Forward sweep: core optimization $i out of $(d-1)")
			# Define V as solution of K V= λ V for smallest λ
			λ,V = K_eigmin_mals(G[i],H[i],tt_opt.ttv_vec[i],tt_opt.ttv_vec[i+1];it_solver=it_solver,maxiter=linsolv_maxiter,tol=linsolv_tol)
			println("Eigenvalue: $λ")
			E = vcat(E,λ)
			tt_opt = right_core_move_mals(tt_opt,i,V,tol,rmax_schedule[i_schedule])
			r_hist = vcat(r_hist,maximum(tt_opt.ttv_rks))
			# Update G[i+1],G_b[i+1]
			Gi = @view(G[i][:,1:tt_opt.ttv_rks[i],:,1:tt_opt.ttv_rks[i],:])
			Gip = @view(G[i+1][:,1:tt_opt.ttv_rks[i+1],:,1:tt_opt.ttv_rks[i+1],:])
			update_G!(tt_opt.ttv_vec[i],A.tto_vec[i+1],Gi,Gip)
		end

		# Second half sweep
		for i = d-1:(-1):1
			println("Backward sweep: core optimization $(d-i) out of $(d-1)")
			# Define V as solution of K*x=P2b in x
			λ,V = K_eigmin_mals(G[i],H[i],tt_opt.ttv_vec[i],tt_opt.ttv_vec[i+1];it_solver=it_solver,maxiter=linsolv_maxiter,tol=linsolv_tol)
			println("Eigenvalue: $λ")
			E = vcat(E,λ)
			tt_opt = left_core_move_mals(tt_opt,i,V,tol,rmax_schedule[i_schedule])
			r_hist = vcat(r_hist,maximum(tt_opt.ttv_rks))
			# Update H[i-1]
			if i > 1
				Hi = @view(H[i][:,:,1:tt_opt.ttv_rks[i+2],:,1:tt_opt.ttv_rks[i+2]])
				Him = @view(H[i-1][:,:,1:tt_opt.ttv_rks[i+1],:,1:tt_opt.ttv_rks[i+1]])
				updateH_mals!(tt_opt.ttv_vec[i+1], A.tto_vec[i], Hi, Him)
			end
		end
	end
	return E::Array{Float64,1}, tt_opt::ttvector{T}, r_hist::Array{Int,1}
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using Nextion
using Nextion.Event: TouchEvent, CurrentPageIDHeadEvent, PositionHeadEvent, SleepPositionHeadEvent, StringHeadEvent, NumberHeadEvent
using Nextion.Event.Touch
using Test

@testset "Events" begin

    @testset "TouchEvent" begin

        @testset "Constants" begin
            @test Integer(Touch.Press) == 0x01
            @test Integer(Touch.Release) == 0x00
            @test Touch.code(0x01) == Touch.Press
            @test Touch.code(0x00) == Touch.Release
        end

        @testset "Event" begin
            msg = [0x65, 0x00, 0x02, 0x01, 0xff, 0xff, 0xff]
            evt = TouchEvent(msg)
            @test evt.code == Return.Code.EVENT_TOUCH_HEAD
            @test evt.pid == PageID(0x00)
            @test evt.cid == ComponentID(0x02)
            @test evt.tevts == Touch.Press  # touch event state
        end
    end

    @testset "CurrentPageIDHeadEvent" begin
        msg = [0x66, 0x02, 0xff, 0xff, 0xff]
        evt = CurrentPageIDHeadEvent(msg)
        @test evt.code == Return.Code.CURRENT_PAGE_ID_HEAD
        @test evt.pid == PageID(0x02)
    end

    @testset "PositionHeadEvent" begin
        msg = [0x67, 0x00, 0x7a, 0x00, 0x1e, 0x01, 0xff, 0xff, 0xff]
        evt = PositionHeadEvent(msg)
        @test evt.code == Return.Code.EVENT_POSITION_HEAD
        @test evt.x == UInt16(122)
        @test evt.y == UInt16(30)
        @test evt.tevts == Touch.Press
    end

    @testset "SleepPositionHeadEvent" begin
        msg = [0x68, 0x00, 0x7a, 0x00, 0x1e, 0x01, 0xff, 0xff, 0xff]
        evt = SleepPositionHeadEvent(msg)
        @test evt.code == Return.Code.EVENT_SLEEP_POSITION_HEAD
        @test evt.x == UInt16(122)
        @test evt.y == UInt16(30)
        @test evt.tevts == Touch.Press
    end

    @testset "StringHeadEvent" begin
        msg = [0x70, 0x61, 0x62, 0x63, 0xff, 0xff, 0xff]
        evt = StringHeadEvent(msg)
        @test evt.code == Return.Code.STRING_HEAD
        @test evt.value == "abc"
    end

    @testset "NumberHeadEvent" begin
        msg = [0x71, 0x66, 0x00, 0x00, 0x00, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.code == Return.Code.NUMBER_HEAD
        @test evt.value == 102

        msg = [0x71, 0x66, 0x01, 0x00, 0x00, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0x00000166 # 102 + 256

        msg = [0x71, 0x01, 0xff, 0x00, 0x00, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0x0000ff01 # 65281

        msg = [0x71, 0x01, 0x00, 0xff, 0x00, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0x00ff0001 # 16711681

        msg = [0x71, 0x01, 0x00, 0x00, 0xff, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0xff000001

        msg = [0x71, 0xff, 0xff, 0xff, 0x7f, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == typemax(Int32)

        msg = [0x71, 0xfe, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0xfffffffe
        @test evt.signedvalue == -2

        msg = [0x71, 0x00, 0x00, 0x00, 0x80, 0xff, 0xff, 0xff]
        evt = NumberHeadEvent(msg)
        @test evt.value == 0x80000000
        @test evt.signedvalue == typemin(Int32)

    end

end


text:
#
# This function computes the time auto-correlation function of a vector
#

using LinearAlgebra

function tcorr(x_in :: Vector{Float64})

  x = copy(x_in)
  n = size(x)[1]

  # Centralizing at the average

  xav = sum(x)/n
  @. x = x - xav

  # Computing the time-dependent correlation function

  tcorr = zeros(Float64,n)
  xnorm1 = 0.
  xnorm2 = 0.
  @inbounds for dt in n-1:-1:0
    xnorm1 = xnorm1 + x[n-dt]^2
    xnorm2 = xnorm2 + x[dt+1]^2
    if xnorm1 != 0 || xnorm2 != 0
      tcorr[dt+1] = LinearAlgebra.dot(@view(x[1:n-dt]),@view(x[dt+1:n]))
      tcorr[dt+1] = tcorr[dt+1] / sqrt(xnorm1*xnorm2)
    else
      tcorr[dt+1] = 0.
    end
  end

  lags = collect(0:n-1)
  return lags, tcorr

end

#
# This function computes the time correlation function of two vectors
#

function tcorr(x_in :: Vector{Float64}, y_in :: Vector{Float64})

  x = copy(x_in)
  y = copy(y_in)

  n = size(x)[1]
  if size(y)[1] != n
    error(" In tcorr both vectors are expected to have the same dimensions. ")
  end

  # Centralizing at the average

  xav = 0.
  yav = 0.
  xav = sum(x)/n
  yav = sum(y)/n
  @. x = x - xav
  @. y = y - yav

  # Computing correlation

  tcorr = zeros(Float64,n)
  xnorm = 0.
  ynorm = 0.
  @inbounds for dt in n-1:-1:0
    xnorm = xnorm + x[n-dt]^2
    ynorm = ynorm + y[dt+1]^2
    if xnorm1 != 0 || xnorm2 != 0
      tcorr[dt+1] = LinearAlgebra.dot(@view(x[1:n-dt]),@view(y[dt+1:n]))
      tcorr[dt+1] = tcorr[dt+1] / sqrt(xnorm*ynorm)
    else
      tcorr[dt+1] = 0.
    end
  end

  lags = collect(0:n-1)
  return lags, tcorr

end


text:
using Test

@testset "GraphBuilding" begin
    adj, els = build_graph(
        abspath(@__DIR__, "..", "test_data", "strucs", "mp-195.cif"),
        use_voronoi = true,
    )
    wm_true = [0.0 1.0 1.0 1.0; 1.0 0.0 1.0 1.0; 1.0 1.0 0.0 1.0; 1.0 1.0 1.0 0.0]
    @test adj == wm_true
    @test els == ["Ho", "Pt", "Pt", "Pt"]

    adj, els = build_graph(
        abspath(@__DIR__, "..", "test_data", "strucs", "mp-195.cif");
        use_voronoi = false,
    )
    @test adj == wm_true
    @test els == ["Ho", "Pt", "Pt", "Pt"]

    # tests for some other file formats
    info = Tuple{Matrix,Vector{String}}[]
    for fp in ["mp-195.poscar", "mp-195.traj", "mp-195.xyz"]
        push!(info, build_graph(abspath(@__DIR__, "..", "test_data", "strucs", fp)))
    end
    for t in info
        @test t[1] == wm_true
        @test t[2] == els
    end

    # test for nonperiodic system
    @test_logs (
        :warn,
        "Voronoi edge weights are not supported if any direction in the structure is nonperiodic. Using cutoff weights method...",
    ) build_graph(
        abspath(@__DIR__, "..", "test_data", "strucs", "methane.xyz"),
        use_voronoi = true,
    )
    adj, els = build_graph(abspath(@__DIR__, "..", "test_data", "strucs", "methane.xyz"))
    @test all(isapprox.(adj[2:5, 1], 1.0, atol = 1e-4))
    @test all(isapprox.(adj[3:2, 2], 0.375, atol = 1e-5))
end


text:
using Documenter, VideoIO

makedocs(
    modules = [VideoIO],
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"],
        analytics = "UA-143027902-2",
        canonical = "https://juliaio.github.io/VideoIO.jl/stable/",
    ),
    sitename = "VideoIO.jl",
    pages = Any[
        "Introduction"=>"index.md",
        "Reading Videos"=>"reading.md",
        "Writing Videos"=>"writing.md",
        "Utilities"=>"utilities.md",
        "Low Level Functionality"=>"lowlevel.md",
        "Index"=>"functionindex.md",
    ],
)
deploydocs(repo = "github.com/JuliaIO/VideoIO.jl.git", push_preview = true)


text:
@testset "Genus" begin
  Qx, x = QQ["x"]
  K, a = NumberField(x^2 - 2, "a")
  OK = maximal_order(K)
  Kt, t  = K["t"]

  E1, b1 = NumberField(t^2 - a, "b1") # ramified at 2
  E2, b2 = NumberField(t^2 - 5, "b2") # unramified at 2

  p = prime_decomposition(OK, 2)[1][1]
  # ramified & dyadic
  g = genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :det)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  g = genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  @test_throws ErrorException genus(HermLat, E1, p, [(0, 1, 1, 1), (2, 2, -1, 0)], type = :det)
  @test_throws ErrorException genus(HermLat, E1, p, [(0, 1, 1, 1), (2, 2, -1, 0)], type = :disc)
  @test_throws ErrorException genus(HermLat, E1, p, [(0, 1, 1, 0), (2, 2, -1, 1)], type = :bla)

  # unramified & dyadic
  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :det)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  @test_throws ErrorException genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :bla)

  # ramified & non-dyadic
  p = prime_decomposition(OK, 5)[1][1]
  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :det)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  g = genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :disc)
  @test sprint(show, "text/plain", g) isa String
  @test sprint(show, g) isa String
  @test_throws ErrorException genus(HermLat, E2, p, [(0, 1, 1), (2, 2, -1)], type = :bla)
end


text:
name(dataset::DatasetName) = lowercase(String(nameof(typeof(dataset))))
url(dataset::DatasetName) = error("URL address not specified for $dataset dataset.")
checksum(::DatasetName) = ""
size(::DatasetName) = error("Size not specified for $dataset dataset.")
problem(::DatasetName) = error("Type of problem not specified for $dataset dataset.")
message(::DatasetName) = ""
function prep(::DatasetName) end

"""
    registering(dsName::DatasetName)

Create a registration block for DataDeps package.
"""
function registering(dsName::DatasetName)
    DataDeps.register(DataDep(
        name(dsName),
        """
            Dataset: $(name(dsName))
            Website: $(url(dsName))
            $(message(dsName))
        """,
        url(dsName),
        checksum(dsName),
        post_fetch_method = prep(dsName)
    ))
end

"""
    extract(path)

Extracts file using `DataDeps.unpack` function and returns path to the new extracted file.
"""
function extract(path)
    newPath = splitext(path)[1]
    DataDeps.unpack(path)
    return newPath
end

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
using SpectralGaussianSimulation
using GeoStatsBase
using Variography
using Plots; gr(size=(900,300))
using VisualRegressionTests
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

@testset "SpectralGaussianSimulation.jl" begin
  if visualtests
    problem = SimulationProblem(RegularGrid{Float64}(100,100), :z => Float64, 3)

    Random.seed!(2019)
    solver = SpecGaussSim(:z => (variogram=GaussianVariogram(range=10.),))
    solution = solve(problem, solver)
    @plottest plot(solution) joinpath(datadir,"isotropic.png") !istravis

    Random.seed!(2019)
    solver = SpecGaussSim(:z => (variogram=GaussianVariogram(distance=Ellipsoidal([20.,5.],[0.])),))
    solution = solve(problem, solver)
    @plottest plot(solution) joinpath(datadir,"anisotropic.png") !istravis
  end
end


text:
using ArgParse;
using Distributions;
using LinearAlgebra;
using Logging;
using DelimitedFiles;
using ProfileView;
using DecFP;
using Quadmath;

include(joinpath("inc", "eap_chain.jl"));
include(joinpath("inc", "average.jl"));
include(joinpath("inc", "acceptance.jl"));

function flip_n!(ψs::Vector)
  ψs[1] += π;
  ψs[2] = π - ψs[2];
end

s = ArgParseSettings();
@add_arg_table! s begin
  "--E0", "-e"
    help = "magnitude of electric field"
    arg_type = Float64
    default = 0.0
  "--chain-type", "-T"
    help = "chain type (dielectric|polar)"
    arg_type = String
    default = "dielectric"
  "--K1", "-J"
    help = "dipole susceptibility along the monomer axis (dielectric chain)"
    arg_type = Float64
    default = 1.0
  "--K2", "-K"
    help = "dipole susceptibility orthogonal to the monomer axis (dielectric chain)"
    arg_type = Float64
    default = 0.0
  "--mu", "-m"
    arg_type = Float64
    default = 1e-2
    help = "dipole magnitude (electret chain)"
  "--energy-type", "-u"
    help = "energy type (noninteracting|interacting)"
    arg_type = String
    default = "noninteracting"
  "--kT", "-k"
    help = "dimensionless temperature"
    arg_type = Float64
    default = 1.0
  "--ensemble-type", "-E"
    help = "ensemble type (force|end-to-end)"
    arg_type = String
    default = "force"
  "--Fz", "-F"
    help = "force in the z-direction (direction of E-field; force ensemble)"
    arg_type = Float64
    default = 0.0
  "--Fx", "-G"
    help = "force in the x-direction (force ensemble)"
    arg_type = Float64
    default = 0.0
  "--rz", "-z"
    help = "end-to-end vector in the z-direction (direction of E-field; etoe ensemble)"
    arg_type = Float64
    default = 0.0
  "--rx", "-x"
    help = "end-to-end vector in the x-direction (etoe ensemble)"
    arg_type = Float64
    default = 0.0
  "--mlen", "-b"
    help = "monomer length"
    arg_type = Float64
    default = 1.0
  "--num-monomers", "-n"
    help = "number of monomers"
    arg_type = Int
    default = 100
  "--num-steps", "-N"
    help = "number of steps"
    arg_type = Int
    default = convert(Int, 1e5)
  "--num-inits", "-M"
    help = "number of random initializations"
    arg_type = Int
    default = 1
  "--force-init", "-I"
    help = "force each random initialization (false to use metro.)"
    action = :store_true
  "--phi-step", "-p"
    help = "maximum ϕ step length"
    arg_type = Float64;
    default = 3*π / 8;
  "--do-flips"
    help = "trial moves with flipping monomers"
    action = :store_true
  "--theta-step", "-q"
    help = "maximum θ step length"
    arg_type = Float64;
    default = 3*π / 16;
  "--chain-frac-step", "-f"
    help = "fraction of monomers to step (end-to-end ensemble)"
    arg_type = Float64
    default = 0.15
  "--step-adjust-lb", "-L"
    help = "adjust step sizes if acc. ratio below this threshold"
    arg_type = Float64
    default = 0.15
  "--step-adjust-ub", "-U"
    help = "adjust step sizes if acc. ratio above this threshold"
    arg_type = Float64
    default = 0.55
  "--step-adjust-scale", "-A"
    help = "scale factor for adjusting step sizes (> 1.0)"
    arg_type = Float64
    default = 1.1
  "--steps-per-adjust", "-S"
    help = "steps between step size adjustments"
    arg_type = Int
    default = 2500
  "--acc", "-a"
    help = "acceptance function (metropolis|kawasaki)"
    arg_type = String
    default = "metropolis"
  "--umbrella-sampling", "-B"
    help = "use umbrella sampling (w/ electrostatic weight function)"
    action = :store_true
  "--update-freq"
    help = "update frequency (seconds)"
    arg_type = Float64;
    default = 15.0;
  "--verbose", "-v"
    help = "verbosity level: 0-nothing, 1-errors, 2-warnings, 3-info"
    arg_type = Int
    default = 3
  "--prefix", "-P"
    help = "prefix for output files"
    arg_type = String
    default = "eap-mcmc"
  "--postfix", "-Q"
    help = "postfix for output files"
    arg_type = String
    default = ""
  "--stepout", "-s"
    help = "steps between storing microstates"
    arg_type = Int
    default = 500
  "--numeric-type"
    help = "numerical data type for averaging (float64|float128|dec128|big)"
    arg_type = String
    default = "float64"
  "--profile", "-Z"
    help = "profile the program"
    action = :store_true
end

pargs = parse_args(s);

if pargs["verbose"] == 3
  global_logger(ConsoleLogger(stderr, Logging.Info));
elseif pargs["verbose"] == 2
  global_logger(ConsoleLogger(stderr, Logging.Warn));
elseif pargs["verbose"] == 1
  global_logger(ConsoleLogger(stderr, Logging.Error));
else
  global_logger(Logging.NullLogger());
end

if pargs["ensemble-type"] == "end-to-end"
  @warn "'end-to-end' ensemble is an experimental option; it has not been validated.";
end

function mcmc(nsteps::Int, pargs)
  ϕstep, θstep = pargs["phi-step"], pargs["theta-step"];
  dϕ_dist = Uniform(-ϕstep, ϕstep);
  dθ_dist = Uniform(-θstep, θstep);
  chain = EAPChain(pargs);
  chain.U = U(chain);
  wf = AntiDipoleWeightFunction(chain);
  acceptor = if pargs["acc"] == "metropolis"
    Metropolis(
               chain,
               (pargs["umbrella-sampling"]) ? wf : WeightlessFunction();
              );
  else
    error("'$(pargs["acc"])' acceptance criteria has not yet been implemented.");
  end
  numeric_type = if pargs["numeric-type"] == "float64"
    Float64;
  elseif pargs["numeric-type"] == "float128"
    Float128;
  elseif pargs["numeric-type"] == "dec128"
    Dec128;
  elseif pargs["numeric-type"] == "big"
    BigFloat;
  else
    error("numeric-type '$(pargs["numeric-type"])' not understood");
    exit(-1);
  end

  force_ensemble_flag = pargs["ensemble-type"] == "force";
  r0 = [pargs["rx"]; 0.0; pargs["rz"]];

  Avg, avgcons, avgconss = if pargs["umbrella-sampling"]
    (UmbrellaAverager,
      ((acc, chain) -> begin;
        UmbrellaAverager(
                         StandardAverager(
                                          numeric_type(0.0*acc(chain)),
                                          numeric_type(0.0),
                                          acc
                                         ),
                         wf
                        );
     end),
     ((acc, chain) -> begin;
        UmbrellaAverager(
                         StandardAverager(
                                          Vector{numeric_type}(0.0*acc(chain)),
                                          numeric_type(0.0),
                                          acc
                                         ),
                         wf
                        );
     end)
     );
  else
    (StandardAverager,
      ((acc, chain) -> begin;
             StandardAverager(
                              numeric_type,
                              acc, chain
                             )
     end),
     ((acc, chain) -> begin;
             StandardAverager(
                              numeric_type,
                              Vector{numeric_type},
                              acc, chain
                             );
     end)
     );
  end
  scalar_averagers = (Avg{numeric_type,numeric_type}[
                       (avgcons(chain -> dot(end_to_end(chain),
                                             end_to_end(chain)), chain)),
                       (avgcons(chain -> dot(chain_μ(chain), chain_μ(chain)),
                                chain)),
                       (avgcons(chain -> chain.U, chain)),
                       (avgcons(chain -> chain.U*chain.U, chain)),
                      ]);
  vector_averagers = (Avg{Vector{numeric_type},numeric_type}[
                       avgconss(end_to_end, chain),
                       avgconss(chain -> map(x -> x*x, end_to_end(chain)), chain),
                       avgconss(chain_μ, chain),
                       avgconss(chain -> map(x -> x*x, chain_μ(chain)), chain)
                      ]);
  outfile = open("$(pargs["prefix"])_trajectory.csv", "w");
  writedlm(outfile, ["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], ',');
  rollfile = open("$(pargs["prefix"])_rolling.csv", "w");
  writedlm(rollfile, ["step" "r1" "r2" "r3" "r1sq" "r2sq" "r3sq" "rsq" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psq" "U" "Usq"], ',');

  start = time();
  last_update = start;
  nacc = 0;
  nacc_total = 0;
  natt = 0;
  for init=1:pargs["num-inits"]

    if !force_ensemble_flag
      for i=1:10
        if move!(chain, 1, 0, 0, r0; frac_mv=1.0) # start with admissible chain
          break;
        end
      end
    end

    for step=1:nsteps
      idx = rand(1:n(chain));
      dϕ = rand(dϕ_dist);
      dθ = (((pargs["do-flips"] && rand(Bool)) ? π - 2*chain.θs[idx] : 0)
            + rand(dθ_dist));
      trial_chain = EAPChain(chain);
      successful = if force_ensemble_flag
        move!(trial_chain, idx, dϕ, dθ);
      else
        move!(trial_chain, idx, dϕ, dθ, r0; frac_mv=pargs["chain-frac-step"]);
      end
      if successful && acceptor(trial_chain, rand())
        chain = trial_chain;
        nacc += 1;
        nacc_total += 1;
      end
      natt += 1;

      if time() - last_update > pargs["update-freq"]
        @info "elapsed: $(time() - start)";
        @info "init:    $init / $(pargs["num-inits"])";
        @info "step:    $step / $nsteps";
        last_update = time();
      end

      if (
          pargs["step-adjust-scale"] != 1.0 &&
          step % pargs["steps-per-adjust"] == 0
         ) # adjust step size?
        acc_ratio = nacc / natt;
        if (acc_ratio > pargs["step-adjust-ub"] &&
            ϕstep != π && θstep != π/2)
          @info "acceptance ratio is high; increasing step size";
          nacc = 0;
          natt = 0;
          ϕstep = min(π, ϕstep*pargs["step-adjust-scale"]);
          θstep = min(π/2, θstep*pargs["step-adjust-scale"]);
        elseif acc_ratio < pargs["step-adjust-lb"]
          @info "acceptance ratio is low; decreasing step size";
          nacc = 0;
          natt = 0;
          ϕstep /= pargs["step-adjust-scale"];
          θstep /= pargs["step-adjust-scale"];
        end
        dϕ_dist = Uniform(-ϕstep, ϕstep);
        dθ_dist = Uniform(-θstep, θstep);
      end

      #for callback in callbacks
      #  callback(chain, step, false, false, pargs);
      #end
      foreach(a -> record!(a, chain), scalar_averagers);
      foreach(a -> record!(a, chain), vector_averagers);
      if step % pargs["stepout"] == 0
        writedlm(outfile,
                 hcat(step, transpose(chain.r),
                      transpose(chain_μ(chain)), chain.U),
                 ',');
        writedlm(rollfile,
                 hcat(
                      step,
                      transpose(get_avg(vector_averagers[1])), # r
                      transpose(get_avg(vector_averagers[2])), # rj2
                      get_avg(scalar_averagers[1]), # r2
                      transpose(get_avg(vector_averagers[3])), # p
                      transpose(get_avg(vector_averagers[4])), # pj2
                      get_avg(scalar_averagers[2]), # p2
                      get_avg(scalar_averagers[3]), # U
                      get_avg(scalar_averagers[4])  # U2
                     ),
                 ',');

      end

    end # steps

    # reinitialize polymer chain
    new_chain = EAPChain(pargs);
    new_chain.U = U(new_chain);
    if (
        pargs["force-init"] ||
        metropolis_acc(chain.kT, new_chain.U - chain.U,
                       prod(chain.sθs), prod(new_chain.sθs), rand())
       )
      chain = new_chain;
    end

  end

  ar = nacc_total / (pargs["num-inits"]*pargs["num-steps"]);
  @info "total time elapsed: $(time() - start)";
  @info "acceptance rate: $ar";

  #for callback in callbacks
  #  callback(chain, pargs["num-steps"], false, true, pargs);
  #end
  close(outfile);
  close(rollfile);

  return (scalar_averagers, vector_averagers, ar);
end

(sas, vas, ar) = if pargs["profile"]
  @info "Profiling the mcmc code...";
  mcmc(5, pargs); # run first to compile code
  @profview mcmc(pargs["num-steps"], pargs);
  println("Press ENTER to continue...");
  readline(stdin);
  exit(1);
else
  mcmc(pargs["num-steps"], pargs);
end
total_steps = pargs["num-steps"] * pargs["num-inits"];

println("<r>    =   $(get_avg(vas[1]))");
println("<r/nb> =   $(get_avg(vas[1]) / (pargs["mlen"]*pargs["num-monomers"]))");
println("<rj2>  =   $(get_avg(vas[2]))");
println("<r2>   =   $(get_avg(sas[1]))");
println("<p>    =   $(get_avg(vas[3]))");
println("<pj2>  =   $(get_avg(vas[4]))");
println("<p2>   =   $(get_avg(sas[2]))");
println("<U>    =   $(get_avg(sas[3]))");
println("<U2>   =   $(get_avg(sas[4]))");
println("AR     =   $ar");

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
using Test, Unmarked
using Random, StatsModels, StatsFuns, StatsBase, DataFrames, CategoricalArrays

include("test_data.jl")
include("test_design.jl")
include("test_links.jl")
include("test_formula.jl")
include("test_occu.jl")
include("test_nmix.jl")
include("test_categorical.jl")

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
# Python xlrd와 성능 비교

using EzExcel
using Benchmark


text:
module Example3

import MultipleShootingPEM
ms = MultipleShootingPEM
using LaTeXStrings
using ProgressMeter
using JLD2

############### Nonlinear system ###############
GL = 9.8/0.3
KA = 2.
M = 3.
D = 0.01
function nl_sys(x, u; gl=GL, ka=KA, m=M, δ=D)
    x1 = x[1] + δ*x[2]
    x2 = -δ*gl*sin(x[1]) + (1 - δ*ka/m)*x[2] + δ*1/m*u[1]
    return [x1, x2]
end

function derivative_x(x, u; gl=GL, ka=KA, m=M, δ=D)
    J = [1 δ;-δ*gl*cos(x[1]) (1 - δ*ka/m)]
    return J
end

function derivative_gl(x, u; gl=GL, ka=KA, m=M, δ=D)
    J = [0; -δ*sin(x[1])]
end

function derivative_ka(x, u; gl=GL, ka=KA, m=M, δ=D)
    J = [0; -δ*1/m*x[2]]
end

################ Controller ###############
GAIN = 4000
ZERO1 = 0.99
ZERO2 = 0.98
POLE1 = 0.02
CTRL = [GAIN, -GAIN*(ZERO1+ZERO2), GAIN*ZERO1*ZERO2, (1+POLE1), -POLE1]

################ Data Generator ##############
# Autonomous system
function pendulum(;N=2000, σv=0.,σw=0., seed=1, ampl=0, rep=20)
    srand(seed)
    if ampl != 0
        u = repeat(ampl*randn(Int(N//rep)), inner=rep)
        N = length(u)
    else
        u = zeros(N)
    end
    x0 = [π/2, 0]
    y = zeros(N)
    x = zeros(N+1, 2)
    x[1, :] = x0
    v = σv*randn(N)
    w = σw*randn(N)
    for k in 1:N
        x[k+1, :] = nl_sys(x[k, :], u[k]+ w[k])
        y[k] = x[k+1, 1] + v[k]
    end
    return Dict("u" => u, "y" => y, "x" => x, "v"=> v, "w" => w)
end

function inverted_pendulum(;N=2000, σv=0.,σw=0., ampl=0, rep=20, seed=1, transient=10000)
    srand(seed)
    Δr = repeat(ampl*randn(Int((N+transient)//rep)), inner=rep)
    total_len = length(Δr)
    n = length(Δr)
    u = 10*ones(total_len)
    x0 = [π, 0]
    x = zeros(total_len, 2)
    x[4, :] = x0
    y = π * ones(total_len)
    r = π * ones(total_len) + Δr
    e = zeros(total_len)
    v = σv*randn(total_len)
    w = σw*randn(total_len)
    for k in 4:total_len-1
        x[k+1, :] = nl_sys(x[k, :], u[k]+ w[k])
        y[k] = x[k+1, 1] + v[k]
        e[k] = r[k] - y[k]
        u[k+1] = dot(CTRL, [e[k], e[k-1], e[k-2], u[k], u[k-1]])
    end
    u, y, x = u[transient:end], y[transient:end], x[transient:end, :]
    return Dict("u" => u, "y" => y, "x" => x, "v" => v,
                "w" => w, "e" => e, "r" => r)
end

function data_generator(; t="pendulum", kwargs...)
    if t == "pendulum"
        return pendulum(; kwargs...)
    elseif t == "inverted_pendulum"
        return inverted_pendulum(; kwargs...)
    end
end

################ Output error model ##############
function output_error_model(data_dict; sim_len=1)
    u = data_dict["u"]
    y = data_dict["y"]
    xx = data_dict["x"]
    k0 = 2
    function g(x_next, dx, dθ, x, k, θ)
        gl = θ[1]
        ka = θ[2]
        k = k + k0
        x_next .= nl_sys(x, u[k]; gl=gl, ka=ka)
        dx .= derivative_x(x, u; gl=gl, ka=ka)
        dθ[:, 1] .= derivative_gl(x, u[k]; gl=gl, ka=ka)
        dθ[:, 2] .= derivative_ka(x, u[k]; gl=gl, ka=ka)
        return
    end

    function h(z, dx, dθ, x, k, θ)
        z .= x[1]
        dx .= [1 0]
        dθ .= [0 0]
        return
    end

    N = length(y)
    k0_list = collect(1:sim_len:N-k0-1)
    x0_ideal_list = [xx[k, :] for k in k0_list+k0]
    x0_list = [[y[k-1], (y[k] - y[k-1])/D] for k in k0_list+k0]
    y_list = [[yi] for yi in y[k0:end-1]]
    θ = [0, 0]
    return (ms.MultipleShooting(g, h, x0_ideal_list, y_list, k0_list, θ),
            ms.OptimizationProblem(g, h, x0_list, y_list, k0_list, θ))
end

function pem(data_dict; model="output_error", kwargs...)
    if model=="output_error"
        return output_error_model(data_dict; kwargs...)
    end
end

################ Compute cost function on grid of points ##############
function grid_cost_funtion(multiple_shoot; gl=(10, 60), ka=(0.5, 10), npoints=(50, 50))
    n_experiments = npoints[1]*npoints[2]
    if n_experiments == 0
        return Dict()
    end
    g_l_range = linspace(gl[1],  gl[2], npoints[1])
    ka_range = linspace(ka[1], ka[2], npoints[2])
    cost = Matrix{Float64}(npoints...)
    i = 1
    p = Progress(n_experiments, 1)
    for g_l in g_l_range
        j = 1
        for ka in ka_range
            x0_list = [oss.x0 for oss in multiple_shoot.simulations]
            ms.new_simulation!(multiple_shoot, x0_list, [g_l, ka])
            cost[i, j] = ms.cost_function(multiple_shoot)
            next!(p)
            j += 1
        end
        i += 1
    end

    return Dict("gl" => g_l_range, "ka" => ka_range, "cost" => cost')
end

################ Solve optimization problem for grid of initial points ##############
function solve_grid(opt; gl=(10, 60), ka=(0.5, 10), npoints=(3, 3))
    n_experiments = npoints[1]*npoints[2]
    if n_experiments == 0
        return Dict()
    end
    g_l_range = linspace(gl[1],  gl[2], npoints[1])
    ka_range = linspace(ka[1], ka[2], npoints[2])
    θ_list = Array{Float64}[]
    cost_list = Float64[]
    res_list = []
    MSE_list = []
    θ0_list = [[g_l, ka] for g_l in g_l_range for ka in ka_range]
    p = Progress(n_experiments, 1)
    for θ0 in θ0_list
        x0_list = [oss.x0 for oss in opt.ms.simulations]
        θ_ext = zeros(2 + length(x0_list)*2)
        ms.build_extended_vector!(θ_ext,  θ0, x0_list, 2, length(x0_list), 2)
        ms.update_simulation!(opt, θ_ext)
        res = ms.solve(opt, options=Dict("gtol" => 1e-10,
                                         "xtol" => 1e-10,
                                         "maxiter" => 2000,
                                         "initial_tr_radius" => 0.01,
                                         "verbose" => 2))
        θ_est = res["x"][1:2]
        cost_est = res["fun"]
        res = delete!(res, "jac")

        push!(θ_list, θ_est)
        push!(cost_list, cost_est)
        push!(res_list, res)
        next!(p)
    end
    return Dict("theta"=> θ_list, "cost"=> cost_list, "info" => res_list)
end

############ Runner ##############
function runner(options_dict)
    println("Generating data...")
    data_dict = data_generator(;options_dict["data_generator"]...)
    println("Done!")

    println("Creating prediction model...")
    multiple_shoot, opt = pem(data_dict; options_dict["pem"]...)
    println("Done!")

    println("Computing cost function in a grid...")
    grid_cost_fun = grid_cost_funtion(multiple_shoot; options_dict["grid_cost_funtion"]...)
    println("Done!")

    println("Solving for diferent initial conditions")
    solutions =  solve_grid(opt; options_dict["solve_grid"]...)
    println("Done!")
    Dict("data" => data_dict, "options"=>options_dict,
         "grid" => grid_cost_fun, "solutions" => solutions)
end

end # module


text:
function MLJBase.fit(model::EvoTreeRegressor, verbosity::Int, X, y)
    Xmatrix = MLJBase.matrix(X)
    fitresult = grow_gbtree(Xmatrix, y, model, verbosity = verbosity)
    cache = (Xmatrix, deepcopy(model))
    report = nothing
    return fitresult, cache, report
end

function MLJBase.update(model::EvoTreeRegressor, verbosity,
                        old_fitresult, old_cache, X, y)

    Xmatrix, old_model = old_cache
    δnrounds = model.nrounds - old_model.nrounds

    # We only continue computation from where we left off if: (i) The
    # number of iterations has not decreased; and (ii) All other
    # hyperparameters are unchanged, with the exception of the
    # learning rate (to allow for externally controlled adaptive
    # learning rate).

    okay_to_continue =
        δnrounds >= 0 &&
        model.loss == old_model.loss &&
        model.λ == old_model.λ &&
        model.γ == old_model.γ &&
        model.max_depth  == old_model.max_depth &&
        model.min_weight == old_model.min_weight &&
        model.rowsample ==  old_model.rowsample &&
        model.colsample ==  old_model.colsample &&
        model.nbins ==  old_model.nbins &&
        model.α ==  old_model.α &&
        model.metric ==  old_model.metric

    if okay_to_continue
        # Grow_gbtree! gets nrounds from old_fitresult.params, which is
        # indentical with model (and which `update` should not modify
        # permanently). For updating, we temporarily change this value
        # to δnrounds:
        old_nrounds = model.nrounds
        model.nrounds = δnrounds

        fitresult = grow_gbtree!(old_fitresult, Xmatrix, y;
                                 verbosity=verbosity)

        # now restore model:
        model.nrounds = old_nrounds
    else
        fitresult = grow_gbtree(Xmatrix, y, model, verbosity = verbosity)
    end

    cache = (Xmatrix, deepcopy(model))
    report = nothing

    return fitresult, cache, report
end

function MLJBase.predict(model::EvoTreeRegressor, fitresult, Xnew)
    Xmatrix = MLJBase.matrix(Xnew)
    pred = predict(fitresult, Xmatrix)
    return pred
end

# shared metadata
const EvoTypes = Union{EvoTreeRegressor}
MLJBase.input_scitype(::Type{<:EvoTreeRegressor}) = MLJBase.Table(MLJBase.Continuous)
# MLJBase.target_scitype(::Type{<:EvoTreeRegressor}) = MLJBase.Continuous

MLJBase.load_path(::Type{<:EvoTreeRegressor}) = "EvoTrees.MLJ.EvoTreeRegressor"
MLJBase.package_name(::Type{<:EvoTypes}) = "EvoTrees"
MLJBase.package_uuid(::Type{<:EvoTypes}) = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
MLJBase.package_url(::Type{<:EvoTypes}) = "https://github.com/Evovest/EvoTrees.jl"
MLJBase.is_pure_julia(::Type{<:EvoTypes}) = true

# function MLJ.clean!(model::EvoTreeRegressor)
#     warning = ""
#     if model.nrounds < 1
#         warning *= "Need nrounds ≥ 1. Resetting nrounds=1. "
#         model.nrounds = 1
#     end
#     if model.λ < 0
#         warning *= "Need λ ≥ 0. Resetting λ=0. "
#         model.λ = 0.0
#     end
#     if model.γ < 0
#         warning *= "Need γ ≥ 0. Resetting γ=0. "
#         model.γ = 0.0
#     end
#     if model.η <= 0
#         warning *= "Need η > 0. Resetting η=0.001. "
#         model.η = 0.001
#     end
#     if model.max_depth < 1
#         warning *= "Need max_depth ≥ 0. Resetting max_depth=0. "
#         model.max_depth = 1
#     end
#     if model.min_weight < 0
#         warning *= "Need min_weight ≥ 0. Resetting min_weight=0. "
#         model.min_weight = 0.0
#     end
#     if model.rowsample < 0
#         warning *= "Need rowsample ≥ 0. Resetting rowsample=0. "
#         model.rowsample = 0.0
#     end
#     if model.rowsample > 1
#         warning *= "Need rowsample <= 1. Resetting rowsample=1. "
#         model.rowsample = 1.0
#     end
#     if model.colsample < 0
#         warning *= "Need colsample ≥ 0. Resetting colsample=0. "
#         model.colsample = 0.0
#     end
#     if model.colsample > 1
#         warning *= "Need colsample <= 1. Resetting colsample=1. "
#         model.colsample = 1.0
#     end
#     if model.nbins > 250
#         warning *= "Need nbins <= 250. Resetting nbins=250. "
#         model.nbins = 250
#     end
#     return warning
# end


text:
export add_apply
export add_plus_c, add_minus_c, add_times_c, add_divide_c, add_diff_c
export add_threshold_c, add_setnz_c
export add_minimum_c, add_maximum_c, add_one_zero_maximum_c
export add_or_c, add_nand_c, add_nor_c, add_xor_c, add_xnor_c
export add_log_c

function add_apply(mgr::Ptr{Manager}, op::Ref, f::Ptr{Node}, g::Ptr{Node})
    res = ccall((:Cudd_addApply, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Cvoid}, Ptr{Node}, Ptr{Node}), mgr, op[], f, g)
    if res == C_NULL # Fails to compute a result
        throw(OutOfMemoryError())
    end
    return res
end

add_apply_func_names = [(:add_plus, :Cudd_addPlus),
                        (:add_minus, :Cudd_addMinus),
                        (:add_times, :Cudd_addTimes),
                        (:add_divide, :Cudd_addDivide),
                        (:add_diff, :Cudd_addDiff),
                        (:add_threshold, :Cudd_addThreshold),
                        (:add_setNZ, :Cudd_addSetNZ),
                        (:add_minimum, :Cudd_addMinimum),
                        (:add_maximum, :Cudd_addMaximum),
                        (:add_one_zero_maximum, :Cudd_addOneZeroMaximum),
                        (:add_agreement, :Cudd_addAgreement),
                        (:add_or, :Cudd_addOr),
                        (:add_nand, :Cudd_addNand),
                        (:add_nor, :Cudd_addNor),
                        (:add_xor, :Cudd_addXor),
                        (:add_xnor, :Cudd_addXnor),
                        (:add_log, :Cudd_addLog)]

for (name, Cudd_name) in add_apply_func_names
    @eval $name(mgr::Ptr{Manager}, f_pp::Ptr{Ptr{Node}}, g_pp::Ptr{Ptr{Node}}) =
                ccall(($(QuoteNode(Cudd_name)), libcudd),
                      Ptr{Node}, (Ptr{Manager}, Ptr{Ptr{Node}}, Ptr{Ptr{Node}}), mgr, f_pp, g_pp)

    # these Refs should get filled in __init__()
    @eval const $(Symbol(name, "_c")) = Ref{Ptr{Cvoid}}()
end

function __init__()
    for (name, _) in add_apply_func_names
        @eval $(Symbol(name, "_c"))[] = Compat.@cfunction($name, Ptr{Node}, (Ptr{Manager}, Ptr{Ptr{Node}}, Ptr{Ptr{Node}}))
    end
end


text:
d = 2
data = BanditDataset(d)
X = rand(2, 2)
a = [1, 2]
r = rand(1, 2)
append_data!(data, X, a, r)
@test data.X == X
@test data.a == a
@test data.r == r
@test arm_data(data, 1) == (X[:, 1:1], r[:, 1:1])
@test arm_data(data, 2) == (X[:, 2:2], r[:, 2:2])
@test_throws(ArgumentError("num_actions must be positive"), RandomPolicy(-1))
@test_throws(ArgumentError("d must be positive"), BanditDataset(-1))
@test_throws(ArgumentError("a must be positive"), arm_data(data, -1))
@test_throws(
    DimensionMismatch("X does not match the dimension of the dataset"),
    append_data!(data, X[1:1, :], a, r)
)


text:
import Base: size, getindex, setindex!, similar, ndims
#' in the future: Base.show(io::IO, tenmat::Tenmat) = print(io, "...")

abstract AbstractTenmat{T} # <: AbstractArray
##### functions to be implemented to make it sub of AbstractArray
# AbstractArray{T,N} where T is the type of object returned by integer indexing (A[1,...,1], when A is not empty)
# 	and N should be the length of the tuple returned by size()

"""
	Tenmat{T}(matrix::Array{T,2}, rowindex::Vector{Int}, colindex::Vector{Int}, tensorsize::Tuple{Vararg{Int}})

Matrix which stores the indices to convert it back to tensor.

## Example

```jldoctest
matrix = Array(reshape(1:48, 6, 8)) # some 6 x 8 matrix
rowindex = [3,1]
colindex = [2,4]
tensorsize = (2,4,3,2)
tenmat = Tenmat(matrix, rowindex, colindex, tensorsize)

# output

TensorMatrices_lemon.Tenmat{Int64}([1 7 … 37 43; 2 8 … 38 44; … ; 5 11 … 41 47; 6 12 … 42 48],[3,1],[2,4],(2,4,3,2))
```
"""
type Tenmat{T} <: AbstractTenmat{T}
	# e.g. Tenmat(A,[1,2],[3,4],(4,3,2,4)), where A is 12x8 matrix
	"converted matrix from tensor"
	matrix::Array{T,2}
	"Tensor indices to map to row"
	rowindex::Vector{Int}
	"Tensor indices to map to coloumn"
	colindex::Vector{Int}
	tensorsize::Tuple{Vararg{Int}} # tensor size

	function Tenmat{T}(matrix::Array{T,2}, rowindex::Vector{Int}, colindex::Vector{Int}, tensorsize::Tuple{Vararg{Int}})
		# the condition for the arguments:
		# the dimension of the mat matches to the tensize
		# -> the product of tensize[rowindex] == mat's row dimension
		# -> also for the column
		(length(rowindex) + length(colindex) == length(tensorsize)) &&
		(productWithIndex(tensorsize, rowindex) == size(matrix, 1)) &&
		(productWithIndex(tensorsize, colindex) == size(matrix, 2)) ?
 		new(matrix, rowindex, colindex, tensorsize)   :
		error("dimension does not match")
	end
end

Tenmat{T}(matrix::Matrix{T}, rowindex::Vector{Int}, colindex::Vector{Int}, tensorsize::Tuple{Vararg{Int}})  = Tenmat{T}(matrix, rowindex, colindex, tensorsize)

"""
	Tenmat{T}(matrix::Matrix{T})

Set a matrix as a Tenmat.
"""
Tenmat{T}(matrix::Matrix{T}) = Tenmat{T}(matrix, [1], [2], size(matrix))

"""
	Tenmat{T,N}(tensor::Array{T,N}, rowindex::Vector{Int}, colindex::Vector{Int})

Set a tensor as a Tenmat.

`Tenmat.matrix` will return the converted matrix.

## Example

```jldoctest basicMethods
A = Array(reshape(1:3024, 4,7,2,9,6)) # tensor, size of [4,7,2,9,6]
rowindex = [5,2] # 5th and 2nd indices maps to the row
colindex = [3,4,1] # 3rd, 4th and 1st indices maps to the column
tenmat = Tenmat(A,rowindex,colindex)

# output

TensorMatrices_lemon.Tenmat{Int64}([1 29 … 452 480; 505 533 … 956 984; … ; 2041 2069 … 2492 2520; 2545 2573 … 2996 3024],[5,2],[3,4,1],(4,7,2,9,6))
```
"""
Tenmat{T,N}(tensor::Array{T,N}, rowindex::Vector{Int}, colindex::Vector{Int}) = tensor2tenmat(tensor, rowindex, colindex)

##### define functions with macro
# tenmat_with_other_args_non_preserving = [:size, :getindex, : setindex!

"""
	size(tenmat::Tenmat, args...)

Return the size of the tenmat.matrix.

```jldoctest basicMethods
julia> size(tenmat)
(42,72)

julia> size(tenmat,2)
72
```
"""
size(tenmat::Tenmat, args...) = size(tenmat.matrix, args...)

"""
	setindex!(tenmat::Tenmat, args...)

Expands `setindex!` by tenmat.matrix.

```jldoctest basicMethods
julia> setindex!(tenmat, 100, 2, 6);
```
"""
setindex!(tenmat::Tenmat, args...) = setindex!(tenmat.matrix, args...)

"""
	getindex(tenmat::Tenmat, args...)

Expands `getindex` by tenmat.matrix.

```jldoctest basicMethods
julia> getindex(tenmat, 2, 6)
100
```
"""
getindex(tenmat::Tenmat, args...) = getindex(tenmat.matrix, args...)

"""
	similar

similar

## Examples

```jldoctest
A = rand(4,7,6)
tenmat = Tenmat(A,[3,1],[2])
similar(tenmat, Int)
similar(tenmat, Int, (5,2,7))
similar(tenmat, Int, [1,3],[2])
similar(tenmat, Int, [4,2,5], [1,3], (2,4,3,1,6))
similar(tenmat)
similar(tenmat, (5,2,7))
similar(tenmat, [1,3],[2])
B = similar(tenmat, [4,2,5], [1,3], (2,4,3,1,6))
size(B)

# output

(24,6)

```
"""
similar{T}(tenmat::Tenmat{T}) = Tenmat{T}(similar(tenmat.matrix), tenmat.rowindex, tenmat.colindex, tenmat.tensorsize)
similar{T}(tenmat::Tenmat{T}, S::Type) =
	similar(tenmat, S, tensorSize(tenmat))
similar{N}(tenmat::Tenmat, T::Type, tensorsize::Tuple{Vararg{Int64,N}}) =
	similar(tenmat, T, tenmat.rowindex, tenmat.colindex, tensorsize)
similar{T}(tenmat::Tenmat{T}, S::Type, rowindex::Vector{Int}, colindex::Vector{Int}) =
	similar(tenmat, S, rowindex, colindex, tensorSize(tenmat))
similar{T}(tenmat::Tenmat{T}, S::Type, rowindex::Vector{Int}, colindex::Vector{Int}, tensorsize::Tuple{Vararg{Int}}) =
	Tenmat(similar(tenmat2tensor(tenmat), T, tensorsize), rowindex, colindex)

similar{N}(tenmat::Tenmat, tensorsize::Tuple{Vararg{Int64,N}}) =
	similar(tenmat, tenmat.rowindex, tenmat.colindex, tensorsize)
similar{T}(tenmat::Tenmat{T}, rowindex::Vector{Int}, colindex::Vector{Int}) = similar(tenmat, rowindex, colindex, tensorSize(tenmat))
similar{T}(tenmat::Tenmat{T}, rowindex::Vector{Int}, colindex::Vector{Int}, tensorsize::Tuple{Vararg{Int}}) =
	Tenmat(similar(tenmat2tensor(tenmat), tensorsize), rowindex, colindex)

"""
	ndims(tenmat::Tenmat)

Return the original dimension of tenmat,
from `tensorsize`.

```jldoctest tenmat2
A = Array(reshape(1:3024, 4,7,2,9,6))
tenmat = Tenmat(A,[5,2],[3,4,1])
ndims(tenmat)

# output

5
```
"""
ndims(tenmat::Tenmat) = length(tensorSize(tenmat))

"""
	tensorSize(tenmat::Tenmat)

Return the tensorsize.

```jldoctest tenmat2
julia> tensorSize(tenmat)
(4,7,2,9,6)
```
"""
function tensorSize(tenmat::Tenmat)
	return tenmat.tensorsize
end

"""
	getMatrix(tenmat::Tenmat)

Return the matrix.

```jldoctest
A = Array(reshape(1:18, 3, 2, 3))
tenmat = Tenmat(A,[3,1],[2])
getMatrix(tenmat)

# output

9×2 Array{Int64,2}:
  1   4
  7  10
 13  16
  2   5
  8  11
 14  17
  3   6
  9  12
 15  18
```
"""
function getMatrix{T}(tenmat::Tenmat{T})
	return tenmat.matrix
end

"""
	productWithIndex{T}(tuple::Tuple{Vararg{T}}, index = collect(1:length(tuple))::Vector{Int})

Get a product of tuple elements with index number in index vector
"""
function productWithIndex{T}(tuple::Tuple{Vararg{T}}, index = collect(1:length(tuple))::Vector{Int})
	prod = 1
	for i in index
		prod *= tuple[i]
	end
	return prod
end

"""
	tensor2tenmat{T,N}(tensor::Array{T,N}, rowindex::Vector{Int}, colindex::Vector{Int})

Convert tensor to tenmat.

## examples

```jldoctest
julia> A = rand(3,5,3,2,6);

julia> tenmatA = tensor2tenmat(A,[3,5,2],[4,1]);

julia> size(getMatrix(tenmatA))
(90,6)

julia> tensorSize(tenmatA)
(3,5,3,2,6)

julia> B = Array(reshape(1:10,2,5));

julia> tenmatB = tensor2tenmat(B,[2],[1]);

julia> tenmatB.matrix == transpose(B)
true

julia> tensorSize(tenmatB)
(2,5)

```
"""
function tensor2tenmat{T,N}(tensor::Array{T,N}, rowindex::Vector{Int}, colindex::Vector{Int})
#=
julia> a = 5; b = 7; c = 8;
julia> D = reshape(1:a*b*c,a,b,c)
julia> tenmatD = tensor2tenmat(D,[3],[1,2])
julia> tenmatD.tensorsize == (a,b,c)
julia> size(tenmatD.matrix) == (c,a*b)
julia> for i = 1:a, j = 1:b, k = 1:c
julia> 	println(D[i,j,k] == tenmatD.matrix[k,(j-1)*a+i])
julia> end
julia>
julia> a = 5; b = 7; c = 8; d = 11
julia> D = reshape(1:a*b*c*d,a,b,c,d)
julia> tenmatD = tensor2tenmat(D,[3],[1,4,2])
julia> tenmatD.tensorsize == (a,b,c,d)
julia> size(tenmatD.matrix) == (c,a*d*b)
julia> for i = 1:a, j = 1:b, k = 1:c, l = 1:d
julia> 	println(D[i,j,k,l] == tenmatD.matrix[k,(j-1)*a*d+(l-1)*a+i])
julia> end

=#
	tensorsize = size(tensor)
	(N == length(rowindex) + length(colindex)) || error("index does not match")

	permten = permutedims(tensor,vcat(rowindex,colindex))
	rowdim = productWithIndex(size(tensor),rowindex)
	coldim = productWithIndex(size(tensor),colindex)

	matrix = reshape(permten,rowdim,coldim)
	tenmat = Tenmat(matrix,rowindex,colindex,tensorsize)

	return tenmat
end

"""
	tenmat2tensor{T}(tenmat::Tenmat{T})

Convert tenmat back to the tensor.

## examples
```jldoctest
julia> A = rand(3,5,3,2,8);

julia> tenmatA = tensor2tenmat(A,[3,5,2],[4,1]);

julia> back2A = tenmat2tensor(tenmatA);

julia> A == back2A
true
```

```jldoctest
B = rand(8,6,7,5,8,5,3,8)
tenmatB = tensor2tenmat(B,[1,7,3,2,6],[5,8,4])
back2B = tenmat2tensor(tenmatB)
check = true
for a=1:8, b=1:6, c=1:7, d=1:5, e=1:8, f=1:5, g=1:3, h=1:8
	back2B[a,b,c,d,e,f,g,h] == B[a,b,c,d,e,f,g,h] || check * false
end
check

# output

true
```
"""
function tenmat2tensor{T}(tenmat::Tenmat{T})
	indices = vcat(tenmat.rowindex, tenmat.colindex)
	tensor = reshape(tenmat.matrix, [tenmat.tensorsize[i] for i in indices]...)
	tensor = ipermutedims(tensor, indices)
	return tensor
end


text:
using Statistics
using Makie

function ppc_dens_overlay(y, yrep; kwargs...)
    f = Figure()
    Axis(f[1, 1])
    Makie.density!(y, color=:transparent, strokecolor = :black, strokewidth=2)

    for yr in yrep
        Makie.density!(yr, color=:transparent, strokecolor=(:blue, 0.1), strokewidth=1)
    end

    return f
end

function ppc_rootogram(y, yrep; prob=0.9, kwargs...)
    alpha = (1 - prob) / 2
    probs = [alpha, 1 - alpha]
    ymax = maximum([y; flatten(yrep)])

    f = Figure()
    Axis(f[1, 1])
    hist!(y)

    return f
end

function ppc_stat(y, yrep; kwargs...)
    fun = get(kwargs, :stat, mean)
    fig, ax, po = hist(fun.(yrep))
    vlines!(ax, [fun(y)], color=:black, linewidth=3)

    return fig
end

function ppc_intervals(y, yrep; kwargs...)
end

function ppc_stat_grouped(y, yrep; kwargs...)
end

max_stars_repo_name: kskyten/BayesPlot

max_issues_repo_name: kskyten/BayesPlot

max_forks_repo_name: kskyten/BayesPlot


text:
@test_throws DimensionMismatch splitobs((X, rand(149)), at=0.7)

@testset "typestability" begin
    @testset "Int" begin
        @test typeof(@inferred(splitobs(10, at=0.))) <: NTuple{2}
        @test eltype(@inferred(splitobs(10, at=0.))) <: UnitRange
        @test typeof(@inferred(splitobs(10, at=1.))) <: NTuple{2}
        @test eltype(@inferred(splitobs(10, at=1.))) <: UnitRange
        @test typeof(@inferred(splitobs(10, at=0.7))) <: NTuple{2}
        @test eltype(@inferred(splitobs(10, at=0.7))) <: UnitRange
        @test typeof(@inferred(splitobs(10, at=0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(10, at=(0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(10, at=0.5))) <: UnitRange
        @test eltype(@inferred(splitobs(10, at=(0.5,0.2)))) <: UnitRange
        @test eltype(@inferred(splitobs(10, at=0.5))) <: UnitRange
        @test eltype(@inferred(splitobs(10, at=(0.,0.2)))) <: UnitRange
    end
    for var in vars
        @test typeof(@inferred(splitobs(var, at=0.7))) <: NTuple{2}
        @test eltype(@inferred(splitobs(var, at=0.7))) <: SubArray
        @test typeof(@inferred(splitobs(var, at=0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(var, at=(0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(var, at=0.5))) <: SubArray
        @test eltype(@inferred(splitobs(var, at=(0.5,0.2)))) <: SubArray
    end
    for tup in tuples
        @test typeof(@inferred(splitobs(tup, at=0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(tup, at=(0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(tup, at=0.5))) <: Tuple
        @test eltype(@inferred(splitobs(tup, at=(0.5,0.2)))) <: Tuple
    end
end

@testset "Int" begin
    @test splitobs(10, at=0.7) == (1:7,8:10)
    @test splitobs(10, at=0.5) == (1:5,6:10)
    @test splitobs(10, at=(0.5,0.3)) == (1:5,6:8,9:10)
    @test splitobs(150, at=0.7) == splitobs(150, at=0.7)
    @test splitobs(150, at=0.5) == splitobs(150, at=0.5)
    @test splitobs(150, at=(0.5,0.2)) == splitobs(150, at=(0.5,0.2))
    @test numobs.(splitobs(150, at=0.7)) == (105,45)
    @test numobs.(splitobs(150, at=(.2,.3))) == (30,45,75)
    @test numobs.(splitobs(150, at=(.1,.2,.3))) == (15,30,45,60)
    # tests if all obs are still present and none duplicated
    @test sum(sum.(getobs.(splitobs(150, at=0.7)))) == 11325
    @test sum(sum.(splitobs(150,at=.1))) == 11325
    @test sum(sum.(splitobs(150,at=(.2,.1)))) == 11325
    @test sum(sum.(splitobs(150,at=(.1,.4,.2)))) == 11325
    @test sum.(splitobs(150, at=0.7)) == (5565, 5760)
end

@testset "Array, SparseArray, and SubArray" begin
    for var in (Xs, ys, vars...)
        @test numobs.(splitobs(var, at=10)) == (10,5)
        @test numobs.(splitobs(var, at=0.7)) == (10, 5)
        @test numobs.(splitobs(var, at=(.2,.3))) == (3,4, 8)
        @test numobs.(splitobs(var, at=(.11,.2,.3))) == (2,3,4,6)
        @test numobs.(splitobs(var, at=(2,3,4))) == (2,3,4,6)
    end
    # tests if all obs are still present and none duplicated
    @test sum(sum.(splitobs(X1,at=.1))) == 1200
    @test sum(sum.(splitobs(X1,at=(.2,.1)))) == 1200
    @test sum(sum.(splitobs(X1,at=(.1,.4,.2)))) == 1200
end

@testset "Tuple of Array, SparseArray, and SubArray" begin
    for tup in ((Xs,ys), (X,ys), (Xs,y), (Xs,Xs), (XX,X,ys), (X,yv), (Xv,y), tuples...)
        @test_throws MethodError splitobs(tup..., at=0.5)
        @test all(map(x->(typeof(x)<:Tuple), splitobs(tup,at=0.5)))
        @test numobs.(splitobs(tup, at=0.6)) == (9, 6)
        @test numobs.(splitobs(tup, at=(.2,.3))) == (3,4,8)
    end
end

# @testset "ObsView" begin
#     A = ObsView(X)
#     b, c = @inferred splitobs(A, .7)
#     @test b isa ObsView
#     @test c isa ObsView
#     @test b == A[1:105]
#     @test c == A[106:end]
# end

# @testset "BatchView" begin
#     A = BatchView(X, 5)
#     b, c = @inferred splitobs(A, .6)
#     @test b isa BatchView
#     @test c isa BatchView
#     @test b == A[1:18]
#     @test c == A[19:end]
# end


text:
using MultilayerPerceptron

using MLDatasets: Iris
using Random: randperm

function iris_example()
    features = begin
        features = Iris.features()
        normalize!(features)
        features
    end
    (n_features, n_samples) = size(features)
    n_classes = 0
    labels = begin
        labels = Iris.labels()
        label_dict = Dict(label => i for (i, label) in enumerate(Set(labels)))
        n_classes = length(label_dict)
        map(l -> label_dict[l], labels)
    end
    oh_labels = one_hot(labels, n_classes)

    # Prepares training and test set
    n_test = n_samples ÷ 3
    n_train = n_samples - n_test
    ind_train, ind_test = begin
        ind = randperm(n_samples)
        sort(ind[1:n_train]), sort(ind[(n_train + 1):end])
    end
    @assert length(ind_test) == n_test

    train_x = features[:, ind_train]
    train_y = labels[ind_train]
    oh_train_y = oh_labels[:, ind_train]

    test_x = features[:, ind_test]
    test_y = labels[ind_test]
    oh_test_y = oh_labels[:, ind_test]

    # The model contains a single output layer, which is equivalent to logistic regression.
    output_layer = Layer((n_classes, n_features))

    mlp = MLP([output_layer])
    fit!(mlp, train_x, oh_train_y;
        loss_fn=logit_cross_entropy,
        learning_rate=0.2, batch_size=5, n_epochs=200)

    oh_pred_y = softmax(predict(mlp, test_x))
    pred_y = inv_one_hot(oh_pred_y)

    n_correct = sum(pred_y .== test_y)
    println("Accuracy = $(n_correct / n_test)")

    confusion_matrix = oh_test_y * oh_pred_y'
    println("Confusion matrix =")
    display(Float16.(confusion_matrix))

    # The accuracy is usually above 90% (it varies from 70% - ~100%).
end

iris_example()

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
#= This module is based on Numerical Recipes. =#
module orderedTableSearch
    @inline function locate(x, grid, btm=0, up=0)
        # result is same as searchsortedlast
        #grid must be monotonically increasing.

    if up == 0
        up = length(grid) + 1
    end

    N = length(grid)

    ascnd = (grid[end] > grid[1])

        # evaluation is out of range
    if x < grid[1]
        return 0
    elseif x > grid[end]
        return N
    end

    mid = 0

    @inbounds while up - btm > 1
        mid = div((up+btm),2)#floor(Int, (up + btm) / 2)
        if  ascnd == (x >= grid[mid])
            btm = mid
        else
            up = mid
        end
    end

    return btm
end

    @inline function hunt(x, grid, init_btm::Int)
    N = length(grid)
    ascnd = (grid[end] > grid[1])
    btm = init_btm
    up = 0

    if init_btm < 1 || btm > N
        return locate(x, grid)
    else
        inc = 1 # increment
        if (x > grid[btm]) == ascnd
            @inbounds while true
                up = btm + inc
                if up > N - 1
                    up = N
                    break
                elseif (x > grid[up]) == ascnd
                    btm = up
                    inc += inc
                else
                    break
                end
            end
        else
            up = btm
            @inbounds while true
                btm = up - inc
                if btm < 1
                    btm = 1
                    break
                elseif (x < grid[btm]) == ascnd
                    up = btm
                    inc += inc
                else
                    break
                end

            end
        end
    end
    return locate(x, grid, btm, up)
end
end


text:
#!/usr/bin/env julia

using OrderedCollections,JSON
filename = @__FILE__

max=3000

# Maximum possible values e.g. 2^i2m < max
i2m = floor(Int,log(2,max))
i3m = floor(Int,log(3,max))
i5m = floor(Int,log(5,max))

# i2 start from 1 because we need even numbers
vals = [2^i2*3^i3*5^i5 for i2 =1:i2m for i3 =0:i3m for i5=0:i5m]
vals = sort(filter(x -> x <= max, vals))

r = OrderedCollections.OrderedDict()
r[raw"$comment"] = "File created by $filename"
r["type"] = "integer"
r["enum"]  = vals

io = open("nlon_nlat.schema.json","w")
JSON.print(io, r, 4)
close(io)


text:
using DelimitedFiles
using LinearAlgebra
using Plots
gr()
include("JuliaAnalysisPlotter.jl")

const plotID = ["Matlab",
"Julia-1-BLAS-Thread", "Julia-4-BLAS-Threads","Julia-MKL",
"Julia-SIMD-1-BLAS-Thread","Julia-SIMD-4-BLAS-Threads","Julia-MKL-SIMD"]
#=
choose among:
"Matlab",
"Julia-1-BLAS-Thread", "Julia-4-BLAS-Threads","Julia-MKL",
"Julia-SIMD-1-BLAS-Thread","Julia-SIMD-4-BLAS-Threads","Julia-MKL-SIMD"
=#

JuliaAnalysisPlotter(plotID)

max_stars_repo_name: aminya/Julia-Matlab-Benchmark

max_issues_repo_name: aminya/Julia-Matlab-Benchmark

max_forks_repo_name: aminya/Julia-Matlab-Benchmark