text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-

"""Pretty-print tabular data."""

_none_type = Nothing
_int_type = Int
_float_type = Float64
_text_type = UTF16String
_binary_type = Uint8

type Line
    start
    hline
    sep
    stop
end

type DataRow
    start
    sep
    stop
end

list(xs) = xs

# A table structure is suppposed to be:
#
#     --- lineabove ---------
#         headerrow
#     --- linebelowheader ---
#         datarow
#     --- linebewteenrows ---
#     ... (more datarows) ...
#     --- linebewteenrows ---
#         last datarow
#     --- linebelow ---------
#
# TableFormat's line* elements can be
#
#   - either None, if the element is not used,
#   - or a Line tuple,
#   - or a function: [col_widths], [col_alignments] -> string.
#
# TableFormat's *row elements can be
#
#   - either None, if the element is not used,
#   - or a DataRow tuple,
#   - or a function: [cell_values], [col_widths], [col_alignments] -> string.
#
# padding (an integer) is the amount of white space around data values.
#
# with_header_hide:
#
#   - either None, to display all table elements unconditionally,
#   - or a list of elements not to be displayed if the table has column headers.
#

type TableFormat
	lineabove
	linebelowheader
	linebetweenrows
	linebelow
	headerrow
	datarow
	padding
	with_header_hide

	function TableFormat(;
		lineabove        = "",
		linebelowheader  = "",
		linebetweenrows  = "",
		linebelow        = "",
		headerrow        = "",
		datarow          = "",
		padding          = "",
		with_header_hide = "",
		)

		return new(lineabove       ,linebelowheader ,linebetweenrows ,linebelow       ,headerrow       ,datarow         ,padding         ,with_header_hide)
	end
end

function _pipe_segment_with_colons(align, colwidth)
    """Return a segment of a horizontal line with optional colons which
    indicate column's alignment (as in `pipe` output format)."""
    w = colwidth
    if align in ["right", "decimal"]
        return ("-" ^ (w - 1)) * ":"
    elseif align == "center"
        return ":" * ("-" ^ (w - 2)) * ":"
    elseif align == "left"
        return ":" * ("-" ^ (w - 1))
    else
        return "-" ^ w
	end
end

function _pipe_line_with_colons(colwidths, colaligns)
    """Return a horizontal line with optional colons to indicate column's
    alignment (as in `pipe` output format)."""
    segments = [_pipe_segment_with_colons(a, w) for (a, w) in zip(colaligns, colwidths)]
    return "|" * "|".join(segments) * "|"
end

function mediawiki_row_with_attrs(separator, cell_values, colwidths, colaligns)
    alignment = { "left" =>    "",
                  "right" =>   "align=\"right\"| ",
                  "center" =>  "align=\"center\"| ",
                  "decimal" => "align=\"right\"| " }
    # hard-coded padding _around_ align attribute && value together
    # rather than padding parameter which affects only the value
    values_with_attrs = [(" " * get(alignment, a, "") * c * " ") for (c, a) in zip(cell_values, colaligns)]
    colsep = separator*2
    return (separator * rstrip(join(colsep,values_with_attrs)))
end

function _latex_line_begin_tabular(colwidths, colaligns)
    alignment = { "left" => "l", "right" => "r", "center" => "c", "decimal" => "r" }
    tabular_columns_fmt = join([alignment.get(a, "l") for (a) in colaligns], "")
    return "\\begin{tabular}{" * tabular_columns_fmt * "}\n\hline"
end

_table_formats = {"simple" =>
                  TableFormat(lineabove=Line("", "-", "  ", ""),
                              linebelowheader=Line("", "-", "  ", ""),
                              linebetweenrows=None,
                              linebelow=Line("", "-", "  ", ""),
                              headerrow=DataRow("", "  ", ""),
                              datarow=DataRow("", "  ", ""),
                              padding=0,
                              with_header_hide=["lineabove", "linebelow"]),
                  "plain" =>
                  TableFormat(lineabove=None, linebelowheader=None,
                              linebetweenrows=None, linebelow=None,
                              headerrow=DataRow("", "  ", ""),
                              datarow=DataRow("", "  ", ""),
                              padding=0, with_header_hide=None),
                  "grid" =>
                  TableFormat(lineabove=Line("+", "-", "+", "+"),
                              linebelowheader=Line("+", "=", "+", "+"),
                              linebetweenrows=Line("+", "-", "+", "+"),
                              linebelow=Line("+", "-", "+", "+"),
                              headerrow=DataRow("|", "|", "|"),
                              datarow=DataRow("|", "|", "|"),
                              padding=1, with_header_hide=None),
                  "pipe" =>
                  TableFormat(lineabove=_pipe_line_with_colons,
                              linebelowheader=_pipe_line_with_colons,
                              linebetweenrows=None,
                              linebelow=None,
                              headerrow=DataRow("|", "|", "|"),
                              datarow=DataRow("|", "|", "|"),
                              padding=1,
                              with_header_hide=["lineabove"]),
                  "orgtbl" =>
                  TableFormat(lineabove=None,
                              linebelowheader=Line("|", "-", "+", "|"),
                              linebetweenrows=None,
                              linebelow=None,
                              headerrow=DataRow("|", "|", "|"),
                              datarow=DataRow("|", "|", "|"),
                              padding=1, with_header_hide=None),
                  "rst" =>
                  TableFormat(lineabove=Line("", "=", "  ", ""),
                              linebelowheader=Line("", "=", "  ", ""),
                              linebetweenrows=None,
                              linebelow=Line("", "=", "  ", ""),
                              headerrow=DataRow("", "  ", ""),
                              datarow=DataRow("", "  ", ""),
                              padding=0, with_header_hide=None),
                  "mediawiki" =>
                  TableFormat(lineabove=Line("{| class=\"wikitable\" style=\"text-align: left;\"",
                                             "", "", "\n|+ <!-- caption -->\n|-"),
                              linebelowheader=Line("|-", "", "", ""),
                              linebetweenrows=Line("|-", "", "", ""),
                              linebelow=Line("|}", "", "", ""),
                              headerrow= (y...) -> _mediawiki_row_with_attrs("!", y...),
                              datarow= (y...) -> _mediawiki_row_with_attrs("|", y...),
                              padding=0, with_header_hide=None),
                  "latex" =>
                  TableFormat(lineabove=_latex_line_begin_tabular,
                              linebelowheader=Line("\\hline", "", "", ""),
                              linebetweenrows=None,
                              linebelow=Line("\\hline\n\\end{tabular}", "", "", ""),
                              headerrow=DataRow("", "&", "\\\\"),
                              datarow=DataRow("", "&", "\\\\"),
                              padding=1, with_header_hide=None),
                  "tsv" =>
                  TableFormat(lineabove=None, linebelowheader=None,
                              linebetweenrows=None, linebelow=None,
                              headerrow=DataRow("", "\t", ""),
                              datarow=DataRow("", "\t", ""),
                              padding=0, with_header_hide=None)}

tabulate_formats = list(sort(keys(_table_formats)))

_invisible_codes = r"\x1b\[\d*m"  # ANSI color codes
_invisible_codes_bytes = r"\x1b\[\d*m"  # ANSI color codes

function imple_separated_format(separator)
    """Construct a simple TableFormat with columns separated by a separator.

    >>> tsv = simple_separated_format("\\t") ; \
        tabulate([["foo", 1], ["spam", 23]], tablefmt=tsv) == "foo \\t 1\\nspam\\t23"
    true

    """
    return TableFormat(None, None, None, None,
                       headerrow=DataRow("", separator, ""),
                       datarow=DataRow("", separator, ""),
                       padding=0, with_header_hide=None)
end

function isconvertible(conv, string)
    try
        n = conv(string)
        return true
    catch ValueError
        return false
	end
end

function isnumber(string)
    """
    >>> _isnumber("123.45")
    true
    >>> _isnumber("123")
    true
    >>> _isnumber("spam")
    false
    """
    return _isconvertible(float, string)
end

function isint(string)
    """
    >>> _isint("123")
    true
    >>> _isint("123.45")
    false
    """
    return isa(string, _int_type) || (isa(string, _binary_type) || isa(string, _text_type)) && _isconvertible(int, string)
end

function gettype(string; has_invisible=true)

	if (has_invisible == true) && (isa(string, _text_type) || isa(string, _binary_type))
		string = _strip_invisible(string)
	end

	if string === Nothing
		return _none_type
	elseif hasattr(string, "isoformat")  # datetime.datetime, date, && time
        return _text_type
    elseif _isint(string)
        return int
    elseif _isnumber(string)
        return float
    elseif isa(string, _binary_type)
        return _binary_type
    else
        return _text_type
	end
end

function afterpoint(string)
    if _isnumber(string)
        if _isint(string)
            return -1
        else
            pos = rfind(string,".")
            pos =  pos < 1 ? rfind(lower(string), "e") : pos
			return pos >= 1 ? length(string) - pos - 1 :  -1  # no point
			end
    else
        return -1  # not a number
	end
end

function padleft(width, s, has_invisible=true)
    """Flush right.

    >>> _padleft(6, "\u044f\u0439\u0446\u0430") == "  \u044f\u0439\u0446\u0430"
    true

    """
    iwidth = has_invisible ? width * len(s) - len(_strip_invisible(s)) : width
    fmt = "{0:>%ds}" % iwidth
    return fmt.format(s)
end

function padright(width, s, has_invisible=true)
    """Flush left.

    >>> _padright(6, "\u044f\u0439\u0446\u0430") == "\u044f\u0439\u0446\u0430  "
    true

    """
    iwidth = has_invisible ? width * len(s) - len(_strip_invisible(s)) : width
    fmt = "{0:<%ds}" % iwidth
    return fmt.format(s)
end

function padboth(width, s, has_invisible=true)
    """Center string.

    >>> _padboth(6, "\u044f\u0439\u0446\u0430") == " \u044f\u0439\u0446\u0430 "
    true

    """
    iwidth = has_invisible ? width * len(s) - len(_strip_invisible(s)) : width
    fmt = "{0:^%ds}" % iwidth
    return fmt.format(s)
end

function strip_invisible(s)
    "Remove invisible ANSI color codes."
    if isa(s, _text_type)
        return re.sub(_invisible_codes, "", s)
    else   # a bytestring
        return re.sub(_invisible_codes_bytes, "", s)
	end
end

function visible_width(s)
    """Visible width of a printed string. ANSI color codes are removed.

    >>> _visible_width("\x1b[31mhello\x1b[0m"), _visible_width("world")
    (5, 5)

    """
    if isa(s, _text_type) || isa(s, _binary_type)
        return len(_strip_invisible(s))
    else
        return len(_text_type(s))
	end
