text:
############################################################
# NOTE: dropout layer only takes one input blob. We
# intentionally make this restriction because as an
# in-place layer, if has to be put directly on top
# of the layer that is producing the input blob for
# drop-out layer. If drop-out layer depends on multiple
# blobs coming from different layers, it might be
# difficult or even impossible to decide a correct
# layer layout.
############################################################
@defstruct DropoutLayer Layer (
  name :: AbstractString = "dropout",
  auto_scale :: Bool = true,
  (ratio :: AbstractFloat = 0.5, 0 < ratio < 1),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 1),
)
@characterize_layer(DropoutLayer,
  can_do_bp  => true,
  is_inplace => true
)

type DropoutLayerState{T} <: LayerState
  layer      :: DropoutLayer
  rand_vals  :: Blob

  # for convenience
  ratio      :: T   # layer.ratio
  scale      :: T   # 1 / (1 - layer.ratio)

  etc        :: Any
end

function setup_etc(backend::CPUBackend, layer::DropoutLayer, inputs::Vector{Blob})
  return nothing
end
function setup(backend::Backend, layer::DropoutLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type = eltype(inputs[1])
  rand_vals = make_blob(backend, data_type, size(inputs[1]))

  scale = layer.auto_scale ? 1.0/(1-layer.ratio) : 1

  etc = setup_etc(backend, layer, inputs)
  return DropoutLayerState(layer, rand_vals,
      convert(data_type, layer.ratio), convert(data_type, scale), etc)
end
function destroy_etc(backend::CPUBackend, state::DropoutLayerState)
  # do nothing
end
function shutdown(backend::Backend, state::DropoutLayerState)
  destroy(state.rand_vals)
  destroy_etc(backend, state)
end

function dropout_forward{T}(input::Array{T}, rand_vals::Array{T}, ratio::T, scale::T)
  len = length(input)
  @simd for i = 1:len
    @inbounds input[i] = input[i] * (rand_vals[i] > ratio) * scale
  end
end
function forward(backend::CPUBackend, state::DropoutLayerState, inputs::Vector{Blob})
  rand!(state.rand_vals.data)
  dropout_forward(inputs[1].data, state.rand_vals.data, state.ratio, state.scale)
end

function dropout_backward{T}(grad::Array{T}, rand_vals::Array{T}, ratio::T, scale::T)
  len = length(grad)
  @simd for i = 1:len
    @inbounds grad[i] = grad[i] * (rand_vals[i] > ratio) * scale
  end
end
function backward(backend::CPUBackend, state::DropoutLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  if !isa(diffs[1], NullBlob)
    dropout_backward(diffs[1].data, state.rand_vals.data, state.ratio, state.scale)
  end
end


text:
# ---
# title: 391. Perfect Rectangle
# id: problem391
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Line Sweep
# link: <https://leetcode.com/problems/perfect-rectangle/description/>
# hidden: true
# ---
#
# Given N axis-aligned rectangles where N > 0, determine if they all together
# form an exact cover of a rectangular region.
#
# Each rectangle is represented as a bottom-left point and a top-right point.
# For example, a unit square is represented as [1,1,2,2]. (coordinate of bottom-
# left point is (1, 1) and top-right point is (2, 2)).
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/rectangle_perfect.gif)
#
# **Example 1:**
#
#
#
#     rectangles = [
#       [1,1,3,3],
#       [3,1,4,2],
#       [3,2,4,4],
#       [1,3,2,4],
#       [2,3,3,4]
#     ]
#
#     Return true. All 5 rectangles together form an exact cover of a rectangular region.
#
#
#
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/rectangle_separated.gif)
#
# **Example 2:**
#
#
#
#     rectangles = [
#       [1,1,2,3],
#       [1,3,2,4],
#       [3,1,4,2],
#       [3,2,4,4]
#     ]
#
#     Return false. Because there is a gap between the two rectangular regions.
#
#
#
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/rectangle_hole.gif)
#
# **Example 3:**
#
#
#
#     rectangles = [
#       [1,1,3,3],
#       [3,1,4,2],
#       [1,3,2,4],
#       [3,2,4,4]
#     ]
#
#     Return false. Because there is a gap in the top center.
#
#
#
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/rectangle_intersect.gif)
#
# **Example 4:**
#
#
#
#     rectangles = [
#       [1,1,3,3],
#       [3,1,4,2],
#       [1,3,2,4],
#       [2,2,4,4]
#     ]
#
#     Return false. Because two of the rectangles overlap with each other.
#
#
#
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
module TextWrapTest

using TextWrap
using Compat

if VERSION < v"0.7.0-DEV.1995"
    using Base.Test
else
    using Test
end

text = """
    Julia is a high-level, high-performance dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a sophisticated
    compiler, distributed parallel execution, numerical accuracy, and an
    extensive mathematical function library."""

@test wrap(text) == """
    Julia is a high-level, high-performance dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a sophisticated
    compiler, distributed parallel execution, numerical accuracy, and an
    extensive mathematical function library."""

@test wrap(text, width=30) == """
    Julia is a high-level, high-
    performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=30, fix_sentence_endings=true, break_on_hyphens=false) == """
    Julia is a high-level,
    high-performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments.  It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=30, initial_indent=2, subsequent_indent=0) == """
      Julia is a high-level, high-
    performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=32, initial_indent=">   ", subsequent_indent="> ") == """
    >   Julia is a high-level, high-
    > performance dynamic
    > programming language for
    > technical computing, with
    > syntax that is familiar to
    > users of other technical
    > computing environments. It
    > provides a sophisticated
    > compiler, distributed parallel
    > execution, numerical accuracy,
    > and an extensive mathematical
    > function library."""

tabtext = "aaaaaaa\tbbbbbbb\t\ncccccc\tddddd\teeee  \tfff\tgg\th\n"

@test wrap(tabtext, width=20, replace_whitespace=true,  expand_tabs=true) ==
    "aaaaaaa bbbbbbb\ncccccc  ddddd   eeee\nfff     gg      h"
@test wrap(tabtext, width=20, replace_whitespace=false, expand_tabs=true) ==
    "aaaaaaa bbbbbbb\ncccccc  ddddd   eeee\nfff     gg      h"
@test wrap(tabtext, width=20, replace_whitespace=true,  expand_tabs=false) ==
    "aaaaaaa bbbbbbb\ncccccc ddddd eeee\nfff gg h"
@test wrap(tabtext, width=20, replace_whitespace=false, expand_tabs=false) ==
    "aaaaaaa\tbbbbbbb\ncccccc\tddddd\teeee\nfff\tgg\th"

longwordstext = """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=172) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=171) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτρ
    ιμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων. The longest Hebrew word is the 25-letter-
    long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=80) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest
    English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft,
    was named the longest published word in the German language by the 1996 Guinness
    Book of World Records, but longer words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσε
    λαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφ
    οφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels)
    וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our
    encyclopedias will have....'"""

@test wrap(longwordstext, width=79) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest
    English word that appears in a major dictionary. A 79 letter word, Donaudampfsc
    hiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named
    the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c.
    392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρανιολε
    ιψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλ
    εκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων. The longest
    Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will
    have....'"""

@test wrap(longwordstext, width=45) == """
    The 45-letter word
    pneumonoultramicroscopicsilicovolcanoconiosis
    is the longest English word that appears in a
    major dictionary. A 79 letter word, Donaudamp
    fschiffahrtselektrizitätenhauptbetriebswerkba
    uunterbeamtengesellschaft, was named the
    longest published word in the German language
    by the 1996 Guinness Book of World Records,
    but longer words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes coined
    the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρ
    ανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκ
    ατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυ
    ονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγ
    ανοπτερύγων. The longest Hebrew word is the
    25-letter-long (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means 'and
    when our encyclopedias will have....'"""

@test wrap(longwordstext, width=44) == """
    The 45-letter word pneumonoultramicroscopics
    ilicovolcanoconiosis is the longest English
    word that appears in a major dictionary. A
    79 letter word, Donaudampfschiffahrtselektri
    zitätenhauptbetriebswerkbauunterbeamtengesel
    lschaft, was named the longest published
    word in the German language by the 1996
    Guinness Book of World Records, but longer
    words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes
    coined the 171-letter word λοπαδοτεμαχοσελαχ
    ογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαρα
    βομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπερισ
    τεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιρ
    αιοβαφητραγανοπτερύγων. The longest Hebrew
    word is the 25-letter-long (including
    vowels) וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means
    'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=19) == """
    The 45-letter word
    pneumonoultramicros
    copicsilicovolcanoc
    oniosis is the
    longest English
    word that appears
    in a major
    dictionary. A 79
    letter word, Donaud
    ampfschiffahrtselek
    trizitätenhauptbetr
    iebswerkbauunterbea
    mtengesellschaft,
    was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes coined
    the 171-letter word
    λοπαδοτεμαχοσελαχογ
    αλεοκρανιολειψανοδρ
    ιμυποτριμματοσιλφιο
    καραβομελιτοκατακεχ
    υμενοκιχλεπικοσσυφο
    φαττοπεριστεραλεκτρ
    υονοπτοκεφαλλιοκιγκ
    λοπελειολαγῳοσιραιο
    βαφητραγανοπτερύγων
    . The longest
    Hebrew word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklope
    dioténu), which
    means 'and when our
    encyclopedias will
    have....'"""

@test wrap(longwordstext, width=18) == """
    The 45-letter word
    pneumonoultramicro
    scopicsilicovolcan
    oconiosis is the
    longest English
    word that appears
    in a major
    dictionary. A 79
    letter word, Donau
    dampfschiffahrtsel
    ektrizitätenhauptb
    etriebswerkbauunte
    rbeamtengesellscha
    ft, was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes
    coined the
    171-letter word λο
    παδοτεμαχοσελαχογα
    λεοκρανιολειψανοδρ
    ιμυποτριμματοσιλφι
    οκαραβομελιτοκατακ
    εχυμενοκιχλεπικοσσ
    υφοφαττοπεριστεραλ
    εκτρυονοπτοκεφαλλι
    οκιγκλοπελειολαγῳο
    σιραιοβαφητραγανοπ
    τερύγων. The
    longest Hebrew
    word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינ
    ו (u'chshelentsikl
    opedioténu), which
    means 'and when
    our encyclopedias
    will have....'"""

@test wrap(longwordstext, width=18, break_long_words=false) == """
    The 45-letter word
    pneumonoultramicroscopicsilicovolcanoconiosis
    is the longest
    English word that
    appears in a major
    dictionary. A 79
    letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft,
    was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes
    coined the
    171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew
    word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu),
    which means 'and
    when our
    encyclopedias will
    have....'"""