end

function align_column(strings, alignment, minwidth=0, has_invisible=true)
    """[string] -> [padded_string]

    >>> list(map(str,_align_column(["12.345", "-1234.5", "1.23", "1234.5", "1e+234", "1.0e234"], "decimal")))
    ["   12.345  ", "-1234.5    ", "    1.23   ", " 1234.5    ", "    1e+234 ", "    1.0e234"]

    >>> list(map(str,_align_column(["123.4", "56.7890"], None)))
    ["123.4", "56.7890"]

    """
    if alignment == "right"
        strings = [s.strip() for (s) in strings]
        padfn = _padleft
    elseif alignment == "center"
        strings = [s.strip() for (s) in strings]
        padfn = _padboth
    elseif alignment == "decimal"
        decimals = [_afterpoint(s) for (s) in strings]
        maxdecimals = max(decimals)
        strings = [s * (maxdecimals - decs) * " "
                   for (s, decs) in zip(strings, decimals)]
        padfn = _padleft
    elseif not alignment
        return strings
    else
        strings = [s.strip() for (s) in strings]
        padfn = _padright
	end

    if has_invisible
        width_fn = _visible_width
    else
        width_fn = len
	end

    maxwidth = max(max(map(width_fn, strings)), minwidth)
    padded_strings = [padfn(maxwidth, s, has_invisible) for (s) in strings]
    return padded_strings
end

function more_generic(type1, type2)
    types = { _none_type: 0, int: 1, float: 2, _binary_type: 3, _text_type: 4 }
    invtypes = { 4: _text_type, 3: _binary_type, 2: float, 1: int, 0: _none_type }
    moregeneric = max(types.get(type1, 4), types.get(type2, 4))
    return invtypes[moregeneric]
end

function column_type(strings, has_invisible=true)
    """The least generic type all column values are convertible to.

    >>> _column_type(["1", "2"]) is _int_type
    true
    >>> _column_type(["1", "2.3"]) is _float_type
    true
    >>> _column_type(["1", "2.3", "four"]) is _text_type
    true
    >>> _column_type(["four", "\u043f\u044f\u0442\u044c"]) is _text_type
    true
    >>> _column_type([None, "brux"]) is _text_type
    true
    >>> _column_type([1, 2, None]) is _int_type
    true
    >>> import datetime as dt
    >>> _column_type([dt.datetime(1991,2,19), dt.time(17,35)]) is _text_type
    true

    """
    types = [_type(s, has_invisible) for (s) in strings ]
    return reduce(_more_generic, types, int)
end

function format(val, valtype, floatfmt, missingval="")
    """Format a value accoding to its type.

    Unicode is supported:

    >>> hrow = ["\u0431\u0443\u043a\u0432\u0430", "\u0446\u0438\u0444\u0440\u0430"] ; \
        tbl = [["\u0430\u0437", 2], ["\u0431\u0443\u043a\u0438", 4]] ; \
        good_result = "\\u0431\\u0443\\u043a\\u0432\\u0430      \\u0446\\u0438\\u0444\\u0440\\u0430\\n-------  -------\\n\\u0430\\u0437             2\\n\\u0431\\u0443\\u043a\\u0438           4" ; \
        tabulate(tbl, headers=hrow) == good_result
    true

    """
    if val === Nothing
        return missingval
	end

    if valtype === _int_type || valtype === _text_type
        return "{0}".format(val)
    elseif valtype === _binary_type
        try
            return _text_type(val, "ascii")
        catch TypeError
            return _text_type(val)
		end
    elseif valtype === _float_type
        return format(float(val), floatfmt)
    else
        return "{0}".format(val)
	end
end

function align_header(header, alignment, width)
    if alignment == "left"
        return _padright(width, header)
    elseif alignment == "center"
        return _padboth(width, header)
    elseif not alignment
        return "{0}".format(header)
    else
        return _padleft(width, header)
	end
end

function normalize_tabular_data(tabular_data, headers)
    """Transform a supported data type to a list of lists, && a list of headers.

    Supported tabular data types:

    * list-of-lists or another iterable of iterables

    * list of named tuples (usually used with headers="keys")

    * list of dicts (usually used with headers="keys")

    * list of OrderedDicts (usually used with headers="keys")

    * 2D NumPy arrays

    * NumPy record arrays (usually used with headers="keys")

    * dict of iterables (usually used with headers="keys")

    * p&&as.DataFrame (usually used with headers="keys")

    The first row can be used as headers if headers="firstrow",
    column indices can be used as headers if headers="keys".

    """

    if hasattr(tabular_data, "keys") && hasattr(tabular_data, "values")
        # dict-like && p&&as.DataFrame?
        if hasattr(tabular_data.values, "__call__")
            # likely a conventional dict
            keys = tabular_data.keys()
            rows = list(izip_longest(tabular_data.values()...))  # columns have to be transposed
        elseif hasattr(tabular_data, "index")
            # values is a property, has .index => it's likely a p&&as.DataFrame (p&&as 0.11.0)
            keys = tabular_data.keys()
            vals = tabular_data.values  # values matrix doesn't need to be transposed
            names = tabular_data.index
            rows = [[v]*list(row) for (v,row) in zip(names, vals)]
        else
            except("tabular data doesn't appear to be a dict or a DataFrame")
		end

        if headers == "keys"
            headers = list(map(_text_type,keys))  # headers should be strings
		end
    else   # it's a usual an iterable of iterables, or a NumPy array
        rows = list(tabular_data)

		if (headers == "keys" &&
            hasattr(tabular_data, "dtype") &&
            getattr(tabular_data.dtype, "names"))
            # numpy record array
            headers = tabular_data.dtype.names
		elseif (headers == "keys"
              && len(rows) > 0
              && isa(rows[0], tuple)
              && hasattr(rows[0], "_fields"))
            # namedtuple
            headers = list(map(_text_type, rows[0]._fields))
        elseif (len(rows) > 0 && isa(rows[0], dict))
            # dict or OrderedDict
            uniq_keys = set() # implements hashed lookup
            keys = [] # storage for set
            if headers == "firstrow"
                firstdict = len(rows) > 0 ? rows[0] : {}
                keys.extend(firstdict.keys())
                uniq_keys.update(keys)
                rows = rows[1:]
			end

            for (row) in rows
                for (k) in row.keys()
                    #Save unique items in input order
                    if k not in uniq_keys
                        keys.append(k)
                        uniq_keys.add(k)
					end
				end
			end

            if headers == "keys"
                headers = keys
            elseif headers == "firstrow" && len(rows) > 0
                headers = [firstdict.get(k, k) for (k) in keys]
                headers = list(map(_text_type, headers))
			end

			rows = [[row.get(k) for (k) in keys] for (row) in rows]

        elseif headers == "keys" && len(rows) > 0
            # keys are column indices
            headers = list(map(_text_type, range(len(rows[0]))))
		end
    # take headers from the first row if necessary
    if headers == "firstrow" && len(rows) > 0
        headers = list(map(_text_type, rows[0])) # headers should be strings
        rows = rows[1:]
	end

    headers = list(map(_text_type,headers))
    rows = list(map(list,rows))

    # pad with empty headers for initial columns if necessary
    if headers && len(rows) > 0
       nhs = len(headers)
       ncols = len(rows[0])
       if nhs < ncols
           headers = [""]^(ncols - nhs) * headers
	   end
   end
    return rows, headers
end

function abulate(tabular_data, headers=[], tablefmt="simple",
             floatfmt="g", numalign="decimal", stralign="left",
             missingval="")
    """Format a fixed width table for pretty printing.

    >>> print(tabulate([[1, 2.34], [-56, "8.999"], ["2", "10001"]]))
    ---  ---------
      1      2.34
    -56      8.999
      2  10001
    ---  ---------

    The first required argument (`tabular_data`) can be a
    list-of-lists (or another iterable of iterables), a list of named
    tuples, a dictionary of iterables, an iterable of dictionaries,
    a two-dimensional NumPy array, NumPy record array, or a P&&as'
    dataframe.

    Table headers
    -------------

    To print nice column headers, supply the second argument (`headers`):

      - `headers` can be an explicit list of column headers
      - if `headers="firstrow"`, then the first row of data is used
      - if `headers="keys"`, then dictionary keys or column indices are used

    Otherwise a headerless table is produced.

    If the number of headers is less than the number of columns, they
    are supposed to be names of the last columns. This is consistent
    with the plain-text format of R && P&&as' dataframes.

    >>> print(tabulate([["sex","age"],["Alice","F",24],["Bob","M",19]],
    ...       headers="firstrow"))
           sex      age
    -----  -----  -----
    Alice  F         24
    Bob    M         19

    Column alignment
    ----------------

    `tabulate` tries to detect column types automatically, && aligns
    the values properly. By default it aligns decimal points of the
    numbers (or flushes integer numbers to the right), && flushes
    everything else to the left. Possible column alignments
    (`numalign`, `stralign`) are: "right", "center", "left", "decimal"
    (only for `numalign`), && None (to disable alignment).

    Table formats
    -------------

    `floatfmt` is a format specification used for columns which
    contain numeric data with a decimal point.

    `None` values are replaced with a `missingval` string:

    >>> print(tabulate([["spam", 1, None],
    ...                 ["eggs", 42, 3.14],
    ...                 ["other", None, 2.7]], missingval="?"))
    -----  --  ----
    spam    1  ?
    eggs   42  3.14
    other   ?  2.7
    -----  --  ----

    Various plain-text table formats (`tablefmt`) are supported:
    "plain", "simple", "grid", "pipe", "orgtbl", "rst", "mediawiki",
    && "latex". Variable `tabulate_formats` contains the list of
    currently supported formats.

    "plain" format doesn't use any pseudographics to draw tables,
    it separates columns with a double space:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                 ["strings", "numbers"], "plain"))
    strings      numbers
    spam         41.9999
    eggs        451

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="plain"))
    spam   41.9999
    eggs  451

    "simple" format is like P&&oc simple_tables:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                 ["strings", "numbers"], "simple"))
    strings      numbers
    ---------  ---------
    spam         41.9999
    eggs        451

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="simple"))
    ----  --------
    spam   41.9999
    eggs  451
    ----  --------

    "grid" is similar to tables produced by Emacs table.el package or
    P&&oc grid_tables:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                ["strings", "numbers"], "grid"))
    +-----------+-----------+
    | strings   |   numbers |
    +===========+===========+
    | spam      |   41.9999 |
    +-----------+-----------+
    | eggs      |  451      |
    +-----------+-----------+

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="grid"))
    +------+----------+
    | spam |  41.9999 |
    +------+----------+
    | eggs | 451      |
    +------+----------+

    "pipe" is like tables in PHP Markdown Extra extension or P&&oc
    pipe_tables:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                ["strings", "numbers"], "pipe"))
    | strings   |   numbers |
    |:----------|----------:|
    | spam      |   41.9999 |
    | eggs      |  451      |

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="pipe"))
    |:-----|---------:|
    | spam |  41.9999 |
    | eggs | 451      |

    "orgtbl" is like tables in Emacs org-mode && orgtbl-mode. They
    are slightly different from "pipe" format by not using colons to
    define column alignment, && using a "+" sign to indicate line
    intersections:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                ["strings", "numbers"], "orgtbl"))
    | strings   |   numbers |
    |-----------+-----------|
    | spam      |   41.9999 |
    | eggs      |  451      |

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="orgtbl"))
    | spam |  41.9999 |
    | eggs | 451      |

    "rst" is like a simple table format from reStructuredText; please
    note that reStructuredText accepts also "grid" tables:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]],
    ...                ["strings", "numbers"], "rst"))
    =========  =========
    strings      numbers
    =========  =========
    spam         41.9999
    eggs        451
    =========  =========

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="rst"))
    ====  ========
    spam   41.9999
    eggs  451
    ====  ========

    "mediawiki" produces a table markup used in Wikipedia && on other
    MediaWiki-based sites:

    >>> print(tabulate([["strings", "numbers"], ["spam", 41.9999], ["eggs", "451.0"]],
    ...                headers="firstrow", tablefmt="mediawiki"))
    {| class="wikitable" style="text-align: left;"
    |+ <!-- caption -->
    |-
    ! strings   !! align="right"|   numbers
    |-
    | spam      || align="right"|   41.9999
    |-
    | eggs      || align="right"|  451
    |}

    "latex" produces a tabular environment of LaTeX document markup:

    >>> print(tabulate([["spam", 41.9999], ["eggs", "451.0"]], tablefmt="latex"))
    \\begin{tabular}{lr}
    \\hline
     spam &  41.9999 \\\\
     eggs & 451      \\\\
    \\hline
    \\end{tabular}

    """

    list_of_lists, headers = _normalize_tabular_data(tabular_data, headers)

    # optimization: look for ANSI control codes once,
    # enable smart width functions only if a control code is found
    plain_text = "\n".join(["\t".join(map(_text_type, headers))] + \
                            ["\t".join(map(_text_type, row)) for (row) in list_of_lists])
    has_invisible = re.search(_invisible_codes, plain_text)
    if has_invisible
        width_fn = _visible_width
    else
        width_fn = len
	end
    # format rows && columns, convert numeric values to strings
    cols = list(zip(list_of_lists...))
    coltypes = list(map(_column_type, cols))
    cols = [[_format(v, ct, floatfmt, missingval) for (v) in c]
             for (c,ct) in zip(cols, coltypes)]

    # align columns
    aligns = [ ct in [int,float] ? numalign : stralign for (ct) in coltypes]
    minwidths = headers ? [width_fn(h)+2 for (h) in headers] : [0]*len(cols)
    cols = [_align_column(c, a, minw, has_invisible)
            for (c, a, minw) in zip(cols, aligns, minwidths)]

    if headers
        # align headers && add headers
        minwidths = [max(minw, width_fn(c[0])) for (minw, c) in zip(minwidths, cols)]
        headers = [_align_header(h, a, minw) for (h, a, minw) in zip(headers, aligns, minwidths)]
        rows = list(zip(cols...))
    else
        minwidths = [width_fn(c[0]) for (c) in cols]
        rows = list(zip(cols...))
	end

    if not isa(tablefmt, TableFormat)
        tablefmt = _table_formats.get(tablefmt, _table_formats["simple"])
	end

    return _format_table(tablefmt, headers, rows, minwidths, aligns)
end

function build_simple_row(padded_cells, rowfmt)
    "Format row according to DataRow format without padding."
    (beginl, sep, endl) = rowfmt
    return rstrip(beginl * join(padded_cells, sep) * endl)
end

function build_row(padded_cells, colwidths, colaligns, rowfmt)
    "Return a string which represents a row of data cells."
    if not rowfmt
        return None
    if hasattr(rowfmt, "__call__")
        return rowfmt(padded_cells, colwidths, colaligns)
    else
        return _build_simple_row(padded_cells, rowfmt)
	end
end

function build_line(colwidths, colaligns, linefmt)
    "Return a string which represents a horizontal line."
    if not linefmt
        return None
	end

    if hasattr(linefmt, "__call__")
        return linefmt(colwidths, colaligns)
    else
        beginl, fill, sep, endl = linefmt
        cells = [fill*w for (w) in colwidths]
        return _build_simple_row(cells, (beginl, sep, endl))
	end
end

function pad_row(cells, padding)
    if cells
        pad = " "*padding
        padded_cells = [pad * cell * pad for (cell) in cells]
        return padded_cells
    else
        return cells
	end
end

function format_table(fmt, headers, rows, colwidths, colaligns)
    """Produce a plain-text representation of the table."""
    lines = []
    hidden = (headers && fmt.with_header_hide) ? fmt.with_header_hide : []
    pad = fmt.padding
    headerrow = fmt.headerrow

    padded_widths = [(w * 2^pad) for (w) in colwidths]
    padded_headers = _pad_row(headers, pad)
    padded_rows = [_pad_row(row, pad) for (row) in rows]

    if fmt.lineabove && "lineabove" not in hidden
        lines.append(_build_line(padded_widths, colaligns, fmt.lineabove))
	end

    if padded_headers
        lines.append(_build_row(padded_headers, padded_widths, colaligns, headerrow))
        if fmt.linebelowheader && "linebelowheader" not in hidden
            lines.append(_build_line(padded_widths, colaligns, fmt.linebelowheader))
		end
	end

    if padded_rows && fmt.linebetweenrows && "linebetweenrows" not in hidden
        # initial rows with a line below
        for (row) in padded_rows[:-1]
            lines.append(_build_row(row, padded_widths, colaligns, fmt.datarow))
            lines.append(_build_line(padded_widths, colaligns, fmt.linebetweenrows))
		end
        # the last row without a line below
        lines.append(_build_row(padded_rows[-1], padded_widths, colaligns, fmt.datarow))
    else
        for (row) in padded_rows
            lines.append(_build_row(row, padded_widths, colaligns, fmt.datarow))
		end
	end

    if fmt.linebelow && "linebelow" not in hidden
        lines.append(_build_line(padded_widths, colaligns, fmt.linebelow))
	end

    return "\n".join(lines)
end


text:
export extract_params, energy, delta_energy, dms_silico

######################################################
"""
    read_par_BM(path::AbstractString, q::Integer = 21)
    read_par_BM(path::GZipStream, q::Integer = 21)

	Reads the parameters of a Potts model in format

	J i j a b
	...
	h i a

	and returns them in tensor format.
	J[a, b, i, j] is such that  1 == "A", ... "21" == "-",
	same for h[a, i].

    "path": full path of the parameter file.
    "q": number of symbols allowed (21 for protein sequences, gaps included).

"""

#-----------------------------------------------

function read_par_BM(path::AbstractString, q::Integer = 21)
    params = readdlm(path,' ', use_mmap = true)[:, 2:6]
    l_file = size(params, 1)
    N = Integer(((q - 2) + sqrt( (q-2)^2 + 8*l_file))/(2*q))
	J = Array{Float64}(undef, q, q, N, N)
	h = Array{Float64}(undef, q, N)
	n_J = Int(q*q*N*(N-1)/2)
	n_h = q*N

	for k in 1:n_J
		i, j, a, b, par_j = params[k, :]
		i += 1
		j += 1
		a == 0 && (a = 21)
		b == 0 && (b = 21)
		J[a, b, i, j] = par_j
	end

	for l in (n_J + 1): n_h + n_J
		i, a, par_h = params[l, :]
		i += 1
		a == 0 && (a = 21)
		h[a, i] = par_h
	end

	return h, J
end

function read_par_BM(path::GZipStream, q::Integer = 21)
    params = readdlm(path,' ', use_mmap = true)[:, 2:6]
    l_file = size(params, 1)
    N = Integer(((q - 2) + sqrt( (q-2)^2 + 8*l_file))/(2*q))
	J = Array{Float64}(undef, q, q, N, N)
	h = Array{Float64}(undef, q, N)
	n_J = Int(q*q*N*(N-1)/2)
	n_h = q*N

	for k in 1:n_J
		i, j, a, b, par_j = params[k, :]
		i += 1
		j += 1
		a == 0 && (a = 21)
		b == 0 && (b = 21)
		J[a, b, i, j] = par_j
	end

	for l in (n_J + 1): n_h + n_J
		i, a, par_h = params[l, :]
		i += 1
		a == 0 && (a = 21)
		h[a, i] = par_h
	end

	return h, J
end

######################################################
"""
    function extract_params(path_par::AbstractString; q::Integer = 21)
    function extract_params(g::DCAgraph)

    (path_par, kwargs..) --> h, J

    USE:
    Given the path of bmDCA parameters,
    return fields "h" and coupling "J"

    INPUT:
    "path_par": path of DCA parameters

    "q": number of symbols in the MSA

    OUTPUT:

    h, J

"""

#-----------------------------------------------

function extract_params(path_par::AbstractString; q::Integer = 21)
	!isfile(path_par) && error("Error: the file \"$(path_par)\" does not exist. Please check the spelling or the folder path.")
	file = GZip.open(path_par)
	h, J = read_par_BM(file, q)
	h = set_max_field_to_0(h)
	J = symmetrize_J(J)
	return h, J
end

#This function needs DCATools to import redefinition of .* and other operators.
#function extract_params(g::DCAgraph)
#    N = g.L
 #   q = g.q
 #   J = Array{Float64}(undef, q, q, N, N)
 #   for i in 1:N
 #       for j in 1:N
 #           for a in 1:q
 #               aa = amino_plus(a)
 #               for b in 1:q
 #                   bb = amino_plus(b)
 #                   J[a, b, i, j ]= g.J[(i .- 1)*q .+ aa, (j .- 1)*q .+ bb]
 #               end
 #           end
 #       end
 #   end
 #   h = Array{Float64}(undef, q, N)
 #   for i in 1:N
 #       for a in 1:q
 #           aa = amino_plus(a)
 #           h[a, i] = g.h[(i .- 1)*g.q .+ aa]
 #   end
 #   end
 #   return SeqEvol.set_max_field_to_0(h), J