tmpf = tempname()
try
    open(tmpf, "w") do f
        print_wrapped(f, text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library."""

    open(tmpf, "w") do f
        println_wrapped(f, text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library.
        """

    open(tmpf, "w") do f
        println_wrapped(f, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == "\n"

    open(tmpf, "w") do f
        bk_stdout = stdout
        try
            redirect_stdout(f)
            print_wrapped(text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
        finally
            redirect_stdout(bk_stdout)
        end
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library."""
finally
    isfile(tmpf) && rm(tmpf)
end

@test_throws ErrorException wrap("", initial_indent=10, width=10)
@test_throws ErrorException wrap("", subsequent_indent=10, width=10)
@test_throws ErrorException wrap("", initial_indent="~~~~~~~~~~", width=10)
@test_throws ErrorException wrap("", subsequent_indent="~~~~~~~~~~", width=10)

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 43c7192c-aa56-4d45-9f47-40964f0e14b6
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ c9e98b87-7187-4a38-acab-1a93519160c2
md"""
# Collections, broadcasting, and Methods
Using built-in collections is the simplest way to group and organize data.

The values in a `immutable` collection cannot be updated after its creation, while in a `mutable` collection can.

The elements in `sequential` collections are accessed by integer indices, while those in `associative` collection are accessed by keys.
"""

# ╔═╡ e530e669-35e1-4314-a080-71c5faf8f348
md"""
# Sequential collections

General rules for sequential collections:

- 1-based indexing, as in R, MATLAB, and Fortran.
- Elements are accessed by an integer index `seq[i]` or an integer range `seq[1:2:end-1]`.
- `length(seq)` returns the total size
- Splat operator `...` passes the inner contents in the collection as positional function arguments.
- Dot syntax (e.g. a .+ b) performs element-wise  / broadcasting operations.
"""

# ╔═╡ 2be423ec-49c5-4e72-98df-494cd8270501
md"""
## Ranges

`start[:step]:end`

- Immuatable
- Sequencial
- Eheap
- Evenly-spaced numerical sequences
"""

# ╔═╡ 7197fe35-797d-401d-b05c-a8eefa3f9187
1:2:10

# ╔═╡ 09a51e04-df11-4158-819e-d5f908104ebb
length(1:2:10)

# ╔═╡ 273756ef-4bab-4963-b19f-d8077a6a5a6b
with_terminal() do
	dump(1:10)
end

# ╔═╡ 014958a2-b5a1-4d54-a6b9-3641d1e19f06
# linspace() equivalent
range(1, 10, length=10)

# ╔═╡ 18147773-faf2-409a-90a1-610c15b5f23a
# Another linspace() equivalent
LinRange(1.5, 5.5, 9)

# ╔═╡ 9959d91f-7edc-45ed-aaff-9c2492b8a4a3
(1:10)[3:end]

# ╔═╡ 0b76dbb0-b0c7-4556-bf61-1f9bc5839877
# Oftedn used with for loops
with_terminal() do
	for i in 1:10
		println(i)
	end
end

# ╔═╡ 1ed78b50-4f52-49d5-a936-5bc67ad63cf9
md"""
## Tuples

- immutable
- sequential collections.
- efficient for heterogenous data
- stack-allocated
"""

# ╔═╡ 32dc0e41-9819-46cd-ab7b-00b10c7c0c5c
tuple(1, 'a', 3.14)

# ╔═╡ 2669a9f0-edb3-47bd-94e6-0689d24bf5fc
(1, 'a', 3.14)

# ╔═╡ 706ad47a-a1ce-4bfd-901d-408518f69920
t1 = (1, 2, 3)

# ╔═╡ 0b8aa474-052b-4b96-8a31-ce5148b9758d
t1[1]

# ╔═╡ d069122d-9a09-43db-a3f0-adb61c5426bc
t2 = (1, 'a', 3.14)

# ╔═╡ 56f6008a-d4ad-483a-9e63-76b8689634c6
with_terminal() do
	dump(t2)
end

# ╔═╡ 9e7f20b4-dadf-4708-a521-e9c4e28424fa
tuple(t1..., t2...)   # Merging multiple tuples using ... operator

# ╔═╡ f049f6bb-d8a7-4fcd-b352-d063fb6bf838
# Tuples are often used to swap elements
let x = 1, y = 2, z = 3
	x, y, z = z, x, y
	with_terminal() do
		@show x, y, z
	end
end

# ╔═╡ 31595906-9558-49cc-9911-087c33e8359b
# Tuples are used for returning multiple values in a function
extrema([1, 5, 6, 7, -1, -3, 0])

# ╔═╡ 25be6939-069f-41fc-a50e-b6b9f15f5dcb
sincospi(1//2)

# ╔═╡ ddf4322e-0741-40c9-ba08-785c15683e6d
neighbors(x) = x+1, x-1

# ╔═╡ 7c8e8c0e-6400-4a36-8bba-9e377d6063e4
neighbors(0)

# ╔═╡ fd0da8c5-7941-4a0d-b056-290c6c3bc2fc
md"""
## Arrays

[Arrays | Julia docs](https://docs.julialang.org/en/v1/manual/arrays/)

`[seq...]` / `collect(seq)`

Arryas are the bread and butter for scientific computing, similar to numpy's ndarrays.

- Column-major (Fortran style) rather than row-major (C and numpy style)
- Assignments and updating may cuase unwanted editing due to memory sharing.

Some useful functions for arrays:

- `length(A)` the number of elements in A
- `ndims(A)` the number of dimensions of A
- `size(A)` a tuple containing the dimensions of A
- `size(A,n)` the size of A along dimension n
- `eachindex(A)` an efficient iterator for visiting each position in A
"""

# ╔═╡ b8ac6869-52f8-4407-8b31-61a2e67ce140
# 1D array (column vector)
x = [5, 6, 7]

# ╔═╡ c2091825-da38-44ee-bbc1-1fbffed7b058
# np.arange() equivalent
collect(1:10)

# ╔═╡ 0a884ee7-ea0b-483e-a830-c39de038ba6d
zeros(2, 5, 2)

# ╔═╡ 13ca7c56-7665-484e-8feb-566cddaaaec0
ones(2, 5)

# ╔═╡ 5bcd3655-7ea4-4541-8962-52d0defa0abb
# Uninitialized array with the same data type and dims as x
similar(x)

# ╔═╡ dd3eb3bc-ba9b-48e1-b118-2420088ea1ac
# np.zeros_like()
zero(x)

# ╔═╡ 69f9d623-6f2f-4d6e-8083-1dd53d1ce60c
rand(1:6, 10)

# ╔═╡ 23a392e1-6a69-4677-9550-9c7f227b8d02
reshape(1:12, 3, 4)

# ╔═╡ 2d45c850-f4db-47a8-a5e6-24182f06c5d8
# repeat 3x2 times
repeat([1 2; 3 4], 3, 2)

# ╔═╡ e7c74fe1-0cf5-48db-9833-7ec0d7dfabfc
# comprehension
[i^2 for i in 1:10 if i >= 5]

# ╔═╡ e8427218-1e30-4bbf-a0a0-516b35800a82
# a comprehension generating 2x3 array
[x * y for x in 1:2, y in 1:3]

# ╔═╡ 3ab98217-b14d-4292-8540-42c98bfbcd3c
# casting comprehension result element type to Float64
Float64[x^2 for x in 1:4]

# ╔═╡ 0e7e8893-bf68-4cf6-800f-d02820d38d08
# a 1-element tuple containing a vector
tuple([1,2,3])

# ╔═╡ 6ce997e4-796f-4a8b-bae3-e9cd8995ae10
# convert ector to tuple
Tuple([1,2,3])

# ╔═╡ 1f695f9a-0cb2-4c81-9684-3c005413d36f
# 2D array (matrix)
# space is a shorthand for hcat()
# semicolone is a shorthand for vcat()
A = [1 2 3;
     4 5 6]

# ╔═╡ d73fc237-5cb6-4f80-8dd8-6c5dec1c40b2
A[1, 2]

# ╔═╡ 5d62c23a-22d8-4c54-818a-46b900daebcd
length(A)

# ╔═╡ d7d75522-a599-4473-b13b-5f490e10bb04
axes(A)

# ╔═╡ c77bf597-35d7-4045-aac0-4e4f36471f8c
size(A)

# ╔═╡ ea81441b-a7a6-44ec-8740-712046f99698
ndims(A)

# ╔═╡ 0ef49715-6933-4ff7-b01f-b546cbf6138f
transpose(A)

# ╔═╡ ad593149-5905-41d6-b01b-34bdbf73e6a5
A'

# ╔═╡ c6f01603-6444-4421-bcae-a7ba7c612f0e
# Matrix multiplication
with_terminal() do
	@show A x A*x
end

# ╔═╡ 47f22551-defe-4f5a-be80-e0cca177a493
b = A * x

# ╔═╡ 682ffe9f-4419-4bc8-869b-0d9ed78bd19e
# Find x for Ax = b, using left division operator `/`
A\b ≈ x

# ╔═╡ fafda1ec-d538-4c12-ba56-bbd9817da0a2
# Convert A to an (1D) vector, reuses memory
vec(A)

# ╔═╡ b5b39dfc-2944-45aa-9612-c76caafaa8da
# Arrays are mutable (i.e. you can update the contents)
begin
	c = copy(x)  # Make a copy to keep the original
	c[1] = 0
	c
end

# ╔═╡ 2db6bf29-2c81-46ca-bae9-41e8f21c5be2
with_terminal() do
	@show x[1] c[1]
end

# ╔═╡ 10d89d76-b93a-4647-bc83-9be329d83761
md"""
# Associative collections

- `d[key]` accesses values by keys
- `d[key] = value` sets a key-value pair for a mutable dictionary.
- `delete!(d, key)` deletes the kay (and its partner) from a mutable dictionary.
- `keys(d)` returns a series of keys
- `values(d)` returns a series of values
- `pairs(d)` returns a series of (key => value) pairs
- `merge(d1, d2, ...)` return combinations of several dicts. `merge!(d1, d2, ...)` combine several dicts and update the first one.
- `get(d, key, default)` returns the value stored for the given key, or the given default value if no mapping for the key is present.
"""

# ╔═╡ ac0c1c30-2aab-4e50-a622-451b0d63847c
md"""
## Named tuples

`Namedtuples` are tuples with key-value pairs.
"""

# ╔═╡ 37134b88-13d2-4d50-8873-ee96b159c98c
nt = (a=1, b=2, c=4)

# ╔═╡ 94218c2f-1fee-445a-b29d-5f6447d4df82
nt.a == nt[:a]

# ╔═╡ b361fcf7-2656-4dea-b2fd-8d51210da2c2
nt.a == nt[1]

# ╔═╡ 0d8c3b58-37f9-49fa-986f-57efcf62a196
nt[1]

# ╔═╡ 8b00db79-1127-4df4-af43-afabaf0e2ba6
md"""
## Dictionaries

[Dictionaries](https://docs.julialang.org/en/v1/base/collections/#Dictionaries) are mutable mappings of key => value.
"""

# ╔═╡ c81aa4bf-ec76-4c1b-80b8-dd80c12b0582
eng2sp = Dict("one" => "uno", "two" => "dos", "three" => "tres")

# ╔═╡ 0e0c298b-6d1a-4ba5-a4c9-ee9f84e4f160
with_terminal() do
	@show eng2sp["two"] keys(eng2sp) values(eng2sp) get(eng2sp, "one", "N/A") get(eng2sp, "four", "N/A") haskey(eng2sp, "one")

	println()

	for (k ,v) in eng2sp
    	println(k, " => ", v)
	end
end

# ╔═╡ f5f32721-5198-4166-bbe3-0b7f5b12c759
# Creating a dict from an arrya of tuples
Dict([('a', 1), ('c', 3), ('b', 2)])

# ╔═╡ 3f94e6b7-8106-432d-840c-b5398fa0e031
# Creating a Dict via a generator (similar to comprehensions)
Dict(i => i^2 for i = 1:10)

# ╔═╡ 77c3a713-f279-4fbd-91ad-47f85a24830f
Dict(zip("abc", 1:3))

# ╔═╡ 0adcc811-09c4-4566-a5e4-640e00a5ded7
md"""
# Broadcasting (Dot) syntax

[Broadcasting](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting) turns scalar operations into vector ones.
"""

# ╔═╡ f330a6a1-5b29-456b-b936-421a4487ccea
[1, 2, 3] .+ [4, 5, 6]

# ╔═╡ 83e2b4c6-58d0-447f-b2a8-8484c04118c0
sinpi.([0.5, 1.0, 1.5, 2.0])

# ╔═╡ 16003328-ef05-4276-9fff-672789913bf0
let
	ks = (:a, :b, :c)
	vs = (1, 2, 3)
	Dict(ks .=> vs)
end

# ╔═╡ e49bd8bf-ea41-4e7a-a448-dbaa79a51888
# logspace in julia
exp10.(LinRange(-3.0, 3.0, 50))

# ╔═╡ 9f5ffe9b-0687-4925-8fcb-bf4873d8de87
# Make a 9*9 multiplication table
collect(1:9) .* transpose(collect(1:9))

# ╔═╡ df975e94-3dd2-447c-8e02-8bb25926d49c
md"""
# Custom structures
"""

# ╔═╡ fa79434d-8e3d-4acf-a1a7-2eb78d1bde35
begin
	struct Point  # struct or mutable struct
		x
		y
	end

	Point() = Point(0.0, 0.0)
end

# ╔═╡ 9d47fdb6-8fbf-4165-9bda-c65e65872774
p1 = Point(1.0, 2.0)

# ╔═╡ e0f3583a-4b3f-4010-bafa-964fa4fa8fe6
p2 = Point(-3.0, 2.0)

# ╔═╡ 7339f1ac-1ae6-4b96-9cce-b2750d7e82ca
add(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y)

# ╔═╡ eaeddf29-da61-4651-8bf6-eee575a709ae
add(p1, p2)

# ╔═╡ adb15ef5-1a8e-43fc-bffc-4cb73481ab11
md"""
# Methods

[Julia docs | Methods](https://docs.julialang.org/en/v1/manual/methods/#Methods)

You can overload the same function with different argument types / numbers.
"""

# ╔═╡ 623117c3-3c61-496e-9a91-cd497c2d4d01
begin
	func(a::Int) = a + 2
	func(a::AbstractFloat) = a/2
	func(a::Rational) = a/11
	func(a::Complex) = sqrt(a)
	func(a, b::String) = "$a, $b"
end

# ╔═╡ ad345b30-2ed7-4d92-91f8-576784d9fe00
func(1)

# ╔═╡ c0fab2cb-1477-421c-a7a9-353f886ac102
func(3.0)

# ╔═╡ 125a939f-e7a5-4db7-b94d-fdd143b7369b
func(33//4)

# ╔═╡ 6cd6ff60-2b42-4309-8832-35f9b8098827
func(-2 + 0im)

# ╔═╡ b06a9b9c-a1d0-414a-9ccd-320020f1edee
func(true, "it just works and compiles down to optimal code")

# ╔═╡ ae91d145-fc94-4935-970a-c887df151736
md"""
# Package management and dependencies

[Julia docs | Pkg](https://docs.julialang.org/en/v1/stdlib/Pkg/)

In Pluto notebooks it is recommended to use a standalone cell for dependency management. See **Appendix** for details.

In regular environments `Project.toml` and `Manifest.toml` describe the dependencies.

## Install packages

In the Julia script

```julia
using Pkg

# Function form
Pkg.add("Plots")

# Or use Pkg's special strings
pkg"add Plots"
```

In the Julia REPL:

```julia-repl
] add Plots
```

## See what are installed

```julia
using Pkg
Pkg.status()

# Or
pkg"st" # pkg"status"

```

In the Julia REPL:

```julia-repl
] st
```

## Remove packages

```julia
using Pkg
Pkg.remove("Plots")

# Or
pkg"rm Plots"
```

```julia-repl
] rm Plots
```

## Update installed packages

```julia
using Pkg
Pkg.update()

# Or
pkg"up" # pkg"update"
```

## Create / Use a environment

```julia

using Pkg

# Activate environment in the foldername directory
Pkg.activate("foldername")

# Or
pkg"activate foldername"

# Activate the current working directory
Pkg.activate(".")

# Or
pkg"activate ."

```
"""

# ╔═╡ 1cf811e0-cc58-11eb-2c48-33822012ebf2
md"""
# See also

- [Fall in love with Julia](https://nbviewer.jupyter.org/github/schlichtanders/fall-in-love-with-julia/blob/master/01%20introduction%20-%2001%20getting%20started.ipynb)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

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

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

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

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

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

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

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

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

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

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

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

# ╔═╡ Cell order:
# ╠═c9e98b87-7187-4a38-acab-1a93519160c2
# ╠═e530e669-35e1-4314-a080-71c5faf8f348
# ╠═2be423ec-49c5-4e72-98df-494cd8270501
# ╠═7197fe35-797d-401d-b05c-a8eefa3f9187
# ╠═09a51e04-df11-4158-819e-d5f908104ebb
# ╠═273756ef-4bab-4963-b19f-d8077a6a5a6b
# ╠═014958a2-b5a1-4d54-a6b9-3641d1e19f06
# ╠═18147773-faf2-409a-90a1-610c15b5f23a
# ╠═9959d91f-7edc-45ed-aaff-9c2492b8a4a3
# ╠═0b76dbb0-b0c7-4556-bf61-1f9bc5839877
# ╠═1ed78b50-4f52-49d5-a936-5bc67ad63cf9
# ╠═32dc0e41-9819-46cd-ab7b-00b10c7c0c5c
# ╠═2669a9f0-edb3-47bd-94e6-0689d24bf5fc
# ╠═706ad47a-a1ce-4bfd-901d-408518f69920
# ╠═0b8aa474-052b-4b96-8a31-ce5148b9758d
# ╠═d069122d-9a09-43db-a3f0-adb61c5426bc
# ╠═56f6008a-d4ad-483a-9e63-76b8689634c6
# ╠═9e7f20b4-dadf-4708-a521-e9c4e28424fa
# ╠═f049f6bb-d8a7-4fcd-b352-d063fb6bf838
# ╠═31595906-9558-49cc-9911-087c33e8359b
# ╠═25be6939-069f-41fc-a50e-b6b9f15f5dcb
# ╠═ddf4322e-0741-40c9-ba08-785c15683e6d
# ╠═7c8e8c0e-6400-4a36-8bba-9e377d6063e4
# ╠═fd0da8c5-7941-4a0d-b056-290c6c3bc2fc
# ╠═b8ac6869-52f8-4407-8b31-61a2e67ce140
# ╠═c2091825-da38-44ee-bbc1-1fbffed7b058
# ╠═0a884ee7-ea0b-483e-a830-c39de038ba6d
# ╠═13ca7c56-7665-484e-8feb-566cddaaaec0
# ╠═5bcd3655-7ea4-4541-8962-52d0defa0abb
# ╠═dd3eb3bc-ba9b-48e1-b118-2420088ea1ac
# ╠═69f9d623-6f2f-4d6e-8083-1dd53d1ce60c
# ╠═23a392e1-6a69-4677-9550-9c7f227b8d02
# ╠═2d45c850-f4db-47a8-a5e6-24182f06c5d8
# ╠═e7c74fe1-0cf5-48db-9833-7ec0d7dfabfc
# ╠═e8427218-1e30-4bbf-a0a0-516b35800a82
# ╠═3ab98217-b14d-4292-8540-42c98bfbcd3c
# ╠═0e7e8893-bf68-4cf6-800f-d02820d38d08
# ╠═6ce997e4-796f-4a8b-bae3-e9cd8995ae10
# ╠═1f695f9a-0cb2-4c81-9684-3c005413d36f
# ╠═d73fc237-5cb6-4f80-8dd8-6c5dec1c40b2
# ╠═5d62c23a-22d8-4c54-818a-46b900daebcd
# ╠═d7d75522-a599-4473-b13b-5f490e10bb04
# ╠═c77bf597-35d7-4045-aac0-4e4f36471f8c
# ╠═ea81441b-a7a6-44ec-8740-712046f99698
# ╠═0ef49715-6933-4ff7-b01f-b546cbf6138f
# ╠═ad593149-5905-41d6-b01b-34bdbf73e6a5
# ╠═c6f01603-6444-4421-bcae-a7ba7c612f0e
# ╠═47f22551-defe-4f5a-be80-e0cca177a493
# ╠═682ffe9f-4419-4bc8-869b-0d9ed78bd19e
# ╠═fafda1ec-d538-4c12-ba56-bbd9817da0a2
# ╠═b5b39dfc-2944-45aa-9612-c76caafaa8da
# ╠═2db6bf29-2c81-46ca-bae9-41e8f21c5be2
# ╠═10d89d76-b93a-4647-bc83-9be329d83761
# ╠═ac0c1c30-2aab-4e50-a622-451b0d63847c
# ╠═37134b88-13d2-4d50-8873-ee96b159c98c
# ╠═94218c2f-1fee-445a-b29d-5f6447d4df82
# ╠═b361fcf7-2656-4dea-b2fd-8d51210da2c2
# ╠═0d8c3b58-37f9-49fa-986f-57efcf62a196
# ╠═8b00db79-1127-4df4-af43-afabaf0e2ba6
# ╠═c81aa4bf-ec76-4c1b-80b8-dd80c12b0582
# ╠═0e0c298b-6d1a-4ba5-a4c9-ee9f84e4f160
# ╠═f5f32721-5198-4166-bbe3-0b7f5b12c759
# ╠═3f94e6b7-8106-432d-840c-b5398fa0e031
# ╠═77c3a713-f279-4fbd-91ad-47f85a24830f
# ╠═0adcc811-09c4-4566-a5e4-640e00a5ded7
# ╠═f330a6a1-5b29-456b-b936-421a4487ccea
# ╠═83e2b4c6-58d0-447f-b2a8-8484c04118c0
# ╠═16003328-ef05-4276-9fff-672789913bf0
# ╠═e49bd8bf-ea41-4e7a-a448-dbaa79a51888
# ╠═9f5ffe9b-0687-4925-8fcb-bf4873d8de87
# ╠═df975e94-3dd2-447c-8e02-8bb25926d49c
# ╠═fa79434d-8e3d-4acf-a1a7-2eb78d1bde35
# ╠═9d47fdb6-8fbf-4165-9bda-c65e65872774
# ╠═e0f3583a-4b3f-4010-bafa-964fa4fa8fe6
# ╠═7339f1ac-1ae6-4b96-9cce-b2750d7e82ca
# ╠═eaeddf29-da61-4651-8bf6-eee575a709ae
# ╠═adb15ef5-1a8e-43fc-bffc-4cb73481ab11
# ╠═623117c3-3c61-496e-9a91-cd497c2d4d01
# ╠═ad345b30-2ed7-4d92-91f8-576784d9fe00
# ╠═c0fab2cb-1477-421c-a7a9-353f886ac102
# ╠═125a939f-e7a5-4db7-b94d-fdd143b7369b
# ╠═6cd6ff60-2b42-4309-8832-35f9b8098827
# ╠═b06a9b9c-a1d0-414a-9ccd-320020f1edee
# ╠═ae91d145-fc94-4935-970a-c887df151736
# ╠═1cf811e0-cc58-11eb-2c48-33822012ebf2
# ╠═43c7192c-aa56-4d45-9f47-40964f0e14b6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
const line_regex = r"^(\d+)\s(\d+)"

"""
    readnetwork(filename::String, limit::Number = Inf; fromzero::Bool = false)

Read a space separated file into an (undirected) Graph in an efficient way.

# Arguments
* `limit::Number`: Maximum number of lines to read (good for large files)
* `fromzero::Bool`: Whether vertices are counted from zero; if so, correct accordingly
"""
function readnetwork(filename::String, limit::Number = Inf; fromzero::Bool = false)
    graph = Graph()
    vertices = 0
    const correction = Int(fromzero)

    # using grep to exclude comment lines
    open(filename) do file
        for (l, line) in enumerate(eachline(file))
            if l > limit
                break
            end

            if (m = match(line_regex, line)) !== nothing
                raw_v1, raw_v2 = m.captures
                v1 = parse(Int, raw_v1) + correction
                v2 = parse(Int, raw_v2) + correction

                new_vertices = max(v1, v2)
                if new_vertices > vertices
                    @assert add_vertices!(graph, new_vertices - vertices)
                    vertices = new_vertices
                end

                # we explicitely ignore inverse directions, if there are any
                @assert has_edge(graph, v1, v2) || add_edge!(graph, v1, v2)
            end
        end
    end

    return graph
end

logistic(x) = 1 / (1 + exp(-x))
softmax(xs) = let e = exp(xs); e / sum(e) end

"""
    annealing(loss::Function, initial_parameter, initial_temperature, temp_decay, stepsize, stop)

Find a parameter setting minimizing `loss`, using simulated annealing.
"""
function annealing(loss::Function, initial_parameter, initial_temperature,
                   temp_decay, stepsize, stop;
                   debug_callback = (;kwargs...) -> return)
    const k = length(initial_parameter)

    steps = 1
    T = initial_temperature
    x_old = initial_parameter
    f_old = loss(initial_parameter)

    while !stop(steps, f_old)
        direction = normalize(randn(k)) # uniform on hypersphere
        x_new = x_old + direction * stepsize
        f_new = loss(x_new)

        if (updated = rand() <= min(1, exp(-(f_new - f_old) / T)))
            x_old = x_new
            f_old = f_new
        end

        debug_callback(steps = steps, T = T, updated = updated, direction = direction,
                       x_old = x_old, f_old = f_old,
                       f_new = f_new, x_new = x_new)

        T = max(temp_decay * T, eps())
        steps += 1
    end

    return x_old
end


text:
export @D, @DD, @GRAD, @pot

# ===========================================================================
#     implement some fun little macros for easier access
#     to the potentials
# ===========================================================================

# Julia 0.4 version:
# call(pp::Potential, varargs...) = evaluate(pp, varargs...)
# call(pp::Potential, ::Type{Val{:D}}, varargs...) = evaluate_d(pp, varargs...)
# call(pp::Potential, ::Type{Val{:DD}}, varargs...) = evaluate_dd(pp, varargs...)
# call(pp::Potential, ::Type{Val{:GRAD}}, varargs...) = grad(pp, varargs...)

# unfortunately, in 0.5 `call` doesn't take an abstract argument anymore,
# which means that we need to specify for every potential how to
# create this syntactic sugar. This is what `@pot` is for.

"""
Annotate a type with `@pot` to setup the syntax sugar
for `evaluate, evaluate_d, evaluate_dd, grad`.

## Usage:

For example, the declaration
```julia
@pot type LennardJones <: PairPotential
   r0::Float64
end
"documentation for `LennardJones`"
LennardJones
```
creates the following aliases:
```julia
lj = LennardJones(1.0)
lj(args...) = evaluate(lj, args...)
@D lj(args...) = evaluate_d(lj, args...)
@DD lj(args...) = evaluate_dd(lj, args...)
@GRAD lj(args...) = grad(lj, args...)
```

Usage of `@pot` is not restricted to pair potentials, but can be applied to
*any* type.
"""
macro pot(fsig)
   @assert fsig.head == :type
   tname, tparams = t_info(fsig.args[2])
   tname = esc(tname)
   for n = 1:length(tparams)
      tparams[n] = esc(tparams[n])
   end
   sym = esc(:x)
   return quote
      $(esc(fsig))
      ($sym::$tname){$(tparams...)}(args...) = evaluate($sym, args...)
      ($sym::$tname){$(tparams...)}(::Type{Val{:D}}, args...) = evaluate_d($sym, args...)
      ($sym::$tname){$(tparams...)}(::Type{Val{:DD}}, args...) = evaluate_dd($sym, args...)
      ($sym::$tname){$(tparams...)}(::Type{Val{:GRAD}}, args...) = grad($sym, args...)
   end
end

# t_info extracts type name as symbol and type parameters as an array
t_info(ex::Symbol) = (ex, tuple())
t_info(ex::Expr) = ex.head == :(<:) ? t_info(ex.args[1]) : (ex, ex.args[2:end])

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

# next create macros that translate
"""
`@D`: Use to evaluate the derivative of a potential. E.g., to compute the
Lennard-Jones potential,
```julia
lj = LennardJones()
r = 1.0 + rand(10)
ϕ = lj(r)
ϕ' = @D lj(r)
```
see also `@DD`.
"""
macro D(fsig::Expr)
    @assert fsig.head == :call
    insert!(fsig.args, 2, Val{:D})
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    return fsig
end

"`@DD` : analogous to `@D`"
macro DD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:DD})
    return fsig
end

"`@GRAD` : analogous to `@D`, but escapes to `grad`"
macro GRAD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:GRAD})
    return fsig
end

# ==================================
#   Basic Potential Arithmetic

@pot type SumPot{P1, P2} <: PairPotential
   p1::P1
   p2::P2
end
"sum of two pair potentials"
SumPot
import Base.+
+(p1::PairPotential, p2::PairPotential) = SumPot(p1, p2)
evaluate(p::SumPot, r) = p.p1(r) + p.p2(r)
evaluate_d(p::SumPot, r) = (@D p.p1(r)) + (@D p.p2(r))
cutoff(p::SumPot) = max(cutoff(p.p1), cutoff(p.p2))
function Base.print(io::Base.IO, p::SumPot)
   print(io, p.p1)
   print(io, " + ")
   print(io, p.p2)
end

@pot type ProdPot{P1, P2} <: PairPotential
   p1::P1
   p2::P2
end
"product of two pair potentials"
ProdPot
import Base.*
*(p1::PairPotential, p2::PairPotential) = ProdPot(p1, p2)
evaluate(p::ProdPot, r) = p.p1(r) * p.p2(r)
evaluate_d(p::ProdPot, r) = (p.p1(r) * (@D p.p2(r)) + (@D p.p1(r)) * p.p2(r))
cutoff(p::ProdPot) = min(cutoff(p.p1), cutoff(p.p2))
function Base.print(io::Base.IO, p::ProdPot)
   print(io, p.p1)
   print(io, " * ")
   print(io, p.p2)
end

# expand usage of prodpot to be useful for TightBinding.jl
# TODO: make sure this is consistent and expand this to other things!
#       basically, we want to allow that
#       a pair potential can depend on direction as well!
#       in this case, @D is already the gradient and @GRAD remaind undefined?
evaluate{P1,P2}(p::ProdPot{P1,P2}, r, R) = p.p1(r, R) * p.p2(r, R)
evaluate_d(p::ProdPot, r, R) = p.p1(r,R) * (@D p.p2(r,R)) + (@D p.p1(r,R)) * p.p2(r,R)


text:
@coeff function K_tau(reactor::Reactor)
  5.818e-3 * epsilon^2 * kappa * G / tau_P
end


text:
using Test
using ApproxManifoldProducts
using Manifolds

##

@testset "test getManifoldPartial on Euclidean(N)" begin

##

M = TranslationGroup(3)

@test getManifoldPartial(M, [1;2;3])[1] == TranslationGroup(3)
@test getManifoldPartial(M, [2;3])[1] == TranslationGroup(2)

@test getManifoldPartial(M, [1;2;3], zeros(3))[1] == TranslationGroup(3)
@test isapprox( getManifoldPartial(M, [1;2;3], zeros(3))[2], zeros(3) )

@test getManifoldPartial(M, [2;3], zeros(3))[1] == TranslationGroup(2)
@test isapprox( getManifoldPartial(M, [2;3], zeros(3))[2], zeros(2))

##
end

@testset "test getManifoldPartial on Circle()" begin

##

M = Circle()

@test getManifoldPartial(M, [1])[1] == Circle()
@test_throws ErrorException getManifoldPartial(M, [2;])

@test getManifoldPartial(M, [1], [0;])[1] == Circle()
@test getManifoldPartial(M, [1], [0;])[2] == [0]

##
end

@testset "test getManifoldPartial on Rotations(2)" begin

##

M = Rotations(2)

@test getManifoldPartial(M, [1])[1] == Rotations(2)
@test_throws ErrorException getManifoldPartial(M, [2;])

@test getManifoldPartial(M, [1], [1 0; 0 1])[1] == Rotations(2)
@test getManifoldPartial(M, [1], [1 0; 0 1])[2] == [1 0; 0 1]

##
end

@testset "test getManifoldPartial on SpecialEuclidean(2)" begin

##

M = SpecialEuclidean(2)

@test getManifoldPartial(M, [1;2;3])[1] == SpecialEuclidean(2)

@test getManifoldPartial(M, [1;])[1] == TranslationGroup(1)
@test getManifoldPartial(M, [2;])[1] == TranslationGroup(1)
@test getManifoldPartial(M, [1;2])[1] == TranslationGroup(2)

@test getManifoldPartial(M, [3;])[1] == SpecialOrthogonal(2)

@test getManifoldPartial(M, [1;3;])[1] == ProductManifold(TranslationGroup(1), SpecialOrthogonal(2))

repr = ProductRepr([0.0; 0], [1 0; 0 1.0])

@test getManifoldPartial(M, [1;2;3], repr)[1] == SpecialEuclidean(2)
@test getManifoldPartial(M, [1;2;3], repr)[2] == repr

@test getManifoldPartial(M, [1;], repr)[1] == TranslationGroup(1)
@test getManifoldPartial(M, [1;], repr)[2] == [0.0;]

@test getManifoldPartial(M, [2;], repr)[1] == TranslationGroup(1)
@test getManifoldPartial(M, [2;], repr)[2] == [0.0;]

@test getManifoldPartial(M, [1;2], repr)[1] == TranslationGroup(2)
@test getManifoldPartial(M, [1;2], repr)[2] == [0.0;0]

@test getManifoldPartial(M, [3;], repr)[1] == SpecialOrthogonal(2)
@test getManifoldPartial(M, [3;], repr)[2] == repr.parts[2]

@test getManifoldPartial(M, [1;3;], repr)[1] == ProductManifold(TranslationGroup(1), SpecialOrthogonal(2))
r_repr = getManifoldPartial(M, [1;3;], repr)[2]
@test r_repr isa ProductRepr
@test r_repr.parts[1] == [0.0;]
@test r_repr.parts[2] == [1 0; 0 1.0]

##
end

@testset "Reminder, getManifoldPartial on Sphere(2) [TBD]" begin

##

@error "Must fix Sphere(2) partial test"
@test_broken false

##
end

@testset "test getPoints under partial with representation" begin

##

N = 100
M = SpecialEuclidean(2)
u0 = ProductRepr([0.0; 0], [1 0; 0 1.0])

pts = [exp(M, u0, hat(M, u0, [10 .+ randn(2);randn()])) for i in 1:N]

P = manikde!(M, pts)

P12 = marginal(P, [1;2])

p12 = getPoints(P12)

@test length(p12) == N
@test length(p12[1]) == 2
@test_broken P12.manifold isa TranslationGroup(2)

##
end

@testset "test marginal of marginal (partial) helper" begin
##

M = TranslationGroup(3)
pts = [randn(3) for _ in 1:75]

X = manikde!(M, pts, partial=[1;3])

X_ = marginal(X, [3])

ps3 = getPoints(X_)

for (i,pt) in enumerate(pts)
  @test isapprox(ps3[i][1], pt[3])
end

try
  M = TranslationGroup(4)
  # check the constructor when only a few points are available
  X = manikde(M, pts, partial=[1;3;4])
catch
  @test_broken false
end

##
end

@testset "test replace overloads full and partial/marginal" begin
##

N = 10
M = TranslationGroup(3)
pts0 = [zeros(3) for _ in 1:N]
X0 = manikde!(M, pts0, bw=zeros(3))

pts = [randn(3) for _ in 1:N]
X = manikde!(M, pts)

##

X_ = replace(X0, X)
@test isapprox(X_, X)

##

X = manikde!(M, pts, partial=[1;3])
X_ = replace(X0, X)

# check metadata
@test  isapprox( getBW(X_, false)[[1;3],1], getBW(X,  false)[[1;3],1] )
@test !isapprox( getBW(X_, false)[[1;3],1], getBW(X0, false)[[1;3],1] )

@test  isapprox( X_.infoPerCoord[[1;3]], X.infoPerCoord[[1;3]] )

@test !isPartial(X_)

# check points
x0 = getPoints(X0)
x  = getPoints(X, false)
x_ = getPoints(X_)
for (i,pt) in enumerate(x0)
  # partial of X does not replace
  @test isapprox(pt[2], x_[i][2])
  @test isapprox(x[i][[1;3]], x_[i][[1;3]])
end

## must also test replace for partial into different partial

pts3 = [randn(3) for _ in 1:N]
X3 = manikde!(M, pts3, partial=[3;])
# and replace partial/marginal values
X__ = replace(X, X3)

@test isPartial(X__)
@test X__._partial == [1;3]

x__ = getPoints(X__, false)
for (i,pt) in enumerate(x)
  @test isapprox(pt[1:2], x__[i][1:2])
  @test isapprox(pts3[i][3], x__[i][3])
end

## union of two partials over all dimensions should drop the partial status

pts12 = [randn(3) for _ in 1:N]
X12 = manikde!(M, pts12, partial=[1;2], infoPerCoord=2*ones(3))

X_np = replace(X12, X3)

@test !isPartial(X_np)
@test isapprox( X_np.infoPerCoord, [2;2;1] )

##
end

#


text:
################################################################################
#
#     Hecke.jl : Hecke main file
#
# This file is part of Hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2015-2019 Claus Fieker, Tommy Hofmann
#
################################################################################

module Hecke

################################################################################
#
#  Import/export
#
################################################################################

import Base: show, minimum, rand, prod, copy, rand, ceil, round, size, in,
             powermod, ^, getindex, ==, <, >, +, *, /, \, -, !=, getindex,
             setindex!, transpose, getindex, //, div, floor, max, BigFloat,
             precision, first, StepRange, show, one, zero, inv, iseven, isodd,
             convert, angle, abs2, isless, exponent, isfinite, zeros, rem, min,
             numerator, denominator, exp, maximum, intersect, reduce, sqrt

# To make all exported Nemo functions visible to someone using "using Hecke"
# we have to export everything again
# dong it the "import" route, we can pick & choose...

using Requires

using LinearAlgebra, Markdown, InteractiveUtils, Libdl, Distributed, Printf, SparseArrays, Serialization, Random, Pkg, Test

import AbstractAlgebra

import LinearAlgebra: dot, istriu, nullspace

import SparseArrays: nnz

import Serialization: serialize, deserialize

import Random: rand!

import Nemo

exclude = [:Nemo, :AbstractAlgebra, :RealField, :zz, :qq, :factor, :call,
           :factors, :parseint, :strongequal, :window, :xgcd, :rows, :cols,
           :can_solve]

for i in names(Nemo)
  i in exclude && continue
  eval(Meta.parse("import Nemo." * string(i)))
  eval(Expr(:export, i))
end

import Nemo: acb_struct, Ring, Group, Field, NmodRing, nmod, arf_struct,
             elem_to_mat_row!, elem_from_mat_row, gfp_elem, gfp_mat,
             Zmodn_poly, Zmodn_mat, GaloisField, acb_vec, array, acb_vec_clear

export @vprint, @hassert, @vtime, add_verbose_scope, get_verbose_level,
       set_verbose_level, add_assert_scope, get_assert_level, set_assert_level,
       update, show, StepRange, domain, codomain, image, preimage,
       modord, resultant, next_prime, ispower, number_field, factor

###############################################################################
#
#   Library initialisation
#
###############################################################################

const pkgdir = joinpath(dirname(pathof(Hecke)), "..")

#const libhecke = joinpath(pkgdir, "local", "lib", "libhecke")
#
#const libdir = joinpath(pkgdir, "local", "lib")

global const number_field = NumberField

function MaximalOrder
end

global const maximal_order = MaximalOrder

function __init__()

  if myid() == 1
    println("")
    print("Welcome to \n")
    printstyled("
    _    _           _
   | |  | |         | |
   | |__| | ___  ___| | _____
   |  __  |/ _ \\/ __| |/ / _ \\
   | |  | |  __/ (__|   <  __/
   |_|  |_|\\___|\\___|_|\\_\\___|
    ", color = :red)

    println()
    print("Version")
    printstyled(" $VERSION_NUMBER ", color = :green)
    print("... \n ... which comes with absolutely no warranty whatsoever")
    println()
    println("(c) 2015-2019 by Claus Fieker, Tommy Hofmann and Carlo Sircana")
    println()
  else
    println("Hecke $VERSION_NUMBER ...")
  end

  if inNotebook()  # to make toggle work in IJulia
    display("text/html", "\$\\require{action}\$")
  end

  #if "HOSTNAME" in keys(ENV) && ENV["HOSTNAME"] == "juliabox"
  #  push!(Libdl.DL_LOAD_PATH, "/usr/local/lib")
  #elseif Sys.islinux()
  #  push!(Libdl.DL_LOAD_PATH, libdir)
  #  Libdl.dlopen(libhecke)
  #else
  #  push!(Libdl.DL_LOAD_PATH, libdir)
  #end

  t = create_accessors(AnticNumberField, acb_root_ctx, get_handle())
  global _get_nf_conjugate_data_arb = t[1]
  global _set_nf_conjugate_data_arb = t[2]

  t = create_accessors(AnticNumberField, Dict{Int, acb_roots}, get_handle())
  global _get_nf_conjugate_data_arb_roots = t[1]
  global _set_nf_conjugate_data_arb_roots = t[2]

  t = create_accessors(AnticNumberField,
                       Tuple{Int, nf_elem},
                       get_handle())
  global _get_nf_torsion_units = t[1]
  global _set_nf_torsion_units = t[2]

  t = create_accessors(AnticNumberField, NfOrd, get_handle())

  global _get_maximal_order_of_nf = t[1]
  global _set_maximal_order_of_nf = t[2]

  t = create_accessors(NfOrd, ClassGrpCtx, get_handle())

  global _get_ClassGrpCtx_of_order = t[1]
  global _set_ClassGrpCtx_of_order = t[2]

  t = create_accessors(NfOrd, UnitGrpCtx, get_handle())

  global _get_UnitGrpCtx_of_order = t[1]
  global _set_UnitGrpCtx_of_order = t[2]

  t = create_accessors(NfOrd, Array, get_handle())

  global _get_carlos_units_of_order = t[1]
  global _set_carlos_units_of_order = t[2]

  t = create_accessors(AnticNumberField, roots_ctx, get_handle())

  global _get_roots_ctx_of_nf = t[1]
  global _set_roots_ctx_of_nf = t[2]

  t = create_accessors(AnticNumberField, Array, get_handle())

  global _get_cyclotomic_ext_nf = t[1]
  global _set_cyclotomic_ext_nf = t[2]

  t = create_accessors(AnticNumberField, Array, get_handle())

  global _get_automorphisms_nf = t[1]
  global _set_automorphisms_nf = t[2]

  t = create_accessors(AnticNumberField, NfAbsOrd{AnticNumberField, nf_elem}, get_handle())
  global _get_equation_order_of_nf = t[1]
  global _set_equation_order_of_nf = t[2]

  t = create_accessors(SimpleNumField, NfRelOrd, get_handle())

  global _get_maximal_order_of_nf_rel = t[1]
  global _set_maximal_order_of_nf_rel = t[2]

  t = create_accessors(NfOrd, MapClassGrp, get_handle())

  global _get_picard_group = t[1]
  global _set_picard_group = t[2]

  t = create_accessors(NfOrd, MapUnitGrp, get_handle())

  global _get_unit_group_non_maximal = t[1]
  global _set_unit_group_non_maximal = t[2]

  t = create_accessors(AnticNumberField, FacElemMon{AnticNumberField}, get_handle())

  global _get_fac_elem_mon_of_nf = t[1]
  global _set_fac_elem_mon_of_nf = t[2]

  t = create_accessors(NfRel, Array, get_handle())

  global _get_automorphisms_nf_rel = t[1]
  global _set_automorphisms_nf_rel = t[2]

  t = Hecke.create_accessors(AnticNumberField, Dict{Int, Tuple{qAdicRootCtx, Dict{nf_elem, Any}}}, get_handle())
  global _get_nf_conjugate_data_qAdic = t[1]
  global _set_nf_conjugate_data_qAdic = t[2]

  global R = _RealRing()

  # Stuff for elliptic curves
  # polynomial rings Zx = ZZ[x] and _Zxy = ZZ[x,y]
  # will be removed eventually
  #global _Zx = PolynomialRing(FlintZZ, "_x")[1]
  #global _Zxy = PolynomialRing(_Zx, "_y")[1]
  #global _x = gen(_Zx)
  #global _y = gen(_Zxy)

  global flint_rand_ctx = flint_rand_state()

  @require GAP="c863536a-3901-11e9-33e7-d5cd0df7b904" begin
    include("FieldFactory/fields.jl")
    @require Revise="295af30f-e4ad-537b-8983-00126c2a3abe" begin
      import .Revise
      Revise.track(Hecke, "FieldFactory/fields.jl")
      Revise.track(Hecke, "FieldFactory/abelian_layer.jl")
      Revise.track(Hecke, "FieldFactory/brauer.jl")
      Revise.track(Hecke, "FieldFactory/merge.jl")
      Revise.track(Hecke, "FieldFactory/read_write.jl")
    end
  end
end

function _get_maximal_order(K::AnticNumberField)
  return _get_maximal_order_of_nf(K)
end

function _set_maximal_order(K::AnticNumberField, O)
  _set_maximal_order_of_nf(K, O)
end

function _get_nf_equation_order(K::AnticNumberField)
  return _get_equation_order_of_nf(K)::NfAbsOrd{AnticNumberField, nf_elem}
end

function _set_nf_equation_order(K::AnticNumberField, O)
  _set_equation_order_of_nf(K, O)
  return nothing
end

function conjugate_data_arb(K::AnticNumberField)
  try
    c = _get_nf_conjugate_data_arb(K)::acb_root_ctx
    return c
  catch
    c = acb_root_ctx(K.pol)
    _set_nf_conjugate_data_arb(K, c)
    return c::acb_root_ctx
  end
end

function conjugate_data_arb_roots(K::AnticNumberField, p::Int)
  already_set = false
  local c
  try
    c = _get_nf_conjugate_data_arb_roots(K)::Dict{Int, acb_roots}
    already_set = true
  catch
    c = Dict{Int, acb_roots}()
  end

  if already_set && haskey(c, p)
    return c[p]
  end

  #if p > 2^18
  #  Base.show_backtr(STDOUT, backtr())
  #end
  rootc = conjugate_data_arb(K)
  q = rootc.prec
  while q < p
    refine(rootc)
    q = rootc.prec
  end
  @assert p <= q
  rall = deepcopy(rootc.roots)
  rreal = deepcopy(rootc.real_roots)
  rcomplex = deepcopy(rootc.complex_roots)
  for z in rall
    expand!(z, -p)
  end
  for z in rreal
    expand!(z, -p)
  end
  for z in rcomplex
    expand!(z, -p)
  end
  c[p] = acb_roots(p, rall, rreal, rcomplex)
  if !already_set
    _set_nf_conjugate_data_arb_roots(K, c)
  end
  return c[p]::acb_roots
end

function _torsion_units(K::AnticNumberField)
  try
    c = _get_nf_torsion_units(K)::Tuple{Int, nf_elem}
    return c
  catch
    ord, gen = _torsion_units_gen(K)
    _set_nf_torsion_units(K, (ord, gen))
    return ord, gen
  end
end

function torsion_units(K::AnticNumberField)
  ord, g = _torsion_units(K)
  return nf_elem[g^i for i = 0:ord-1]
end

function torsion_units_gen(K::AnticNumberField)
  ar, g = _torsion_units(K)
  return g
end

################################################################################
#
#  Intermediate backwards compatibility
#
################################################################################

trace(x...) = tr(x...)

#lufact(x...) = lu(x...)
#
#lufact!(x...) = lu!(x...)

Base.adjoint(x) = transpose(x)

################################################################################
#
#  Version number
#
################################################################################

global VERSION_NUMBER = v"0.6.5"

################################################################################
#
#  Verbose printing
#
################################################################################

global hecke = Hecke

global VERBOSE_SCOPE = Symbol[]

global VERBOSE_LOOKUP = Dict{Symbol, Int}()

global VERBOSE_PRINT_INDENT = [ 0 ]

function add_verbose_scope(s::Symbol)
  !(s in VERBOSE_SCOPE) && push!(VERBOSE_SCOPE, s)
  nothing
end

function pushindent()
  a = VERBOSE_PRINT_INDENT[1]
  VERBOSE_PRINT_INDENT[1] = a + 1
  nothing
end

function clearindent()
  VERBOSE_PRINT_INDENT[1] = 0
  nothing
end

function popindent()
  a = VERBOSE_PRINT_INDENT[1]
  VERBOSE_PRINT_INDENT[1] = a - 1
  nothing
end

function _global_indent()
  s = "  "^VERBOSE_PRINT_INDENT[1]
  return s
end

macro vprint(args...)
  if length(args) == 2
    quote
      if get_verbose_level($(args[1])) >= 1
        print(_global_indent())
        print($(esc((args[2]))))
      end
    end
  elseif length(args) == 3
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        print(_global_indent())
        print($(esc((args[3]))))
      end
    end
  end
end

macro v_do(args...)
  if length(args) == 2
    quote
      if get_verbose_level($(esc(args[1]))) >= 1
       $(esc(args[2]))
      end
    end
  elseif length(args) == 3
    quote
      if get_verbose_level($(esc(args[1]))) >= $(esc(args[2]))
        $(esc(args[3]))
      end
    end
  end
end

function set_verbose_level(s::Symbol, l::Int)
  !(s in VERBOSE_SCOPE) && error("Not a valid symbol")
  VERBOSE_LOOKUP[s] = l
end

function get_verbose_level(s::Symbol)
  !(s in VERBOSE_SCOPE) && error("Not a valid symbol")
  return get(VERBOSE_LOOKUP, s, 0)::Int
end

################################################################################
#
#  Assertions
#
################################################################################

global ASSERT_SCOPE = Symbol[]

global ASSERT_LOOKUP = Dict{Symbol, Int}()

function add_assert_scope(s::Symbol)
  !(s in ASSERT_SCOPE) && push!(ASSERT_SCOPE, s)
  nothing
end

function set_assert_level(s::Symbol, l::Int)
  !(s in ASSERT_SCOPE) && error("Not a valid symbol")
  ASSERT_LOOKUP[s] = l
end

function get_assert_level(s::Symbol)
  !(s in ASSERT_SCOPE) && error("Not a valid symbol")
  return get(ASSERT_LOOKUP, s, 0)::Int
end

macro hassert(args...)
  if length(args) == 2
    quote
      if get_assert_level($(args[1])) >= 1
        @assert $(esc(args[2]))
      end
    end
  elseif length(args) == 3
    quote
      if get_assert_level($(args[1])) >= $(args[2])
        @assert $(esc(args[3]))
      end
    end
  end
end

######################################################################
# named printing support
######################################################################
# to use:
# in HeckeMap
#   in the show function, start with @show_name(io, map)
# for other objetcs
#   add @declare_other to the struct
#   add @show_name(io, obj) to show
#   optionally, add @show_special(io, obj) as well
# on creation, or whenever, call set_name!(obj, string)
# @show_name will set on printing if bound in the REPL

macro declare_other()
  esc(quote other::Dict{Symbol, Any} end )
end

function set_name!(G, name::String)
  set_special(G, :name => name)
end

function hasspecial(G)
  if isa(G, Map{<:Any, <:Any, HeckeMap, <:Any})
    if isdefined(G.header, :other)
      return true, G.header.other
    else
      return false, nothing
    end
  end
  if !isdefined(G, :other)
    return false, nothing
  else
    return true, G.other
  end
end

import Nemo: get_special, set_special

function get_special(G, s::Symbol)
  fl, D = hasspecial(G)
  fl && return get(D, s, nothing)
  nothing
end

function set_name!(G)
  s = get_special(G, :name)
  s === nothing || return
  sy = find_name(G)
  sy === nothing && return
  set_name!(G, string(sy))
end

function set_special(G, data::Pair{Symbol, <:Any}...)
  if isa(G, Map{<:Any, <:Any, HeckeMap, <:Any})
    if !isdefined(G.header, :other)
      G.header.other = Dict{Symbol, Any}()
    end
    D = G.header.other
  elseif !isdefined(G, :other)
    D = G.other = Dict{Symbol, Any}()
  else
    D = G.other
  end

  for d in data
    push!(D, d)
  end
end

extra_name(G) = nothing

macro show_name(io, O)
  return :( begin
    local i = $(esc(io))
    local o = $(esc(O))
    s = get_special(o, :name)
    if s === nothing
      sy = find_name(o)
      if sy === nothing
        sy = extra_name(o)
      end
      if sy !== nothing
        s = string(sy)
        set_name!(o, s)
      end
    end
    if get(i, :compact, false) &&
       s !== nothing
      print(i, s)
      return
    end
  end )
end

macro show_special(io, O)
  return :( begin
    local i = $(esc(io))
    local o = $(esc(O))
    s = get_special(o, :show)
    if s !== nothing
      s(i, o)
      return
    end
  end )
end

################################################################################
#
#  Custom test function
#
################################################################################

function test_module(x, new::Bool = true)
   julia_exe = Base.julia_cmd()
   if x == "all"
     test_file = joinpath(pkgdir, "test/runtests.jl")
   else
     test_file = joinpath(pkgdir, "test/$x.jl")
   end

   if new
     cmd = "using Test; using Hecke; include(\"$test_file\");"
     @info("spawning ", `$julia_exe -e \"$cmd\"`)
     run(`$julia_exe -e $cmd`)
   else
     @info("Running tests for $x in same session")
     try
       include(test_file)
     catch e
       if isa(e, LoadError)
         println("You need to do \"using Test\"")
       else
         rethrow(e)
       end
     end
   end
end

################################################################################
#
#  Verbose time printing
#
################################################################################

macro vtime(args...)
  if length(args) == 2
    msg = string(args[2])
    quote
      if get_verbose_level($(args[1])) >= 1
        local t0 = time_ns()
        local val = $(esc(args[2]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[2]))
        val2
      end
    end
  elseif length(args) == 3
    msg = string(args[3])
    quote
      if get_verbose_level($(args[1])) >= $(args[2])
        local t0 = time_ns()
        local val = $(esc(args[3]))
        println((time_ns()-t0)/1e9, " @ ", $msg)
        val
      else
        local val2 = $(esc(args[3]))
        val2
      end
    end
  end
end

#usage
# @vtime_add_ellapsed :ClassGroup 2 clg :saturate  s= hnf(a)
# @vtime_add :ClassGroup 2 clg :saturate  0.5
# -> clg.time[:saturate] +=
function _vtime_add(D::Dict, k::Any, v::Any)
  if haskey(D, k)
    D[k] += v
  else
    D[k] = v
  end
end

macro vtime_add(flag, level, var, key, value)
  quote
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, $(esc(value)))
    end
  end
end

macro vtime_add_elapsed(flag, level, var, key, stmt)
  quote
    tm = @elapsed $(esc(stmt))
    if get_verbose_level($flag) >= $level
      _vtime_add($(esc(var)).time, $key, tm)
    end
  end
end

################################################################################
#
#  Exception types
#
################################################################################

mutable struct LowPrecisionCholesky <: Exception end

Base.showerror(io::IO, e::LowPrecisionCholesky) =
    print(io, """
    Negative diagonal in Cholesky decomposition, probably a precision issue""")

mutable struct LowPrecisionLLL <: Exception end

Base.showerror(io::IO, e::LowPrecisionLLL) =
    print(io, """
    Transformation matrix has too large entries relative to precision in LLL""")

# what is this function doing here?
function checkbounds(a::Int, b::Int) nothing; end;

################################################################################
add_assert_scope(:PID_Test)
set_assert_level(:PID_Test, 0)

################################################################################
#
#  "Submodules"
#
################################################################################

include("HeckeTypes.jl")
include("NfRel/Types.jl")
include("AlgAss/Types.jl")
include("Grp.jl")
include("Map.jl")
include("Misc.jl")
include("GrpAb.jl")
include("LinearAlgebra.jl")
include("NfAbs.jl")
include("NumField.jl")
include("NfOrd.jl")
include("Sparse.jl")
include("BigComplex.jl")
include("conjugates.jl")
include("NfRel.jl")
include("analytic.jl")
include("helper.jl")
include("EllCrv.jl")
include("LargeField.jl")
include("RCF.jl")
include("ModAlgAss.jl")
include("AlgAss.jl")
include("AlgAssAbsOrd.jl")
include("AlgAssRelOrd.jl")
include("LocalField.jl")
include("QuadForm.jl")

################################################################################
#
#  Object overloading for map types
#
################################################################################

for T in subtypes(Map(HeckeMap))
  (M::T)(a) = image(M, a)
end

################################################################################
#
#  Verbose printing
#
################################################################################

@doc Markdown.doc"""
    vshow(A) -> Nothing

Prints all fields of $A$.
"""
function vshow(A)
  for i in fieldnames(typeof(A))
    if isdefined(A, i)
      print("$i: ")
      println(getfield(A, i), "\n")
    else
      println("$i: Not definied")
    end
  end
end

################################################################################
#
#  Element types for parent types
#
################################################################################

# Nemo only provides element_types for parent objects

elem_type(::Type{FacElemMon{T}}) where {T} = FacElem{elem_type(T), T}

elem_type(::Type{Generic.ResRing{T}}) where {T} = Generic.Res{T}

################################################################################
#
#  Aliases
#
################################################################################

hasroot(a...) = ispower(a...)  # catch all... needs revisiting:
                               #hasroot(poly) != ispower(poly)....

################################################################################
#
#  Trace function calls
#
################################################################################

_trace_cache = Dict()

function _trace_call(;cache = _trace_cache, depth = 2, print = false)
  s = Base.stacktrace()[4:4 + depth - 1]
  if !haskey(cache, s)
    if print
      println("\n Trace call ... \n ")
      Base.show_backtrace(stdout, s)
    end
    cache[s] = true
  end
end

function _print_traces(;cache = _trace_cache)
  println(cache)
end

################################################################################
#
# examples
#
################################################################################

export example

function example(s::String)
  Base.include(Main, joinpath(dirname(pathof(Hecke)), "..", "examples", s))
end

function revise(s::String)
  s = joinpath(dirname(pathof(Hecke)), "..", "examples", s)
  Main.Revise.track(Main, s)
end

#same (copyed) as in stdlib/v1.0/InteractiveUtils/src/InteractiveUtils.jl
#difference: names(m, all = true) to also see non-exported variables, aka
# caches...

function varinfo(m::Module=Main, pattern::Regex=r"")
    rows =
        Any[ let value = getfield(m, v)
                 Any[string(v),
                     (value===Base || value===Main || value===Core ? "" : format_bytes(summarysize(value))),
                     summary(value)]
             end
             for v in sort!(names(m, all = true)) if isdefined(m, v) && occursin(pattern, string(v)) ]

    pushfirst!(rows, Any["name", "size", "summary"])

    return Markdown.MD(Any[Markdown.Table(rows, Symbol[:l, :r, :l])])
end
varinfo(pat::Regex) = varinfo(Main, pat)

function print_cache(sym::Array{Any, 1})
  for f in sym;
    #if f[2] isa Array || f[2] isa Dict || f[2] isa IdDict;
    try
      print(f[1], " ", length(f[2]), "\n");
    catch e
    end
  end
end

function print_cache()
  print_cache(find_cache(Nemo))
  print_cache(find_cache(Nemo.Generic))
  print_cache(find_cache(Hecke))
end

function find_cache(M::Module)
  sym = []
  for a in collect(names(M, all = true))
    d = Meta.parse("$M.$a")
      try
        z = eval(d);
        if isa(z, AbstractArray) || isa(z, AbstractDict)
          push!(sym, (d, z))
        end
    catch e
    end
  end
  return sym
end

protect = [:(Hecke.ASSERT_LOOKUP), :(Hecke.VERBOSE_LOOKUP),
           :(Hecke.ASSERT_SCOPE), :(Hecke.VERBOSE_SCOPE),
           :(Hecke._euler_phi_inverse_maximum),
           :(Hecke.odlyzko_bound_grh),
           :(Hecke.nC), :(Hecke.B1), #part of ECM
           :(Hecke.VERBOSE_PRINT_INDENT)]

function clear_cache(sym::Array{Any, 1})
  for f in sym;
    if f[1] in protect
      continue
    end
    try
      empty!(f[2])
    catch e
    end
  end
end

function clear_cache()
  clear_cache(find_cache(Nemo))
  clear_cache(find_cache(Nemo.Generic))
  clear_cache(find_cache(Hecke))
end

@inline __get_rounding_mode() = Base.MPFR.rounding_raw(BigFloat)

end # module


text:
module BioformatsLoader
using JavaCall
using LightXML
using ImageMetadata

export
set_id,
OMEXMLReader,
openbytes,
set_series,
getpixeltype,
open_img,
bf_import,
open_stack,
metadata

include("metadata.jl")

const JServiceFactory = @jimport loci.common.services.ServiceFactory
const JService = @jimport loci.common.services.Service
const JIFormatReader = @jimport loci.formats.IFormatReader
const JImageReader = @jimport loci.formats.ImageReader
const JOMEXMLMetadata = @jimport loci.formats.ome.OMEXMLMetadata
const JMetadataStore = @jimport loci.formats.meta.MetadataStore

const datatypes = [Int8,UInt8,Int16,UInt16,Int32,UInt32,Float32,Float64,Bool]

function create_service(classname::String)
    factory = JServiceFactory(())
    cls = classforname(classname)
    service = jcall(factory, "getInstance", JService, (JClass,), cls)
    return convert(JavaObject{Symbol(classname)}, service)
end

struct OMEXMLReader
    reader::JavaObject
    meta::JavaObject
end

function OMEXMLReader()
    service = create_service("loci.formats.services.OMEXMLService")
    meta = jcall(service, "createOMEXMLMetadata", JOMEXMLMetadata, ())
    reader = convert(JIFormatReader, JImageReader(()))

    jcall(reader, "setMetadataStore", Nothing, (JMetadataStore,), meta)
    OMEXMLReader(reader, meta)
end

function set_id(oxr::OMEXMLReader, fname::AbstractString)
    jcall(oxr.reader, "setId", Nothing, (JString,), string(fname))
end

function openbytes(oxr::OMEXMLReader, index::Int)
    jcall(oxr.reader, "openBytes", Vector{jbyte}, (jint,), index)
end

function set_series(oxr::OMEXMLReader, index::Int)
    jcall(oxr.reader, "setSeries", Nothing, (jint,), index)
end

function getpixeltype(oxr::OMEXMLReader)
    id = jcall(oxr.reader, "getPixelType", jint, ())
    return datatypes[id+1]
end

function islittleendian(oxr::OMEXMLReader)
    jcall(oxr.reader, "isLittleEndian", jboolean, ())
end

function open_img(oxr::OMEXMLReader, index::Int)
    sx = jcall(oxr.reader, "getSizeX", jint, ())
    sy = jcall(oxr.reader, "getSizeY", jint, ())
    arr = reinterpret(getpixeltype(oxr), openbytes(oxr, index))
    if islittleendian == 0
        @. arr = ntoh(arr)
    end
    return reshape(arr, (sx, sy))
end

for F in (:open_stack, :bf_import, :metadata)
    @eval begin
        function $F(fname::String, redirect::IO)
            redirect_stdout(() -> $F(fname), redirect)
        end
        function $F(fname::String, show_output::Bool)
            if show_output
                $F(fname)
            else
                STDOLD = stdout
                local ret
                redirect_stdout()
                try
                    ret = $F(fname)
                finally
                    redirect_stdout(STDOLD)
                end
                ret
            end
        end
    end
end

function open_stack(fname::String)
    oxr = OMEXMLReader()
    set_id(oxr, fname)

    SizeX = jcall(oxr.reader, "getSizeX", jint, ())
    SizeY = jcall(oxr.reader, "getSizeY", jint, ())
    SizeC = jcall(oxr.reader, "getSizeC", jint, ())
    SizeZ = jcall(oxr.reader, "getSizeZ", jint, ())
    SizeT = jcall(oxr.reader, "getSizeT", jint, ())
    order = jcall(oxr.reader, "getDimensionOrder", JString, ())

    arr = openbytes(oxr, 0)

    for i in 1:((SizeC*SizeZ*SizeT)-1)
        arr_nd = openbytes(oxr, i)
        append!(arr,arr_nd)
    end

    if islittleendian == 0
        @. arr = ntoh(arr)
    end
    arr = reinterpret(getpixeltype(oxr), arr)

    #Scikit image order of dimension TZXYC (permutedims)
    if order == "XYCZT"
        im = reshape(arr, (SizeX, SizeY, SizeC, SizeZ, SizeT))
        im = permutedims(im, [5,4,1,2,3])
    elseif order == "XYZTC"
        im = reshape(arr, (SizeX, SizeY, SizeZ, SizeT, SizeC))
        im = permutedims(im, [4,3,1,2,5])
    end
end

function bf_import(fname::String)
    oxr = OMEXMLReader()
    set_id(oxr, fname)
    xml = parse_string(jcall(oxr.meta, "dumpXML", JString, ()))
    images = Array{ImageMeta,1}()

    metalst = get_elements_by_tagname(root(xml), "Image")
    for i = 1:length(metalst)
        set_series(oxr, i-1)
        img = open_img(oxr, 0)
        properties = xml_to_dict(metalst[i])
        push!(images, ImageMeta(img, properties))
    end

    return images
end

function metadata(fname::String)
    oxr = OMEXMLReader()
    set_id(oxr, fname)
    xml = parse_string(jcall(oxr.meta, "dumpXML", JString, ()))
    properties = Dict{String, Any}()
    metalst = get_elements_by_tagname(root(xml), "Image")

    SizeC = jcall(oxr.reader, "getSizeC", jint, ())

    for a = attributes(metalst[1])
        push!(properties, "$(name(a))"=>value(a))
    end
    for a = attributes(find_element(metalst[1], "Pixels"))
        try
            push!(properties, "$(name(a))"=> parse(Float64, value(a)))
        catch ex
            if ex isa ArgumentError
                push!(properties, "$(name(a))"=>  value(a))
            end
        end
    end

    for i in 1:SizeC
        for a in attributes(metalst[1]["Pixels"][1]["Channel"][i])
            if name(a)== "EmissionWavelength"
                push!(properties, "channel_$i"=>value(a))
            end
        end
    end
end

const levels = ["ALL", "DEBUG", "ERROR", "FATAL", "INFO", "OFF", "TRACE", "WARN"]
function enableLogging(level::String)
    level in levels || error("enableLogging: level argument must be one of "*join(levels,", "))
    dt = @jimport(loci.common.DebugTools)
    Bool(jcall(dt,"enableLogging",jboolean,(JString,),level))
end

get_bf_path() = joinpath(dirname(@__DIR__), "deps", "bioformats_package.jar")

function init(;memory=1024::Int,log_level::String="ERROR")
    bfpkg_path = get_bf_path()
    JavaCall.init(["-ea", "-Xmx$(memory)M", "-Djava.class.path=$bfpkg_path"])
    enableLogging(log_level) || @warn "Could not enable logging."
    nothing
end

end


text:
module RetroSignalModel

using Catalyst
using Parameters
using DataFrames
using CSV
using DifferentialEquations
using ModelingToolkit
using Random
using LabelledArrays
using FindSteadyStates
using Distributions
using GLPK
using ProgressMeter
using Distributed
import JuMP as jp

include("utils.jl")
include("settings.jl")
include("models/models.jl")
include("data.jl")
include("simulation.jl")
include("paramSearch.jl")

end


text:
abstract type AbstractFeedbackSystem <: AbstractEnv end

# DelayFDI_Plant_BacksteppingControl_StaticAllocator_ControlSystem
struct DelayFDI_Plant_BacksteppingControl_StaticAllocator_ControlSystem <: AbstractFeedbackSystem
    plant::DelayFDI_Plant
    control_system::BacksteppingControl_StaticAllocator_ControlSystem
end

"""
Multicopter + BacksteppingPositionController
"""
function State(env::DelayFDI_Plant_BacksteppingControl_StaticAllocator_ControlSystem)
    @unpack plant, control_system = env
    @unpack multicopter = plant
    @unpack controller = control_system
    @unpack m, g = multicopter
    return function (; args_multicopter=())
        x_plant = State(plant)(; args_multicopter=args_multicopter)
        pos0 = copy(x_plant.multicopter.p)
        x_cs = State(control_system)(pos0, m, g)
        ComponentArray(plant=x_plant, control_system=x_cs)
    end
end

function Dynamics!(env::DelayFDI_Plant_BacksteppingControl_StaticAllocator_ControlSystem)
    @unpack plant, control_system = env
    effectiveness_matrix_functions = EffectivenessMatrixFunction(plant)
    @unpack Λ̂_func = effectiveness_matrix_functions
    @unpack controller, allocator = control_system
    @unpack multicopter = plant
    @unpack m, J, g = multicopter
    @Loggable function dynamics!(dx, x, p, t; pos_cmd=nothing)
        Λ̂ = Λ̂_func(t)
        @unpack p, v, R, ω = x.plant.multicopter
        @unpack ref_model, Td = x.control_system.controller
        xd, vd, ad, ȧd, äd = ref_model.x_0, ref_model.x_1, ref_model.x_2, ref_model.x_3, ref_model.x_4
        command_info = Command(control_system)(p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g, Λ̂)
        @unpack νd, Ṫd, u_cmd = command_info
        @log νd
        @nested_log :plant Dynamics!(plant)(dx.plant, x.plant, (), t; u=u_cmd)
        @nested_log :control_system Dynamics!(control_system)(dx.control_system, x.control_system, (), t; pos_cmd=pos_cmd, Ṫd=Ṫd)
        nothing
    end
end

# DelayFDI_Plant_BacksteppingControl_AdaptiveAllocator_ControlSystem
struct DelayFDI_Plant_BacksteppingControl_AdaptiveAllocator_ControlSystem <: AbstractFeedbackSystem
    plant::DelayFDI_Plant
    control_system::BacksteppingControl_AdaptiveAllocator_ControlSystem
end

function State(env::DelayFDI_Plant_BacksteppingControl_AdaptiveAllocator_ControlSystem)
    @unpack plant, control_system = env
    @unpack multicopter = plant
    @unpack controller = control_system
    @unpack m, g = multicopter
    return function (; args_multicopter=())
        x_plant = State(plant)(; args_multicopter=args_multicopter)
        pos0 = copy(x_plant.multicopter.p)
        x_cs = State(control_system)(pos0, m, g)
        ComponentArray(plant=x_plant, control_system=x_cs)
    end
end

function Dynamics!(env::DelayFDI_Plant_BacksteppingControl_AdaptiveAllocator_ControlSystem)
    @unpack plant, control_system = env
    effectiveness_matrix_functions = EffectivenessMatrixFunction(plant)
    @unpack Λ̂_func = effectiveness_matrix_functions
    @unpack controller, allocator = control_system
    @unpack multicopter = plant
    @unpack m, J, g = multicopter
    @Loggable function dynamics!(dx, x, p, t; pos_cmd=nothing)
        Λ̂ = Λ̂_func(t)
        @unpack p, v, R, ω = x.plant.multicopter
        @unpack ref_model, Td = x.control_system.controller
        Θ̂ = x.control_system.allocator
        xd, vd, ad, ȧd, äd = ref_model.x_0, ref_model.x_1, ref_model.x_2, ref_model.x_3, ref_model.x_4
        command_info = Command(control_system)(p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g, Λ̂, Θ̂)
        @unpack νd, Ṫd, u_cmd, Θ̂̇ = command_info
        @log νd
        @nested_log :plant Dynamics!(plant)(dx.plant, x.plant, (), t; u=u_cmd)
        @nested_log :control_system Dynamics!(control_system)(dx.control_system, x.control_system, (), t; pos_cmd=pos_cmd, Ṫd=Ṫd, Θ̂̇=Θ̂̇)
        nothing
    end
end