#end

######################################################
"""
    set_max_field_to_0(h_old::Array{Float64, 2})

    Takes in input a 2-dimensional tensor, a matrix,
    containing the values of the fiels for every amino acid a
    and position i.
    Returns the same matrix, by removing the biggerst field value
    for every position i, so that the biggest field results 0.
    Help in computing exponentials without going to overflow.

"""

#-----------------------------------------------

function set_max_field_to_0(h_old::Array{Float64, 2})
   h_new = deepcopy(h_old)
   q, N = size(h_old)
   for i in 1:N
        hmax = maximum([h_old[a, i] for a in 1:21])
        h_new[:, i] .-= hmax
   end
   return h_new
end

######################################################
"""
    symmetrize_J(J_old::Array{Float64, 4})

    Takes as input a 4-dimensional tensor,
    containing the values of the couplings for every couple of
    amino acids a, b and positions i, j.
    Returns the same tensor, by symmetrizing it,
    since J[a, b, i, j] = 0 if i < j, or viceversa, depends.

"""

#-----------------------------------------------

function symmetrize_J(J_old::Array{Float64, 4})
   J_s = deepcopy(J_old)
   q,q, N, N = size(J_old)
   for i in 1:N
        for j in i+1:N
            for a in 1:q
                for b in 1:q
                      J_s[b, a, j, i]  = J_old[a, b, i, j]
                end
            end
        end
    end
    return J_s
end

######################################################
"""
    energy(h::Array{Float64,2},
           J::Array{Float64,4},
           S::Array{<:Integer, 1})
    energy(h::Array{Float64,2},
           J::Array{Float64,4},
           S::Array{<:Integer, 2})

    (h, J, S) --> E

    USE:
    Given fields and couplings of a Potts models,
    returns its energy.

    INPUT:
    "h": fields in qxN format

    "J": coulings in qxqxNxN format

    "S": sequence in number format (" 1 2 .. 21 <--> A C -  ")
    "S": MSA in number format (" 1 2 .. 21 <--> A C -  ")

    OUTPUT:
    Energy
"""

#-----------------------------------------------

function energy(h::Array{Float64,2},
                J::Array{Float64,4},
                S::Array{<:Integer, 1})
    q, N = size(h)
    E = 0.0
    @fastmath for i = 1:N
         E -= h[S[i],i]
         @fastmath for j = (i+1):N
             E -= J[S[j],S[i],j,i]
        end
    end
    return E
end

function energy(h::Array{Float64,2},
                J::Array{Float64,4},
                S::Array{<:Integer, 2})
    q, N = size(h)
    M, N = size(S)

    E = zeros(M)
    for i in 1:M
      E[i] = energy(h, J, S[i, :])
    end

    return E
end

######################################################
"""
	delta_energy(h::Array{Float64,2},
                              J::Array{Float64,4},
                              S::Array{<:Integer,1},
                              ref::Array{<:Integer, 1})

    (h, J, S, ref) --> deltaE

    USE:
    Given fields and couplings of a Potts models, and two sequence
    returns the difference in energy.

    INPUT:
    "h": fields in qxN format

    "J": coulings in qxqxNxN format

    "S": sequence in number format (" 1 2 .. 21 <--> A C -  ")

    "ref": sequence in number format

    OUTPUT:

    Energy[S] - Energy[ref]
"""

#-----------------------------------------------

function delta_energy(h::Array{Float64,2},
                              J::Array{Float64,4},
                              S::Array{<:Integer,1},
                              ref::Array{<:Integer, 1})

    q, N = size(h)
    E = 0.0
    index_v = collect(1:N)
    bool_common = Array{Bool}(undef, N)
    ll = 0
    @fastmath for val in S
        ll+=1
        @inbounds bool_common[ll] = (val == ref[ll])
    end
    common = index_v[bool_common]
    non_common = index_v[(!).(bool_common)]

    @fastmath for i in non_common
         @inbounds  E -= (h[S[i],i] - h[ref[i],i])
         for j = 1:N
            if j > i
              @inbounds  E -= (J[S[j],S[i],j,i] - J[ref[j],ref[i],j,i] )
            end
        end
    end

    @fastmath for i in common
         for j in non_common
            if j > i
              @inbounds  E -= (J[S[j],S[i],j,i] - J[ref[j],ref[i],j,i] )
            end
        end
    end

    return E
end

######################################################
"""
    dms_silico(h::Array{Float64,2},
		J::Array{Float64,4},
		wt::Array{<:Integer,1};
		skipgap::Bool = false,
		skipsyn::Bool = false,
		no_df = false)

    (h, J, wt) --> in silico DMS

    USE:
    Given a sequence, computes the Deep Mutational Scanning in silico
    and returns a dataframe with the output.

    INPUT:
    "h": fields in qxN format
    "J": coulings in qxqxNxN format
    "wt": sequence in number format (" 1 2 .. 21 <--> A C -  ")

    optional:
	"skipgap": do not return gapped positions
	"skisyn":  do not return synonimus mutations (ΔE = 0)
	"no_df":   return only a vector

    OUTPUT:
    Dataframe in the following format

    :res    :wt_amino  :wt_var  :ΔE
    ________________________________
      x         x         x      x

"""

#-----------------------------------------------

function dms_silico(h::Array{Float64,2}, J::Array{Float64,4}, wt::Array{<:Integer,1}; skipgap::Bool = false, skipsyn::Bool = false, df = true)
    L = length(wt)

    dms = DataFrame(res = Int64.(zeros(L*20)), wt_amino = Int64.(zeros(L*20)),
        var_amino = Int64.(zeros(L*20)),  ΔE = zeros(L*20) )
    allowmissing!(dms)
    for res in 1:L
        for amino in 1:20
            seq = copy(wt)
            wt_amino = wt[res]
            seq[res] = amino
            ΔE = delta_energy(h, J, seq, wt)
            wt[res] == 21 && (ΔE = missing)
            dms[(res-1)*20 + amino, :] .= [res,wt_amino, amino, ΔE]
        end
    end
    dms = skipgap ? dropmissing!(dms) : dms
    dms = skipsyn ? dms[dms[:, :wt_amino] .!= dms[:, :var_amino], :] : dms
    dms = df ? dms : dms[:, :ΔE]
    return dms
end

ext: jl

max_stars_repo_name: matteobisardi/KitMSA

max_issues_repo_name: matteobisardi/KitMSA

max_forks_repo_name: matteobisardi/KitMSA


text:
module Caesar

# import RoME: initfg # collision on RoME.initfg() since no parameters are given in both RoME and Caesar
import Distributions: Normal
import RoME: getRangeKDEMax2D
import IncrementalInference: getSample, initfg
import DistributedFactorGraphs: getManifold

# handy project consts (not for export)
import IncrementalInference: NothingUnion, InstanceType

using Requires
using Dates

using Manifolds
using StaticArrays

import Rotations as _Rot

# TODO remove
const _Rotations = _Rot

using
  Pkg,
  DelimitedFiles,
  Distributed,
  Statistics,
  LinearAlgebra,
  IncrementalInference,
  TransformUtils,
  CoordinateTransformations,
  JSON,
  JSON2,
  FileIO,
  DataStructures,
  ProgressMeter,
  ImageMagick, # figure out who else is using this and move to requires
  ImageCore,
  DocStringExtensions,
  Unmarshal,
  YAML,
  FFTW,
  TimeZones,
  TensorCast

using Optim

using Reexport

# public API exports
include("ExportAPI.jl")

## ===============================================================================================
# and source files
include("BearingRangeTrackingServer.jl")

include("SlamServer.jl")
include("DataUtils.jl")
include("UserFunctions.jl")

# Configuration
include("config/CaesarConfig.jl")

include("Deprecated.jl")

# Multisession operation
# include("attic/multisession/Multisession.jl")

# SAS-SLAM
include("beamforming/czt.jl")
include("beamforming/CBF.jl")
include("beamforming/MatchedFilter.jl")
include("beamforming/SASBearing2D.jl")
include("beamforming/SASUtils.jl")

# conditional loading for ROS
function __init__()
  # ZMQ server and endpoints
  @require ZMQ="c2297ded-f4af-51ae-bb23-16f91089e4e1" include("zmq/ZmqCaesar.jl")
  @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    @info "Loading Caesar PyCall specific utilities (using PyCall)."
    @eval using .PyCall
    @require RobotOS="22415677-39a4-5241-a37a-00beabbbdae8" include("ros/Utils/RosbagSubscriber.jl")
  end
  @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" include("3rdParty/_PCL/_PCL.jl")
  @require AprilTags="f0fec3d5-a81e-5a6a-8c28-d2b34f3659de" begin
    include("images/apriltags.jl")
    @require ImageDraw="4381153b-2b60-58ae-a1ba-fd683676385f" include("images/AprilTagDrawingTools.jl")
  end
  @require ImageDraw="4381153b-2b60-58ae-a1ba-fd683676385f" include("images/imagedraw.jl")
  @require ImageMagick="6218d12a-5da1-5696-b52f-db25d2ecc6d1" include("images/imagedata.jl")
  @require Images="916415d5-f1e6-5110-898d-aaa5f9f070e0" begin
    include("images/images.jl")
    include("images/ScanMatcherUtils.jl")
    include("images/ScanMatcherPose2.jl")
    include("images/ScatterAlignPose2.jl")

    @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" include("plotting/ScatterAlignPlotting.jl")
  end
  @require Distributed="8ba89e20-285c-5b6f-9357-94700520ee1b" include("images/DistributedUtils.jl")
end

end


text:
import Luna
import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes, Fields
import Luna.PhysData: wlfreq
import Test: @test, @testset

@testset "multi-point" begin
Z = [0, 0.25, 0.5, 1]
P = [0, 1, 0.5, 0]
coren, densityfun = Capillary.gradient(:He, Z, P)
for ii = 1:3
    λ = Z[ii+1] - Z[ii]
    ζ = range(0, λ; length=10)
    z = range(Z[ii], Z[ii+1]; length=10)
    p0, p1 = P[ii], P[ii+1]
    Pζ = @. sqrt(p0^2 + ζ/λ*(p1^2 - p0^2))
    ρζ = PhysData.density.(:He, Pζ)
    ρz = densityfun.(z)
    @test all(isapprox.(ρζ, ρz, rtol=1e-10))
end
end

@testset "field" begin
a = 13e-6
gas = :Ar
pres = 5
τ = 30e-15
λ0 = 800e-9
L = 5e-2

# Common setup
grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 0.5e-12)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

# Constant
dens0 = PhysData.density(gas, pres)
dens(z) = dens0
m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid)
linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(
    grid, dens, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_const = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_const, status_period=10)

# Gradient
coren, densityfun = Capillary.gradient(gas, L, pres, pres)
m = Capillary.MarcatiliMode(a, coren, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(
    grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_grad = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_grad, status_period=10)

# Gradient array
coren, densityfun = Capillary.gradient(gas, [0,L], [pres, pres]);
m = Capillary.MarcatiliMode(a, coren, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(
    grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_grad_array = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_grad_array, status_period=10)

@test all(output_grad.data["Eω"][grid.sidx, :] .≈ output_const.data["Eω"][grid.sidx, :])
@test all(output_grad_array.data["Eω"][grid.sidx, :] .≈ output_const.data["Eω"][grid.sidx, :])
end

@testset "envelope" begin
a = 13e-6
gas = :Ar
pres = 5
τ = 30e-15
λ0 = 800e-9
L = 5e-2

# Common setup
grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 0.5e-12)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)

# Constant
dens0 = PhysData.density(gas, pres)
dens(z) = dens0
m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid);
linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0);
Eω, transform, FT = Luna.setup(grid, dens, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_const = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_const, status_period=10)

# Gradient
coren, densityfun = Capillary.gradient(gas, L, pres, pres)
m = Capillary.MarcatiliMode(a, coren, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(
    grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_grad = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_grad, status_period=10)

# Gradient array
coren, densityfun = Capillary.gradient(gas, [0,L], [pres, pres]);
m = Capillary.MarcatiliMode(a, coren, loss=false)
aeff(z) = Modes.Aeff(m, z=z)
energyfun, energyfunω = Fields.energyfuncs(grid)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(
    grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
output_grad_array = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output_grad_array, status_period=10)

@test all(output_grad.data["Eω"][grid.sidx, :] .≈ output_const.data["Eω"][grid.sidx, :])
@test all(output_grad_array.data["Eω"][grid.sidx, :] .≈ output_const.data["Eω"][grid.sidx, :])
end


text:
# This file is a part of SimilaritySearch.jl

export Performance, StatsKnn, StatsComparison, scores, probe

mutable struct DistCounter{DistType<:PreMetric} <: PreMetric
    dist::DistType
    count::Int
end

function evaluate(D::DistCounter, a, b)
    D.count += 1
    evaluate(D.dist, a, b)
end

"""
    StatsKnn(res::KnnResult)
    StatsKnn(res::AbstractVector)

Computes some basic statistics for the KnnResult or an array of KnnResults
"""
struct StatsKnn
    distancessum::Float64
    nearestdist::Float64
    farthestdist::Float64
    length::Float64

    StatsKnn(res::KnnResult) = new(sum(p[end] for p in res), minimum(res), maximum(res), length(res))
    function StatsKnn(reslist::AbstractVector{T}) where {T<:KnnResult}
        distancessum = 0.0
        nearestdist = 0.0
        farthestdist = 0.0
        len = 0.0

        for res in reslist
            p = StatsKnn(res)
            distancessum += p.distancessum
            nearestdist += p.nearestdist
            farthestdist += p.farthestdist
            len += p.length
        end

        n = length(reslist)
        new(distancessum/n, nearestdist/n, farthestdist/n, len/n)
    end
end

"""
    StatsComparison

Stores the result of the comparison between two indexes, see [`Performance`](@ref).
"""
struct StatsComparison
    macrorecall::Float64
    macroprecision::Float64
    macrof1::Float64
    searchtime::Float64
    evaluations::Float64
    stats::StatsKnn
    goldsearchtime::Float64
    goldstats::StatsKnn
end

"""
    Performance(_goldsearch::AbstractSearchContext, queries::AbstractVector, ksearch::Integer; popnearest=false)

Creates performance comparer for the given set of queries using a gold standard index.

- `_goldsearch`: a gold standard index
- `queries`: a set of queries
- `ksearch`: the number of neighbors to retrieve
- `popnearest`: set as `true` whenever queries are part of the dataset.
"""
struct Performance{DataType<:AbstractVector, KnnResultType<:KnnResult}
    queries::DataType
    ksearch::Int
    popnearest::Bool
    goldreslist::Vector{KnnResultType}
    goldsearchtime::Float64
    goldevaluations::Float64
    goldstats::StatsKnn
end

function Performance(_goldsearch::AbstractSearchContext, queries::AbstractVector, ksearch::Integer; popnearest=false)
    dist = DistCounter(_goldsearch.dist, 0)
    goldsearch = copy(_goldsearch, dist=dist)
    gold, searchtime, evaluations = perf_search_batch(goldsearch, queries, ksearch, popnearest)
    Performance(queries, Int(ksearch), popnearest, gold, searchtime, evaluations, StatsKnn(gold))
end

function perf_search_batch(index::AbstractSearchContext, queries, ksearch::Integer, popnearest::Bool)
    m = length(queries)
    if popnearest
        ksearch += 1
    end
    reslist = [KnnResult(ksearch) for i in 1:m]
    search(index, queries[1], ksearch) # warming step
    evaluations = index.dist.count
    start = time()

    for i in 1:m
        search(index, queries[i], reslist[i])
        popnearest && popfirst!(reslist[i])
    end

    elapsed = time() - start
    reslist, elapsed / m, (index.dist.count - evaluations) / m
end

"""
    probe(perf::Performance, _index::AbstractSearchContext)

Compares the performance of `_index` with the gold standard index of `perf`.
"""
function probe(perf::Performance, _index::AbstractSearchContext)
    index = copy(_index, dist=DistCounter(_index.dist, 0))
    reslist, searchtime, evaluations = perf_search_batch(index, perf.queries, perf.ksearch, perf.popnearest)
    n = length(reslist)
    recall = 0.0
    nearest = 0.0
    precision = 0.0
    f1 = 0.0
    for i in eachindex(reslist)
        p = scores(perf.goldreslist[i], reslist[i])
        recall += p.recall
        precision += p.precision
        f1 += p.f1
    end

    StatsComparison(recall/n, precision/n, f1/n, searchtime, evaluations, StatsKnn(reslist), perf.goldsearchtime, perf.goldstats)
end

"""
    scores(gold, res)

Compute recall and precision scores from the result sets.
"""
function scores(gold::Set, res::Set)
    tp = intersect(gold, res)  # tn
    fp = length(setdiff(res, tp)) # |fp| == |ft| when |res| == |gold|
    fn = length(setdiff(gold, tp))
    _tp = length(tp)
    recall = _tp / (_tp + fn)
    precision = _tp / (_tp + fp)

    (recall=recall, precision=precision, f1=2 * precision * recall / (precision + recall))
end

scores(gold::KnnResult, res::KnnResult) = scores(Set(gold.id), Set(res.id))


text:
const TEST_FLOAT32 = false
const TEST_DOUBLE64 = false
const TEST_STATIC_SIZED = false

using Manifolds
using ManifoldsBase
using ManifoldsBase: number_of_coordinates
import ManifoldsBase: active_traits, merge_traits

using LinearAlgebra
using Distributions
using DoubleFloats
using Quaternions
using Random
using StaticArrays
using Statistics
using StatsBase
using Test
using Graphs
using SimpleWeightedGraphs

function include_test(path)
    @info "Testing $path"
    @time include(path)  # show basic timing, (this will print a newline at end)
end

function our_ambiguities(m=Base)
    ambigs = Test.detect_ambiguities(m)
    modules_we_care_about =
        [Base, LinearAlgebra, Manifolds, ManifoldsBase, StaticArrays, Statistics, StatsBase]
    our_ambigs = filter(ambigs) do (m1, m2)
        we_care = m1.module in modules_we_care_about && m2.module in modules_we_care_about
        return we_care && (m1.module === Manifolds || m2.module === Manifolds)
    end
    return our_ambigs
end

"""
    has_type_in_signature(sig, T)
    Test whether the signature `sig` has an argument of type `T` as one of its paramaters
"""
function has_type_in_signature(sig, T::Type)
    return any(map(Base.unwrap_unionall(sig.sig).parameters) do x
        xw = Base.rewrap_unionall(x, sig.sig)
        return (xw isa Type ? xw : xw.T) <: T
    end)
end


text:
#####
##### Value-Priority pairs
#####

struct Priority{T}
    priority::Int
    val::T
end

Base.isless(a::P, b::P) where {P <: Priority} = isless(a.priority, b.priority)
Base.:(==)(a::P, b::P) where {P <: Priority} = a.priority == b.priority
Base.hash(a::Priority, h::UInt = UInt(0x78907890)) = hash(a.priority, h)
unwrap(x::Priority) = x.val

#####
##### LRU implementation
#####

const MutableMinHeap{T} = DataStructures.MutableBinaryHeap{T,Base.ForwardOrdering}
mutable struct LRU{T}
    heap::MutableMinHeap{Priority{T}}
    handles::Dict{T,Int}
end

LRU{T}() where {T} = LRU(MutableMinHeap{Priority{T}}(), Dict{T,Int}())
Base.isempty(lru::LRU) = isempty(lru.heap)
Base.length(lru::LRU) = length(lru.heap)

function Base.pop!(lru::LRU)
    v = pop!(lru.heap)
    delete!(lru.handles, v.val)
    return v.val
end

function Base.push!(lru::LRU{T}, v::T, priority::Int) where {T}
    # Assert this for now.
    @check !haskey(lru.handles, v)

    # Add this to the heap and record its handle.
    lru.handles[v] = push!(lru.heap, Priority(priority, v))
    return v
end

Base.first(lru::LRU) = unwrap(first(lru.heap))

function update!(lru::LRU{T}, v::T, priority::Int) where {T}
    handle = lru.handles[v]
    DataStructures.update!(lru.heap, handle, Priority(priority, v))
end

function Base.delete!(lru::LRU{T}, v::T) where {T}
    delete!(lru.heap, lru.handles[v])
    delete!(lru.handles, v)
    return lru
end

Base.in(v::T, lru::LRU{T}) where {T} = haskey(lru.handles, v)


text:
using CUSPARSE
using Base.Test

using CUSPARSE: cusparseMatDescr_t, cusparseCreateMatDescr
mdptr = cusparseMatDescr_t[0]
cusparseCreateMatDescr(mdptr)
md = mdptr[1]

using CUSPARSE: cusparseGetMatDiagType, cusparseGetMatFillMode, cusparseGetMatIndexBase, cusparseGetMatType
@show cusparseGetMatDiagType(md)
@show cusparseGetMatFillMode(md)
@show cusparseGetMatIndexBase(md)
@show cusparseGetMatType(md)


text:
function details(; kwargs...)
    resultjson = requestjson("place/details"; kwargs...)
    if resultjson["status"] == "OK"
        return resultjson["result"]
    else
        return []
    end
end

"""
Returns the results of a place details request.

Once you have a place_id or a reference from a Place Search, you can request
more details about a particular establishment or point of interest by initiating
a Place Details request. A Place Details request returns more comprehensive
information about the indicated place such as its complete address, phone
number, user rating and reviews.

https://developers.google.com/places/web-service/details

## Parameters
* `key::String`: Your Google Maps API Key.
* `placeid`: An identifier that uniquely identifies a place, returned from a
    Place Search. For more information about place IDs, see the place ID
    overview https://developers.google.com/places/web-service/place-id.

"""
details(key::String, placeid::String; kwargs...) =
    details(key=key, placeid=placeid; kwargs...)


text:
module SBTestDirection

using SpatialBoundaries
using Test

@test last(SpatialBoundaries._rategradient(1.0,1.0)) == 225.0
@test last(SpatialBoundaries._rategradient(-1.0,-1.0)) == 45.0
@test last(SpatialBoundaries._rategradient(-1.0,0.0)) == 90.0
@test last(SpatialBoundaries._rategradient(1.0,0.0)) == 270.0
@test last(SpatialBoundaries._rategradient(0.0,-1.0)) == 180.0
@test last(SpatialBoundaries._rategradient(0.0,1.0)) == 0.0
@test last(SpatialBoundaries._rategradient(1.0,-1.0)) == 315.0
@test last(SpatialBoundaries._rategradient(-1.0,1.0)) == 135.0

end


text:
abstract type AbstractJSObject end

"""
References objects stored in Javascript.
Maps the following expressions to actual calls on the Javascript side:
```julia
jso = JSObject(name, scope, typ)
# getfield:
x = jso.property # returns a new JSObject
# setfield
jso.property = "newval" # Works with JSObjects or Julia objects as newval
# call:
result = jso.func(args...) # works with julia objects and other JSObjects as args

# constructors are wrapped this way:
scene = jso.new.Scene() # same as in JS: scene = new jso.Scene()
```
"""
mutable struct JSObject <: AbstractJSObject
    # fields are private and not accessible via jsobject.field
    name::Symbol
    session::Session
    typ::Symbol
    # transporting the UUID allows us to have a uuid different from the objectid
    # which will help to better capture === equivalence on the js side.
    uuid::UInt64

    function JSObject(name::Symbol, scope::Session, typ::Symbol)
        obj = new(name, scope, typ)
        setfield!(obj, :uuid, objectid(obj))
        # finalizer(remove_js_reference, obj)
        return obj
    end

    function JSObject(name::Symbol, scope::Session, typ::Symbol, uuid::UInt64)
        obj = new(name, scope, typ, uuid)
        # finalizer(remove_js_reference, obj)
        return obj
    end
end

struct JSGlobal <: AbstractJSObject
    name::Symbol
end

macro jsglobal(name)
    JSGlobal(Symbol(name))
end

function serialize_string(io::IO, g::JSGlobal)
    print(io, g.name)
end

function serialize_string(io::JSONSerializer, jso::JSObject)
    serialize_string(io, Dict(
        :__javascript_type__ => :JSObject,
        :payload => uuidstr(jso)
    ))
end
function serialize_string(io::IO, jso::JSObject)
    serialize_string(io, js"get_heap_object($(uuidstr(jso)))")
end

"""
    JSObject(jso::JSObject, typ::Symbol)

Copy constructor with a new `typ`
"""
function JSObject(jso::JSObject, typ::Symbol)
    jsonew = JSObject(name(jso), session(jso), typ)
    # point new object to old one on the javascript side:
    evaljs(session(jso), js"put_on_heap($(uuidstr(jsonew)), $jso); undefined;")
    return jsonew
end
function JSObject(session::Session, name::Symbol)
    return JSObject(name, session, :variable)
end

# define accessors
for name in (:name, :session, :typ, :uuid)
    @eval $(name)(jso::JSObject) = getfield(jso, $(QuoteNode(name)))
end

"""
Removes an JSObject from the object pool!
"""
function remove_js_reference(jso::JSObject)
    evaljs(session(jso), js"delete $jso")
end

"""
    uuidstr(jso::JSObject)

Returns the uuid as a string
"""
uuidstr(jso::JSObject) = string(uuid(jso))

"""
Overloading getproperty to allow the same semantic as Javascript.
Since there is no `new` keyword in Julia like in JS, we missuse
jsobject.new, to return an instance of jsobject with a new modifier.

So this Javascript:
```js
obj = new Module.Constructor()
```

Will translates to the following Julia code:
```Julia
obj = Module.new.Constructor()
```
"""
function Base.getproperty(jso::AbstractJSObject, field::Symbol)
    if field === :new
        # Create a new instance of jso, with the `new` modifier
        return JSObject(jso, :new)
    else
        result = JSObject(field, session(jso), typ(jso))
        send(
            session(jso),
            type = JSGetIndex,
            object = jso,
            field = field,
            result = uuidstr(result),
        )
        return result
    end
end

function Base.setproperty!(jso::AbstractJSObject, field::Symbol, value)
    send(
        session(jso),
        type = JSSetIndex,
        object = jso,
        value = value,
        field = field
    )
    return value
end

"""
    construct_arguments(args, keyword_arguments)
Constructs the arguments for a JS call.
Can only use either keyword arguments or positional arguments.
"""
function construct_arguments(args, keyword_arguments)
    if isempty(keyword_arguments)
        return args
    elseif isempty(args)
        # tojs isn't recursive bug:
        return keyword_arguments
    else
        # TODO: I'm not actually sure about this :D
        error("""
        Javascript only supports keyword arguments OR arguments.
        Found posititional arguments and keyword arguments
        """)
    end
end

"""
    jsobject(args...; kw_args...)

Call overload for JSObjects.
Only supports keyword arguments OR positional arguments.
"""
function jscall(jso::AbstractJSObject, args, kw_args)
    result = JSObject(:result, session(jso), :call)
    send(
        session(jso),
        type = JSCall,
        func = jso,
        needs_new = getfield(jso, :typ) === :new,
        arguments = construct_arguments(args, kw_args),
        result = uuidstr(result)
    )
    return result
end

(jso::JSObject)(args...; kw_args...) = jscall(jso, args, kw_args)
(jso::JSGlobal)(args...; kw_args...) = jscall(jso, args, kw_args)

struct JSModule <: AbstractJSObject
    session::Session
    mod::JSObject
    document::JSObject
    window::JSObject
    this::JSObject
    display_func # A function that gets called on show with the modules Scope
end
session(x::JSModule) = getfield(x, :session)

function make_renderable!(jsm::JSModule)
    jss = js"""
        function (mod){
            $(object_pool_identifier) = {}
            $(object_pool_identifier)[$(uuidstr(jsm.mod))] = mod
            $(object_pool_identifier)[$(uuidstr(jsm.document))] = document
            $(object_pool_identifier)[$(uuidstr(jsm.window))] = window
            $(object_pool_identifier)[$(uuidstr(jsm.this))] = this
        }
    """
    onimport(jsm.scope, jss)
    return jsm.display_func(jsm.scope)
end

fuse(f, jso::JSObject) = fuse(f, session(jso))


text:
function skipto(str,char)
	for ci in 1:length(str)
		#print(str[ci])
		if str[ci]==char
			return ci
		end
	end
	return 0
end
function depvert(fname)
	text=readall("$fname.jl")
	deploc=search(text,"Union(")
	while !isempty(collect(deploc))
		depend=skipto(text[deploc[end]:end],')')
		#println(depend,text,deploc)
		text=text[1:deploc[end]-1]*"{"*text[deploc[end]+1:end]
		text=text[1:deploc[end]-2+depend]*"}"*text[deploc[end]+depend:end]
		deploc=search(text[deploc[end]+depend:end],"Union(")+deploc[end]+depend-1
	end
	f=open("$fname.jl","w")
	write(f,text)
	close(f)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 71a70a62-28fb-11eb-2092-7d3640b6ddf5
using Pkg, DrWatson

# ╔═╡ 71d2cdaa-28fb-11eb-3826-cf35b99cde1f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ d31428e6-85d2-11eb-396b-c9d6da91e9ae
md" ### Clip-10-01-04s.jl"

# ╔═╡ ed43435a-85d2-11eb-0aea-2ffa1acc2cc5
m = [0 0 10 0 0;
	0 1 8 1 0;
	0 2 6 2 0;
	1 2 4 2 1;
	2 2 2 2 2];

# ╔═╡ c65f5f04-8770-11eb-2407-05dbb86e82b7
begin
	pr_m = float.(deepcopy(m))
	r, c = size(pr_m)
	for row in 1:r
		pr_m[row, :] = [pr_m[row, i] / sum(pr_m[row, :]) for i in 1:c]
	end
	pr_m
end

# ╔═╡ 8b5eb0ee-881a-11eb-022c-adf0d7440e84
function factorial_fraction(m)
	w = []
	for i in 1:size(m, 1)
		n = m[i,:]
		append!(w, [factorial(sum(n))/prod(factorial, n)])
	end
	w
end

# ╔═╡ 2b2c55a8-8821-11eb-0176-91dca1254a8b
ways = factorial_fraction(m)

# ╔═╡ edac53be-881a-11eb-1dbd-a189c57ac336
begin
	N = 10
	logwayspp = (1/N) * log.(factorial_fraction(m))
end

# ╔═╡ 5480e430-876a-11eb-0c63-aff9c5cf17d2
begin
	h = zeros(c)
	for i in 1:r
		s = [pr_m[i, j] == 0 ? 0 : pr_m[i, j] * log(pr_m[i, j]) for j in 1:c]
		h[i] = -sum(s)
	end
	h
end

# ╔═╡ ae0c69d8-85d6-11eb-2db9-01a3f56c9a65
begin
	fig = plot(logwayspp, h, ylims=(-0.1, 1.75),
		xlab="log(ways) per pebble", ylab="Entropy", leg=false)
	scatter!(logwayspp, h)
	for (ind, sym) in enumerate([:A, :B, :C, :D, :E])
		annotate!([(logwayspp[ind], h[ind] + 0.1,
			Plots.text(sym, 10, :red, :right))])
	end
	fig
end

# ╔═╡ 22d43054-8661-11eb-0dd7-fb846edb81c4
md" ### End clip-10-01-04s.jl"

# ╔═╡ Cell order:
# ╟─d31428e6-85d2-11eb-396b-c9d6da91e9ae
# ╠═71a70a62-28fb-11eb-2092-7d3640b6ddf5
# ╠═71d2cdaa-28fb-11eb-3826-cf35b99cde1f
# ╠═ed43435a-85d2-11eb-0aea-2ffa1acc2cc5
# ╠═c65f5f04-8770-11eb-2407-05dbb86e82b7
# ╠═8b5eb0ee-881a-11eb-022c-adf0d7440e84
# ╠═2b2c55a8-8821-11eb-0176-91dca1254a8b
# ╠═edac53be-881a-11eb-1dbd-a189c57ac336
# ╠═5480e430-876a-11eb-0c63-aff9c5cf17d2
# ╠═ae0c69d8-85d6-11eb-2db9-01a3f56c9a65
# ╟─22d43054-8661-11eb-0dd7-fb846edb81c4


text:
# ---
# title: 445. Add Two Numbers II
# id: problem445
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/add-two-numbers-ii/description/>
# hidden: true
# ---
#
# You are given two **non-empty** linked lists representing two non-negative
# integers. The most significant digit comes first and each of their nodes
# contain a single digit. Add the two numbers and return it as a linked list.
#
# You may assume the two numbers do not contain any leading zero, except the
# number 0 itself.
#
# **Follow up:**
# What if you cannot modify the input lists? In other words, reversing the lists
# is not allowed.
#
# **Example:**
#
#
#
#     Input: (7 -> 2 -> 4 -> 3) + (5 -> 6 -> 4)
#     Output: 7 -> 8 -> 0 -> 7
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "compute_JY.jl" begin

end


text:
function get_name(device::D) where D <: PSY.Device
    return device.name
end

"""
This function returns all generation and storage technologies included in the PSY System.
"""
function get_all_techs(sys::PSY.System)
    sys_gens = PSY.get_components(PSY.Generator, sys)
    sys_storage = PSY.get_components(PSY.Storage, sys)
    sys_techs = union(sys_gens, sys_storage)
    return sys_techs
end

"""
This function does nothing if the PSY System is not defined.
"""
function get_system_services(sys::Nothing)
    services = PSY.Service[]
    return services
end

"""
This function returns all Services modeled in the PSY System.
"""
function get_system_services(sys::PSY.System)
    sys_techs = get_all_techs(sys)

    services = PSY.Service[]
    for tech in sys_techs
        for service in PSY.get_services(tech)
            if !in(service, services)
                push!(services, service)
            end
        end
    end

    return services
end

"""
This function does nothing if Service is not of ReserveUp type.
"""
function add_device_services!(services::Vector{PSY.Service},
                              device::PSY.Device,
                              product::P) where P <: Product
    return
end

"""
This function adds ReserveUp service for PSY Devices.
"""
function add_device_services!(services::Vector{PSY.Service},
                              device::PSY.Device,
                              product::OperatingReserve{ReserveUpEMIS})

    device_zone = PSY.get_name(PSY.get_load_zone(PSY.get_bus(device)))

    for service in services
        if typeof(service) == PSY.VariableReserve{PSY.ReserveUp}
            service_zone = "zone_$(split(PSY.get_name(service), "_")[end])"
            if device_zone == service_zone
                PSY.add_service!(device, service)
            end
        end
    end

    return
end

"""
This function adds ReserveDown service for PSY Devices.
"""
function add_device_services!(services::Vector{PSY.Service},
                              device::PSY.Device,
                              product::OperatingReserve{ReserveDownEMIS})

    device_zone = PSY.get_name(PSY.get_load_zone(PSY.get_bus(device)))

    for service in services
        if typeof(service) == PSY.VariableReserve{PSY.ReserveDown}
            service_zone = "zone_$(split(PSY.get_name(service), "_")[end])"
            if device_zone == service_zone
                PSY.add_service!(device, service)
            end
        end
    end

    return
end

"""
This function does nothing if Device is not of RenewableGen type.
"""
function add_device_forecast!(simulation_dir::String,
                              sys::PSY.System,
                              device::D,
                              availability_df::Vector{Float64},
                              start_year::Int64
                              ) where D <: Union{PSY.ThermalGen, PSY.HydroGen, PSY.Storage}
    return
end

"""
This function adds forecast timeseries if Device is of RenewableGen type.
"""
function add_device_forecast!(simulation_dir::String,
                              sys::PSY.System,
                              device::D,
                              availability_raw::Vector{Float64},
                              start_year::Int64) where D <: PSY.RenewableGen

    YearAhead = collect(
    Dates.DateTime("1/1/$(start_year)  0:00:00", "d/m/y  H:M:S"):Dates.Hour(1):Dates.DateTime(
        "31/12/$(start_year)  23:00:00", "d/m/y  H:M:S",),)

        timeseries = TS.TimeArray(YearAhead, availability_raw)

        PSY.add_forecast!(sys, device, IS.Deterministic("get_rating", timeseries))

        load_n_vg_df =  read_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data", "load_n_vg_data.csv"))
        load_n_vg_df[:, get_name(device)] = availability_raw * get_device_size(device) * PSY.get_base_power(sys)

        write_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data"), "load_n_vg_data.csv", load_n_vg_df)

    return
end

"""
This function does nothing if PSY System is not defined.
"""
function update_PSY_timeseries!(sys_UC::Nothing,
                               load_growth::AxisArrays.AxisArray{Float64, 1},
                               days::Int64)
    return
end

"""
This function updates the PSY load and reserve requirment timeseries each year.
"""
function update_PSY_timeseries!(sys::PSY.System,
                               load_growth::AxisArrays.AxisArray{Float64, 1},
                               num_days::Int64)

    num_hours = 24
    initial_time = PSY.get_forecasts_initial_time(sys)

    # update load forecasts.
    nodal_loads = PSY.get_components(PSY.ElectricLoad, sys)
    for load in nodal_loads
        zone = "zone_$(PSY.get_name(PSY.get_area(PSY.get_bus(load))))"
        ts_data= PSY.get_data(PSY.get_forecast(
            PSY.Deterministic,
            load,
            initial_time,
            "get_max_active_power",
            num_days * num_hours,
        ))
        ts_timestamps = TS.timestamp(ts_data)
        ts_values = TS.values(ts_data)

        PSY.remove_forecast!(PSY.Deterministic,
            sys,
            load,
            initial_time,
            "get_max_active_power",
            )

        new_ts = TS.TimeArray(ts_timestamps, ts_values * (1 + load_growth[zone]))
        PSY.add_forecast!(sys, load, IS.Deterministic("get_max_active_power", new_ts))

    end

    # update service requirement forecasts.
    services = get_system_services(sys)
    for service in services
        zone = "zone_$(last(split(PSY.get_name(service), "_")[end], 1))"

        ts_data= PSY.get_data(PSY.get_forecast(
            PSY.Deterministic,
            service,
            initial_time,
            "get_requirement",
            num_days * num_hours,
        ))
        ts_timestamps = TS.timestamp(ts_data)
        ts_values = TS.values(ts_data)

        PSY.remove_forecast!(PSY.Deterministic,
            sys,
            service,
            initial_time,
            "get_requirement",
            )

        new_ts = TS.TimeArray(ts_timestamps, ts_values * (1 + load_growth[zone]))

        PSY.add_forecast!(sys, service, IS.Deterministic("get_requirement", new_ts))
    end

    return
end

"""
This function finds the buses located in each zone.
"""
function find_zonal_bus(zone::String, sys::PSY.System)
    buses = PSY.get_components(PSY.Bus, sys)

    zonal_bus = nothing

    for b in buses
        if "zone_$(PSY.get_name(PSY.get_area(b)))" == zone
            zonal_bus = b
        end
    end

    return zonal_bus

end


text:
using DashBootstrapComponents, DashCoreComponents

dropdown = dbc_formgroup([
    dbc_label("Dropdown", html_for="dropdown"),
    dcc_dropdown(
        id="dropdown",
        options=[
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
    ),
]);

slider = dbc_formgroup([
    dbc_label("Slider", html_for="slider"),
    dcc_slider(id="slider", min=0, max=10, step=0.5, value=3),
]);

range_slider = dbc_formgroup([
    dbc_label("RangeSlider", html_for="range-slider"),
    dcc_rangeslider(id="range-slider", min=0, max=10, value=[3, 7]),
]);

form = dbc_form([dropdown, slider, range_slider]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""
    Bitmap3D(grid::Grid3d)

Create an empty `Bitmap3D` from an grid.
"""
Bitmap3D(grid::Grid3D) = Bitmap3D(falses(size(grid)...), grid)

Base.copy(M::Bitmap3D) = Bitmap3D(copy(M.data), grid)

Base.size(M::Bitmap3D) = size(M.grid)
Base.size(M::Bitmap3D, i::Int64) = Base.size(M.grid, i)
Base.ndims(M::Bitmap3D) = 3
Base.getindex(M::Bitmap3D, i) = getindex(M.data, i)
Base.getindex(M::Bitmap3D, i, j, k) = getindex(M.data, i, j, k)
Base.setindex!(M::Bitmap3D, x::Bool, i) = setindex!(M.data, x, i)
Base.setindex!(M::Bitmap3D, x::Bool, i, j, k) = setindex!(M.data, x, i, j, k)
Base.eachindex(M::Bitmap3D) = eachindex(M.data)

function safe_getindex(M::Bitmap3D, i, j, k)
    if 1 ≤ i ≤ size(M.data, 1) &&
       1 ≤ j ≤ size(M.data, 2) &&
       1 ≤ k ≤ size(M.data, 3)
        return getindex(M, i, j, k)
    end
    return false
end
function safe_setindex!(M::Bitmap3D, x::Bool, i, j, k)
    if 1 ≤ i ≤ size(M.data, 1) &&
       1 ≤ j ≤ size(M.data, 2) &&
       1 ≤ k ≤ size(M.data, 3)
        setindex!(M, x, i, j, k)
    end
end

@recipe function plot(M::Bitmap3D; boundarycolor=:dodgerblue)
    inner, boundary = _scatterdata(M)

    @series begin
        seriescolor --> :dodgerblue
        legend --> false
        st := :scatter3d
        inner
    end

    @series begin
        seriescolor --> boundarycolor
        legend --> false
        st := :scatter3d
        boundary
    end
end

function _scatterdata(D::Bitmap3D)
    xs = Float64[]
    ys = Float64[]
    zs = Float64[]

    bxs = Float64[]
    bys = Float64[]
    bzs = Float64[]
    for l in eachindex(D.data)
        if D[l]
            if !isboundary_point(D.data, l)
                i, j, k = ind2sub(D.data, l)
                push!(xs, D.grid.xrange[i])
                push!(ys, D.grid.yrange[j])
                push!(zs, D.grid.zrange[k])
            else
                i, j, k = ind2sub(D.data, l)
                push!(bxs, D.grid.xrange[i])
                push!(bys, D.grid.yrange[j])
                push!(bzs, D.grid.zrange[k])
            end
        end
    end
    (xs, ys, zs), (bxs, bys, bzs)
end

function isboundary_point(data, l)
    i, j, k = ind2sub(data, l)
    i0 = max(i - 1, 1)
    i1 = min(i + 1, size(data, 1))

    j0 = max(j - 1, 1)
    j1 = min(j + 1, size(data, 2))

    k0 = max(k - 1, 1)
    k1 = min(k + 1, size(data, 3))
    for ii = i0:i1, jj = j0:j1, kk=k0:k1
        if !data[ii, jj, kk]
            return true
        end
    end
    return false
end


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

t4 = time()
using KernelAbstractions

using Tracker

GRAD = :Tracker
_gradient(x...) = Tracker.gradient(x...)

@testset "KernelAbstractions + gradients" begin
    A = (rand(3,4));
    B = (rand(4,5));
    @tullio C[i,k] := A[i,j] * B[j,k]  threads=false  # verbose=2
    @test C ≈ A * B

    @tullio threads=false # else KernelAbstractions CPU kernels not used
    include("gradients.jl")
    @tullio threads=true

    for sy in Tullio.SYMBOLS
        @test !isdefined(@__MODULE__, sy)
    end
end

using CUDA

if is_buildkite
    # If we are on Buildkite, we should assert that we have a CUDA GPU available
    @test CUDA.has_cuda_gpu()
end

if CUDA.has_cuda_gpu()
    @info "===== found a GPU, starting CUDA tests ====="
    @testset "===== CUDA tests on GPU =====" begin
        include("cuda.jl")
    end
end

@info @sprintf("KernelAbstractions tests took %.1f seconds", time()-t4)

@tullio cuda=false


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for utils

############################################################################################################################

@testset "Display" begin
   g = completegraph(10)
   setlabel!(g, map(string, 1:10))

   # Set vertex properties
   setvprop!(g, :, 1 : 10, :p1)
   setvprop!(g, :, 1 : 10, :p2)
   setvprop!(g, :, 1 : 10, :p3)
   setvprop!(g, :, 1 : 10, :p4)

   # Set edge properties
   seteprop!(g, :, 1 : 90, :p1)
   seteprop!(g, :, 1 : 90, :p2)
   seteprop!(g, :, 1 : 90, :p3)
   seteprop!(g, :, 1 : 90, :p4)

   V, E = g
   @test isa(V, VertexDescriptor)
   @test isa(E, EdgeDescriptor)

   ss = IOBuffer()

   # Show
   println(ss, V)
   println(ss, E)

   # Getindex V
   println(ss, V["1"])
   try println(ss, V["11"]); @test false catch @test true end

   # Getindex E
   println(ss, E["1", "2"])
   try println(ss, V["1", "1"]); @test false catch @test true end
end


text:
#typealias MultiindexSet{T} Array{T}
#typ MultiindexSet{T} <: AbstractArray{T, 2}

function multiindex{T<:Integer}(m::Integer,p::T)
    if m==0
        I = zeros(T,1,0)
        return I
    end
    I = zeros(T,0,m)
    I1 = multiindex(m-1, p)
    ord = reshape(sum(I1,2), size(I1,1))
    for q=0:p
        ind=map(o -> (o<=q), ord)
        n=sum(ind)
        Iq = [I1[ind,:] (q-sum(I1[ind,:],2))]
        I = [I; Iq]
    end
    return I
end

function multiindex_order{T<:Integer}(I::AbstractArray{T})
  return vec(sum(I,2))
end

export
    multiindex,
    multiindex_order


text:
## Permutation Population Generators

#random
function Gen_perm(N::Int)
  function rperm(N::Int)
    return randperm(N)
  end
end

# Kmedian epsilon with K minpts
function Gen_Kmedian{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_Kmedian(A,k)
  return Gen_clust(A,ϵ,k)
end

# Kmean epsilon with K minpts
function Gen_Kmean{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_Kmean(A,k)
  return Gen_clust(A,ϵ,k)
end

# Grid estimate (uniform) epsilon with K minpts
function Gen_Kgrid{T <: Real}(A::Matrix{T},k::Int)
  ϵ = eps_grid(A)
  return Gen_clust(A,ϵ,k)
end

# Sample random from K [median,85%] epsilon with K minpts
function Gen_Km85{T <: Real}(A::Matrix{T},k::Int)
  eps = eps_randm85(A,k)
  return Gen_clust(A,eps,k)
end

# Generic cluster generator function
function Gen_clust{T<:Real}(D::DenseMatrix{T}, ϵ::Real, minpts::Int)
  rslt = dbscan(D,ϵ,minpts)
  fillZero!(rslt.assignments)
  cmap = clustmap(rslt.assignments)
  k = length(keys(cmap))
  function clustPop(N::Int)
    tgt = zeros(Int,N)
    p = 1
    for i in randperm(k)
      l = length(cmap[i])
      tgt[p:p+l-1] = cmap[i]
      p += l
    end
    return tgt
  end
  return clustPop
end

#nearest neighbor generator
function Gen_NN{T<:Real}(D::DenseMatrix{T})
  function nnPop(N::Int)
    nearest_neighbor(D,rand(1:N))
  end
  return nnPop
end


text:
# Quelques méthodes ou macros utilitaires personnelles
# (indépendantes de tout projet julia)
# Modif le 21/11/2019

"""
Démarre le chrono si nécessaire (la première fois) et retourne la durée
écoulée depuis ce démarrage (à la milli-seconde prêt).
"""
function ms()
    global MS_START
    if ! @isdefined MS_START
        MS_START = time()
    end
    ms = round(Int, 1000*(time() - MS_START))
    return (ms/1000)
end

function ms_reset()
    global MS_START = time()
end

"""
Chronomètre l'exécution d'une commande :

- affiche la durée depuis le lancement du programme (peut-être long en
  mode interactife !)
- affiche le fichier depuis lequel elcette macro est appelée
- affiche la commande à exécuter
- exécute cette commande
- affiche la durée d'exécution sur la ligne suivante

Exemple d'affichage en interactif :
```
@ms sleep(1)
=>
@ms TODO 9349.527s (mode repl):sleep(1) ...
@ms done 9350.552s (mode repl):sleep(1) en 1.002s
```

Exemples depuis en script :
```
@ms using JuMP
=>
@ms TODO 0.871s (mode repl):using JuMP ...
@ms done 2.198s (mode repl):using JuMP en 1.327s
```
"""
macro ms(cmd)
    cmdstr = string(cmd)
    # info sur le fichier qui a appelé cette macro
    if isinteractive()
        # mode interactif => pas de fichier d'appel
        fileinfo = "(mode repl)"
    else
        # On récupère le nom relatif du fichier qui a appellé cette macro
        fname = string(__source__.file)
        fname = fname[length("$APPDIR")+2 : end]
        # On récupère le numéro de ligne de l'appel
        fline = string(__source__.line)
        fileinfo = "$(fname):$fline"
    end

    quote
        # la notation $(cmd) est réservée dans un environnement quote
        local t0 = time()
        # println("@ms ", round(ms(), digits=3), "s ", $fileinfo, ":", $(cmdstr), " ... ") # ok 21/11/2019
        print("@ms TODO ", round(ms(), digits=3), "s ", $fileinfo, ":", $(cmdstr))
        println(" ... ")
        local val = $(esc(cmd)) # <= EXECUTION DU CALCUL Á CHRONOMÉTRER
        local t1 = time()
        # println("    => fait en ", round(t1-t0, digits=3), "s", " total=", ms(), "s")
        print("@ms done ", round(ms(), digits=3), "s ", $fileinfo, ":", $(cmdstr))
        println(" en ", round(t1-t0, digits=3), "s")
        val
    end
end
# On démarre le chronomètre
ms()


text:
using Primes

@testset "Chineese remainder theorem" begin

    modular_images(a, ms) = map(m -> mod(a, m), ms)

    Ms = [
        [2, 3],
        [3, 7],
        [2^31 - 1, Primes.nextprime(2^31)],
        [Primes.nextprime(BigInt(2)^100), Primes.nextprime(BigInt(2)^31)],
        [Primes.nextprime(BigInt(2)^10000), Primes.nextprime(BigInt(2)^31)]
    ]
    buf, buf1, buf2, buf3, n1, n2 = [BigInt(0) for _ in 1:6]

    for ms in Ms
        M = prod(ms)
        nums = [rand(0:M-1) for _ in 1:1000]

        msinv = [invmod(ms[1], ms[2]), invmod(ms[2], ms[1])]
        for a in nums
            rs = modular_images(a, ms)

            r1, r2 = BigInt.(rs)
            m1, m2 = BigInt.(ms)
            minv1, minv2 = BigInt.(msinv)
            Groebner.CRT!(M, buf, n1, n2, r1, minv1, UInt64(r2), minv2, m1, m2)
            @test buf == a
        end

    end

end


text:
# External packages
using Revise, DelimitedFiles, SymPy

version = "Version0.1"
cdir = pwd();
push!(LOAD_PATH, cdir);

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Julia", cdir);
sourcepath = cdir[1:(ii[1]-1)];
include(sourcepath * "/Installation/setpath.jl");

# MDP packages
using Preprocessing, Gencode, Postprocessing

app = Preprocessing.initializeapp(sourcepath,version);

#using Gencode
# create pde structure and mesh structure
#app = Gencode.initializeapp(version);

app.dim = 3;
app.ncq = 1;
app.nceta = 10;
app.nckappa = 2;
app.ncmu1a = 1;
app.ncmu1b = 1;
app.ncmu2a = 10;
app.ncmu2b = 10;
app.ncmu2c = 10;
app.ncmu3a = 10;
app.ncmu3b = 10;
app.ncmu3c = 10;
app.ncmu4a = 4;
app.ncmu4b = 4;
app.natomtype = 2;
app.potentialfile = "potential.jl";

include(app.potentialfile);  # include the potential file
Gencode.gencode(app); # Generate C++ code

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP