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

using Markdown
using InteractiveUtils

# ╔═╡ b0dff76e-8726-11eb-31c8-ab5878dd9ed4
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ bef082a6-8726-11eb-088f-5ffda6b048e7
md"""# Tuples

This chapter presents one more built-in type, the tuple, and then shows how arrays, dictionaries, and tuples work together. It also introduces a useful feature for variable-length argument arrays, the gather and scatter operators."""

# ╔═╡ d10a4cb8-8726-11eb-27e2-b55de2487e54
md"""## Tuples Are Immutable

A *tuple* is a sequence of values. The values can be of any type, and they are indexed by integers, so in that respect tuples are a lot like arrays. The important difference is that tuples are immutable and that each element can have its own type.

Syntactically, a tuple is a comma-separated list of values:

```julia
julia> t = 'a', 'b', 'c', 'd', 'e'
('a', 'b', 'c', 'd', 'e')
```

Although it is not necessary, it is common to enclose tuples in parentheses:

```julia
julia> t = ('a', 'b', 'c', 'd', 'e')
('a', 'b', 'c', 'd', 'e')
```

To create a tuple with a single element, you have to include a final comma:

```julia
julia> t1 = ('a',)
('a',)
julia> typeof(t1)
Tuple{Char}
```

!!! danger
    A value in parentheses without comma is not a tuple:

    ```julia
    julia> t2 = ('a')
    'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
    julia> typeof(t2)
    Char
    ```

Another way to create a tuple is using the built-in function `tuple`. With no argument, it creates an empty tuple:

```julia
julia> tuple()
()
```

If multiple arguments are provided, the result is a tuple with the given arguments:

```julia
julia> t3 = tuple(1, 'a', pi)
(1, 'a', π = 3.1415926535897...)
```

Because `tuple` is the name of a built-in function, you should avoid using it as a variable name.

Most array operators also work on tuples. The bracket operator indexes an element:

```julia
julia> t = ('a', 'b', 'c', 'd', 'e'); julia> t[1]
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```

And the slice operator selects a range of elements:

```julia
julia> t[2:4]
('b', 'c', 'd')
```

But if you try to modify one of the elements of the tuple, you get an error:

```julia
julia> t[1] = 'A'
ERROR: MethodError: no method matching setindex!(::NTuple{5,Char},::Char, ::Int64)
```

Because tuples are immutable, you can’t modify the elements.

The relational operators work with tuples and other sequences. Julia starts by comparing the first element from each sequence. If they are equal, it goes on to the next elements, and so on until it finds elements that differ. Subsequent elements are not considered (even if they are really big):

```julia
julia>(0, 1, 2) < (0, 3, 4)
true
julia> (0, 1, 2000000) < (0, 3, 4)
true
```
"""

# ╔═╡ c7ffd72e-8727-11eb-0f6d-fb8e9395741c
md"""## Tuple Assignment

It is often useful to swap the values of two variables. With conventional assignments, you have to use a temporary variable. For example, to swap `a` and `b`:

```julia
temp=a
a=b
b=temp
```

This solution is cumbersome; *tuple assignment* is more elegant:

```julia
a, b = b, a
```

The left side is a tuple of variables; the right side is a tuple of expressions. Each value is assigned to its respective variable. All the expressions on the right side are evaluated before any of the assignments.

The number of variables on the left has to be fewer than the number of values on the right:

```julia
julia> (a, b) = (1, 2, 3)
(1, 2, 3)
julia>a, b, c= 1, 2
ERROR: BoundsError: attempt to access (1, 2) at index [3]
```

More generally, the right side can be any kind of sequence (string, array, or tuple). For example, to split an email address into a username and a domain, you could write:

```julia
julia> addr = "julius.caesar@rome"
"julius.caesar@rome"
julia> uname, domain = split(addr, '@')
2-element Array{SubString{String},1}:
 "julius.caesar"
 "rome"
```

The return value from `split` is an array with two elements; the first element is assigned to `uname`, the second to `domain`:

```julia
julia> uname
"julius.caesar"
julia> domain
"rome"
```
"""

# ╔═╡ 574b6fb0-8728-11eb-0eca-51b778a97602
md"""## Tuples as Return Values

Strictly speaking, a function can only return one value, but if the value is a tuple, the effect is the same as returning multiple values. For example, if you want to divide two integers and compute the quotient and remainder, it is inefficient to compute `x ÷ y` and then `x % y`. It is better to compute them both at the same time.
The built-in function divrem takes two arguments and returns a tuple of two values, the quotient and remainder. You can store the result as a tuple:

```julia
julia> t = divrem(7, 3)
(2, 1)
```

Or use tuple assignment to store the elements separately:

```julia
julia> q, r = divrem(7, 3);
julia> @show q r;
q = 2
r = 1
```

Here is an example of a function that returns a tuple:

```julia
minmax(t) = minimum(t), maximum(t)
```

`maximum` and `minimum` are built-in functions that find the largest and smallest ele‐ ments of a sequence. `minmax` computes both and returns a tuple of two values. The built-in function `extrema` is more efficient.
"""

# ╔═╡ a5a354be-8728-11eb-2ab3-4f994c1939c1
md"""## Variable-Length Argument Tuples

Functions can take a variable number of arguments. A parameter name that ends with `...` *gathers* arguments into a tuple. For example, `printall` takes any number of arguments and prints them:

```julia
printall(args...) = println(args)
```

The gather parameter can have any name you like, but `args` is conventional. Here’s how the function works:

```julia
julia> printall(1, 2.0, '3')
(1, 2.0, '3')
```

The complement of gather is *scatter*. If you have a sequence of values and you want to pass it to a function as multiple arguments, you can use the `...` operator. For example, `divrem` takes exactly two arguments; it doesn’t work with a tuple:

```julia
julia> t = (7, 3);

julia> divrem(t)
ERROR: MethodError: no method matching divrem(::Tuple{Int64,Int64})
```

But if you scatter the tuple, it works:

```julia
julia> divrem(t...)
(2, 1)
```

Many of the built-in functions use variable-length argument tuples. For example, `max` and `min` can take any number of arguments:

```julia
julia> max(1, 2, 3)
3
```

But sum does not:

```julia
julia> sum(1, 2, 3)
ERROR: MethodError: no method matching sum(::Int64, ::Int64, ::Int64)
```

In the Julia world, gather is often called “slurp” and scatter “splat.”
"""

# ╔═╡ 19ee0db4-8729-11eb-3ca7-737de22181ff
md"""#### Exercise 12-1

Write a function called `sumall` that takes any number of arguments and returns their sum."""

# ╔═╡ 3d55da4a-8729-11eb-3d2f-37fb30c62793
md"""## Arrays and Tuples

`zip` is a built-in function that takes two or more sequences and returns a collection of tuples where each tuple contains one element from each sequence. The name of the function refers to a zipper, which joins and interleaves two rows of teeth.

This example zips a string and an array:

```julia
julia> s = "abc";

julia> t = [1, 2, 3];
julia> zip(s, t)
Base.Iterators.Zip{Tuple{String,Array{Int64,1}}}(("abc", [1, 2, 3]))
```

The result is a *zip object* that knows how to iterate through the pairs. The most common use of `zip` is in a `for` loop:

```julia
julia> for pair in zip(s, t)
	   	   println(pair)
	   end
('a', 1)
('b', 2)
('c', 3)
```

A zip object is a kind of *iterator*, which is any object that iterates through a sequence. Iterators are similar to arrays in some ways, but unlike arrays, you can’t use an index to select an element from an iterator.

If you want to use array operators and functions, you can use a zip object to make an array:

```julia
julia> collect(zip(s, t))
3-element Array{Tuple{Char,Int64},1}:
 ('a', 1)
 ('b', 2)
 ('c', 3)
```

The result is an array of tuples; in this example, each tuple contains a character from the string and the corresponding element from the array.

If the sequences are not the same length, the result has the length of the shorter one:

```julia
julia> collect(zip("Anne", "Elk"))
3-element Array{Tuple{Char,Char},1}:
 ('A', 'E')
 ('n', 'l')
 ('n', 'k')
```

You can use tuple assignment in a `for` loop to traverse an array of tuples:

```julia
julia> t = [('a', 1), ('b', 2), ('c', 3)];

julia> for (letter, number) in t
           println(number, " ", letter)
       end
1a 2b 3c
```

Each time through the loop, Julia selects the next tuple in the array and assigns the elements to `letter` and `number`. The parentheses around (letter, number) are compulsory.

If you combine `zip`, `for`, and tuple assignment, you get a useful idiom for traversing two (or more) sequences at the same time. For example, `hasmatch` takes two sequences, `t1` and `t2`, and returns true if there is an index `i` such that `t1[i] == t2[i]`:

```julia
function hasmatch(t1, t2)
	for (x, y) in zip(t1, t2)
		if x == y
			return true
		end
	end
	false
end
```

If you need to traverse the elements of a sequence and their indices, you can use the built-in function `enumerate`:

```julia
julia> for (index, element) in enumerate("abc")
           println(index, " ", element)
       end
1 a
2 b
3 c
```

The result from `enumerate` is an enumerate object, which iterates a sequence of pairs; each pair contains an index (starting from 1) and an element from the given sequence.
"""

# ╔═╡ 221a3360-872a-11eb-2ead-8d1358efe435
md"""## Dictionaries and Tuples

Dictionaries can be used as iterators that iterate the key-value pairs. You can use a dictionary in a `for` loop like this:

```julia
julia> d = Dict('a'=>1, 'b'=>2, 'c'=>3);

julia> for (key, value) in d
           println(key, " ", value)
       end
a 1
c 3
b 2
```

As you should expect from a dictionary, the items are in no particular order.

Going in the other direction, you can use an array of tuples to initialize a new dictionary:

```julia
julia> t = [('a', 1), ('c', 3), ('b', 2)];

julia> d = Dict(t)
Dict{Char,Int64} with 3 entries:
  'a' => 1
  'c' => 3
  'b' => 2
```

Combining `Dict` with `zip` yields a concise way to create a dictionary:

```julia
julia> d = Dict(zip("abc", 1:3))
Dict{Char,Int64} with 3 entries:
  'a' => 1
  'c' => 3
  'b' => 2
```

It is common to use tuples as keys in dictionaries. For example, a telephone directory might map from *last-name*, *first-name* pairs to telephone numbers. Assuming that we have defined `last`, `first`, and `number`, we could write:

```julia
directory[last, first] = number
```

The expression in brackets is a tuple. We could use tuple assignment to traverse this
dictionary:

```julia
for ((last, first), number) in directory
	println(first, " ", last, " ", number)
end
```

This loop traverses the key-value pairs in `directory`, which are tuples. It assigns the elements of the key in each tuple to `last` and `first`, and the value to `number`, then prints the name and corresponding telephone number.

There are two ways to represent tuples in a state diagram. The more detailed version shows the indices and elements just as they appear in an array. For example, the tuple `("Cleese", "John")` would appear as in Figure 12-1.
"""

# ╔═╡ ce6beec2-872a-11eb-3b39-5531bfc637cb
Drawing(width=720, height=70) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	rect(x=285, y=10, width=150, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=295, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=315, y1=25, x2=355, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=375, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("Cleese")
	end
	text(x=295, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=315, y1=45, x2=355, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=375, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("John")
	end
end

# ╔═╡ bf393bf4-872b-11eb-0233-cbb084e3ad3f
md"*Figure 12-1. State diagram.*"

# ╔═╡ c580c4b4-872b-11eb-1aaf-b51e08f9f26c
md"""But in a larger diagram you might want to leave out the details. For example, a diagram of the telephone directory might appear as in Figure 12-2."""

# ╔═╡ c608335e-872b-11eb-278e-51296de15fbc
Drawing(width=720, height=150) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	rect(x=175, y=10, width=370, height=130, fill="rgb(242, 242, 242)", stroke="black")
	text(x=355, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Cleese\", \"John\")")
	end
	line(x1=365, y1=25, x2=405, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
	text(x=355, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Chapman\", \"Graham\")")
	end
	line(x1=365, y1=45, x2=405, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
	text(x=355, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Idle\", \"Eric\")")
	end
	line(x1=365, y1=65, x2=405, y2=65, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
	text(x=355, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Gilliam\", \"Terry\")")
	end
	line(x1=365, y1=85, x2=405, y2=85, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
	text(x=355, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Jones\", \"Terry\")")
	end
	line(x1=365, y1=105, x2=405, y2=105, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
	text(x=355, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("(\"Palin\", \"Michael\")")
	end
	line(x1=365, y1=125, x2=405, y2=125, stroke="black", marker_end="url(#arrow)")
	text(x=425, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("08700 100 222")
	end
end

# ╔═╡ 1f4b8a1e-872d-11eb-20b1-014ddea498de
md"*Figure 12-2. State diagram.*"

# ╔═╡ 286fe504-872d-11eb-06d3-e1c2682eff9c
md"""Here the tuples are shown using Julia syntax as a graphical shorthand. The telephone number in the diagram is the complaints line for the BBC, so please don’t call it."""

# ╔═╡ 3260e720-872d-11eb-065c-93bbb15621b9
md"""## Sequences of Sequences

I have focused on arrays of tuples, but almost all of the examples in this chapter also work with arrays of arrays, tuples of tuples, and tuples of arrays. To avoid enumerating the possible combinations, it is sometimes easier to talk about sequences of sequences.

In many contexts, the different kinds of sequences (strings, arrays, and tuples) can be used interchangeably. So how should you choose one over the others?

To start with the obvious, strings are more limited than other sequences because the elements have to be characters. They are also immutable. If you need the ability to change the characters in a string (as opposed to creating a new string), you might want to use an array of characters instead.

Arrays are more common than tuples, mostly because they are mutable. But there are a few cases where you might prefer tuples:

* In some contexts, like a `return` statement, it is syntactically simpler to create a tuple than an array.
* If you are passing a sequence as an argument to a function, using tuples reduces the potential for unexpected behavior due to aliasing.
* For performance reasons. The compiler can specialize on the type.

Because tuples are immutable, they don’t provide functions like `sort!` and `reverse!`, which modify existing arrays. But Julia provides the built-in functions `sort`, which takes an array and returns a new array with the same elements in sorted order, and `reverse`, which takes any sequence and returns a sequence of the same type in reverse order.
"""

# ╔═╡ 6cd09676-872d-11eb-0ffd-3fd772eb7db4
md"""## Debugging

Arrays, dictionaries, and tuples are examples of *data structures*; in this chapter we are starting to see compound data structures, like arrays of tuples, or dictionaries that contain tuples as keys and arrays as values. Compound data structures are useful, but they are prone to what I call *shape errors*; that is, errors caused when a data structure has the wrong type, size, or structure. For example, if you are expecting an array with one integer and I give you a plain old integer (not in an array), it won’t work.

Julia allows you to attach a type to elements of a sequence. How this is done is detailed in Lecture 17. Specifying the type eliminates a lot of shape errors.
"""

# ╔═╡ 96a213d2-872d-11eb-25e4-fdd44a16150d
md"""## Glossary

*tuple*:
An immutable sequence of elements where every element can have its own type.

*tuple assignment*:
An assignment with a sequence on the right side and a tuple of variables on the left. The right side is evaluated and then its elements are assigned to the variables on the left.

*gather*:
The operation of assembling a variable-length argument tuple.

*scatter*:
The operation of treating a sequence as a list of arguments.

*zip object*:
The result of calling a built-in function zip, an object that iterates through a sequence of tuples.

*iterator*:
An object that can iterate through a sequence, but that does not provide array operators and functions.

*data structure*:
A collection of related values, often organized in arrays, dictionaries, tuples, etc.

*shape error*:
An error caused because a value has the wrong shape; that is, the wrong type or size.
"""

# ╔═╡ e7f038d6-872d-11eb-28a3-6307d3e72e85
md"""## Exercises

#### Exercise 12-2

Write a function called `mostfrequent` that takes a string and prints the letters in decreasing order of frequency. Find text samples from several different languages and see how letter frequency varies between languages. Compare your results with the tables at [https://en.wikipedia.org/wiki/Letter_frequencies](https://en.wikipedia.org/wiki/Letter_frequencies).
"""

# ╔═╡ 05a7a222-872e-11eb-3de9-c73ab5da92a0
md"""#### Exercise 12-3

More anagrams!

1. Write a program that reads a word list from a file (see “Reading Word Lists”) and prints all the sets of words that are anagrams. Here is an example of what the output might look like:

   ```julia
   ["deltas", "desalt", "lasted", "salted", "slated", "staled"]
   ["retainers", "ternaries"]
   ["generating", "greatening"]
   ["resmelts", "smelters", "termless"]
   ```

   !!! tip
       You might want to build a dictionary that maps from a collection of letters to an array of words that can be spelled with those letters. The question is, how can you represent the col‐ lection of letters in a way that can be used as a key?

2. Modify the previous program so that it prints the longest array of anagrams first, followed by the second longest, and so on.

3. In Scrabble a “bingo” is when you play all seven tiles in your rack, along with a letter on the board, to form an eight-letter word. What collection of eight letters forms the most possible bingos?
"""

# ╔═╡ 586908a2-872e-11eb-3919-f96ba66ee674
md"""#### Exercise 12-4

Two words form a “metathesis pair” if you can transform one into the other by swap‐ ping two letters; for example, “converse” and “conserve.” Write a program that finds all of the metathesis pairs in *words.txt*.

!!! tip
    Don’t test all pairs of words, and don’t test all possible swaps.

*Credit*: This exercise is inspired by an example at [http://puzzlers.org](http://puzzlers.org).
"""

# ╔═╡ 8a7b57dc-872e-11eb-189c-7f2d7cb2e695
md"""#### Exercise 12-5

Here’s another Car Talk Puzzler:

> What is the longest English word, that remains a valid English word, as you remove its letters one at a time?
>
> Now, letters can be removed from either end, or the middle, but you can’t rearrange any of the letters. Every time you drop a letter, you wind up with another English word. If you do that, you’re eventually going to wind up with one letter and that too is going to be an English word—one that’s found in the dictionary. I want to know what’s the longest word. What’s the word, and how many letters does it have?
>
> I’m going to give you a little modest example: Sprite. Ok? You start off with sprite, you take a letter off, one from the interior of the word, take the r away, and we’re left with the word spite, then we take the e off the end, we’re left with spit, we take the s off, we’re left with pit, it, and I.

Write a program to find all words that can be reduced in this way, and then find the longest one.

!!! tip
    This exercise is a little more challenging than most, so here are some suggestions:

    1. You might want to write a function that takes a word and com‐ putes an array of all the words that can be formed by removing one letter. These are the “children” of the word.

    2. Recursively, a word is reducible if any of its children are reducible. As a base case, you can consider the empty string reducible.

    3. The word list I provided, *words.txt*, doesn’t contain singleletter words. So, you might want to add “I,” “a,” and the empty string.

    4. To improve the performance of your program, you might want to memoize the words that are known to be reducible.
"""

# ╔═╡ Cell order:
# ╟─b0dff76e-8726-11eb-31c8-ab5878dd9ed4
# ╟─bef082a6-8726-11eb-088f-5ffda6b048e7
# ╟─d10a4cb8-8726-11eb-27e2-b55de2487e54
# ╟─c7ffd72e-8727-11eb-0f6d-fb8e9395741c
# ╟─574b6fb0-8728-11eb-0eca-51b778a97602
# ╟─a5a354be-8728-11eb-2ab3-4f994c1939c1
# ╟─19ee0db4-8729-11eb-3ca7-737de22181ff
# ╟─3d55da4a-8729-11eb-3d2f-37fb30c62793
# ╟─221a3360-872a-11eb-2ead-8d1358efe435
# ╟─ce6beec2-872a-11eb-3b39-5531bfc637cb
# ╟─bf393bf4-872b-11eb-0233-cbb084e3ad3f
# ╟─c580c4b4-872b-11eb-1aaf-b51e08f9f26c
# ╟─c608335e-872b-11eb-278e-51296de15fbc
# ╟─1f4b8a1e-872d-11eb-20b1-014ddea498de
# ╟─286fe504-872d-11eb-06d3-e1c2682eff9c
# ╟─3260e720-872d-11eb-065c-93bbb15621b9
# ╟─6cd09676-872d-11eb-0ffd-3fd772eb7db4
# ╟─96a213d2-872d-11eb-25e4-fdd44a16150d
# ╟─e7f038d6-872d-11eb-28a3-6307d3e72e85
# ╟─05a7a222-872e-11eb-3de9-c73ab5da92a0
# ╟─586908a2-872e-11eb-3919-f96ba66ee674
# ╟─8a7b57dc-872e-11eb-189c-7f2d7cb2e695


text:
function delete_checksum_lines(; checksum_file::AbstractString,
                                 stdlib_name::AbstractString)
    lines = readlines(checksum_file)
    f = line -> startswith(strip(line), "$(stdlib_name)-")
    file_contains_matching_lines = any(f, lines)
    filter!(!f, lines)
    rm(checksum_file; force = true, recursive = true)
    open(checksum_file, "w") do io
        for line in lines
            line_stripped = strip(line)
            if !isempty(line_stripped)
                println(io, line_stripped)
            end
        end
    end
    return file_contains_matching_lines
end

ext: jl


text:
module MaxHelpingHandSeekerBarrierColorController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/SeekerBarrierColorController" SeekerBarrierColorController(x::Integer, y::Integer,
    color::String="FFFFFF", particleColor::String="FFFFFF", transparency::Number=0.15, particleTransparency::Number=0.5, persistent::Bool=false, particleDirection::Number=0.0)

const placements = Ahorn.PlacementDict(
    "Seeker Barrier Color Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SeekerBarrierColorController
    )
)

function Ahorn.selection(entity::SeekerBarrierColorController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SeekerBarrierColorController, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.northernLights, -12, -12)

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

utilities for spawning benchmark runs
=#

# reduce printing for worker
Base.eval(Distributed, :(function redirect_worker_output(ident, stream)
    @async while !eof(stream)
        println(readline(stream))
    end
end))

function spawn_step(
    fun::Function,
    fun_name::Symbol,
    time_limit::Real,
    worker::Int,
    )
    @assert nprocs() == 2
    status = :Success
    time_start = time()

    fut = Future()
    @async put!(fut, begin
        try
            remotecall_fetch(fun, worker)
        catch e
            println(fun_name, ": caught error: ", e)
        end
    end)

    killed_proc = false
    while !isready(fut)
        if Sys.free_memory() < free_memory_limit
            killstatus = :KilledMemory
        elseif time() - time_start > time_limit
            killstatus = :KilledTime
        else
            sleep(1)
            continue
        end
        try
            if !isready(fut)
                interrupt()
                sleep(5)
            end
            if !isready(fut) && nprocs() > 1
                w = workers()[end]
                run(`kill -SIGKILL $(remotecall_fetch(getpid, w))`)
                killed_proc = true
                sleep(10)
            end
        catch e
        end
        status = Symbol(fun_name, killstatus)
        println("status: ", status)
        break
    end

    output = (killed_proc || !isready(fut)) ? nothing : fetch(fut)
    if isnothing(output)
        if status == :Success
            status = Symbol(fun_name, :CaughtError)
        end
        output = NamedTuple()
    end
    finalize(fut)
    flush(stdout); flush(stderr)

    return (status, output)
end

function spawn_instance(
    ex_name::String,
    ex_type::Type{<:Examples.ExampleInstanceJuMP},
    compile_inst::Tuple,
    inst_data::Tuple,
    extender::Union{Symbol, Nothing},
    solver::Tuple,
    solve::Bool,
    num_threads::Int,
    )
    worker = addprocs(1, enable_threaded_blas = true,
        exeflags = `--threads $num_threads`)[1]
    @assert nprocs() == 2
    println("loading files")
    @fetchfrom worker begin
        @eval import LinearAlgebra
        LinearAlgebra.BLAS.set_num_threads(num_threads)
        @eval using MosekTools
        include(joinpath(@__DIR__, "../../examples/Examples.jl"))
        @eval using Main.Examples
        flush(stdout); flush(stderr)
        return
    end
    println("running compile instance")
    original_stdout = stdout
    (out_rd, out_wr) = redirect_stdout() # don't print output
    @fetchfrom worker begin
        Examples.run_instance(ex_type, compile_inst, extender, NamedTuple(), solver[2],
            default_options = solver[3], test = false)
        flush(stdout); flush(stderr)
        return
    end
    redirect_stdout(original_stdout)
    close(out_wr)
    println("finished compile instance")

    println("\nsetup model")
    print_memory()

    setup_model_args = (ex_type, inst_data, extender, solver[3], solver[2])
    setup_fun() = @eval begin
        (model, model_stats) = Examples.setup_model($setup_model_args...)
        GC.gc()
        return model_stats
    end
    setup_time = @elapsed (script_status, model_stats) =
        spawn_step(setup_fun, :SetupModel, setup_time_limit, worker)
    setup_killed = (script_status != :Success)
    if setup_killed
        println("setup model failed: $script_status")
        model_stats = NamedTuple()
    end

    if solve && !setup_killed
        println("\nsolve and check")
        print_memory()
        check_fun() = @eval begin
            solve_stats = Examples.solve_check(model, test = false)
            return solve_stats
        end
        check_time = @elapsed (script_status, solve_stats) =
            spawn_step(check_fun, :SolveCheck, check_time_limit, worker)
        check_killed = (script_status != :Success)
        check_killed && println("solve and check failed: $script_status")
    else
        solve_stats = NamedTuple()
        check_time = 0.0
        check_killed = true
    end
    if check_killed
        if script_status == :Success
            @assert !solve
            script_status = :SkippedSolveCheck
        end
        solver_hit_limit = true
    else
        solver_hit_limit = (string(solve_stats.status) == "TimeLimit")
        solver_hit_limit && println("solver hit time limit")
    end

    try
        wait(rmprocs(worker))
    catch e
        @warn("error during process shutdown: ", e)
    end
    @assert nprocs() == 1

    script_status = string(script_status)
    run_perf = (; model_stats..., solve_stats..., setup_time,
        check_time, script_status)
    return (setup_killed, solver_hit_limit, run_perf)
end


text:
struct Point2D{T<:Real}
    x::T
    y::T
end

struct Point3D{T<:Real}
    x::T
    y::T
    z::T
end

"""
    euclidean2DRounded(i1::Int,i2::Int, map::Dict{Int, Point2D})

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xy-plane between the two points rounded to the nearest integer

# Arguments
- `i1::Int`: an index of a point in the map
- `i2::Int`: another index of a point in the map
- `map::Dict{Int, Point2D}`: the map of indexes to points
"""
function euclidean2DRounded(i1::Int,i2::Int, map::Dict{Int, Point2D})
    return euclidean2DRounded(map[i1],map[i2])
end

"""
    euclidean2DRounded(p1::Point2D,p2::Point2D)

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xy-plane between the two points rounded to the nearest integer

# Arguments
- `p1::Point2D`: a point to compare
- `p2::Point2D`: another point to compare
"""
function euclidean2DRounded(p1::Point2D,p2::Point2D)
    return round(sqrt(((p1.x - p2.x)^2)+((p1.y - p2.y)^2)))
end

"""
    euclidean2D(i1::Int,i2::Int, map::Dict{Int, Point2D})

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xy-plane between the two points

# Arguments
- `i1::Int`: an index of a point in the map
- `i2::Int`: another index of a point in the map
- `map::Dict{Int, Point2D}`: the map of indexes to points
"""
function euclidean2D(i1::Int,i2::Int, map::Dict{Int, Point2D})
    return euclidean2D(map[i1],map[i2])
end

"""
    euclidean2D(p1::Point2D,p2::Point2D)

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xy-plane between the two points

# Arguments
- `p1::Point2D`: a point to compare
- `p2::Point2D`: another point to compare
"""
function euclidean2D(p1::Point2D,p2::Point2D)
    return sqrt(((p1.x - p2.x)^2)+((p1.y - p2.y)^2))
end

"""
    euclidean3DRounded(i1::Int,i2::Int, map::Dict{Int, Point3D})

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xyz-plane between the two points rounded to the nearest integer

# Arguments
- `i1::Int`: an index of a point in the map
- `i2::Int`: another index of a point in the map
- `map::Dict{Int, Point3D}`: the map of indexes to points
"""
function euclidean3DRounded(i1::Int,i2::Int, map::Dict{Int, Point3D})
    return euclidean3DRounded(map[i1],map[i2])
end

"""
    euclidean3DRounded(p1::Point3D,p2::Point3D)

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xyz-plane between the two points rounded to the nearest integer

# Arguments
- `p1::Point3D`: a point to compare
- `p2::Point3D`: another point to compare
"""
function euclidean3DRounded(p1::Point3D,p2::Point3D)
    return round(sqrt(((p1.x - p2.x)^2)+((p1.y - p2.y)^2)+((p1.z - p2.z)^2)))
end

"""
    euclidean3D(i1::Int,i2::Int, map::Dict{Int, Point3D})

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xyz-plane between the two points

# Arguments
- `i1::Int`: an index of a point in the map
- `i2::Int`: another index of a point in the map
- `map::Dict{Int, Point3D}`: the map of indexes to points
"""
function euclidean3D(i1::Int,i2::Int, map::Dict{Int, Point3D})
    return euclidean3D(map[i1],map[i2])
end

"""
    euclidean3D(p1::Point3D,p2::Point3D)

Takes in a map of indexes to points and two indexes
Returns the euclidean distance on the xyz-plane between the two points

# Arguments
- `p1::Point3D`: a point to compare
- `p2::Point3D`: another point to compare
"""
function euclidean3D(p1::Point3D,p2::Point3D)
    return sqrt(((p1.x - p2.x)^2)+((p1.y - p2.y)^2)+((p1.z - p2.z)^2))
end


text:
# construct constant uniform distribution (case 1)
function constant_distr()

  # define constant function
  function uniform(x, p)
    return 0.5f0
  end

  # binary samples
  samples = [0, 1]
  return uniform, samples
end

# construct distribution with jump at p_crit (case 2)
function jump_distr(;p_crit=1.0f0)

  # define step function
  function stepfunc(x, p, p_crit)
    if p <= p_crit
      if x == zero(eltype(x))
        return one(eltype(p))
      else
        return zero(eltype(p))
      end
    else
      if x == zero(eltype(x))
        return zero(eltype(p))
      else
        return one(eltype(p))
      end
    end
  end

  # binary samples
  samples = [0, 1]
  return (x, p) -> stepfunc(x, p, p_crit), samples
end

# construct ``tilted'' distribution with jump at p_crit (case 3)
function tilt_distr(p_range; p_crit = 2.0f0)

  # tilt function
  function tilt(p, p_crit)
    a_left = -1.0f0
    b_left = 5.0f0
    a_right = -1.0f0
    b_right = 2.0f0
    if p <= p_crit
      x = a_left*p+b_left
    else
      x = a_right*p+b_right
    end

    return x
  end

  # sample from tilt function
  samples = map(p -> tilt(p, p_crit), p_range)
  function tilt_distr(x ,p, p_range)
    target = samples[findall(x -> x==p, p_range)[1]]
    if target == x
      return one(eltype(p))
    else
      return zero(eltype(p))
    end
  end
  return (x, p) -> tilt_distr(x, p, p_range), samples
end

# boltzmann factor
function boltzmann_factor_energy(energy, T)
  if T == Inf
    return one(eltype(T))
  else
    return exp(-energy/T)
  end
end

# probability associated with energy at given temperature based on Boltzmann distribution
function get_probabilities(p_range, energies)
  probs = zeros(eltype(p_range[1]), length(p_range), size(energies)[1])
  for indx_p in 1:length(p_range)
    p = p_range[indx_p]
    partition_func = zero(eltype(p))
    for i in 1:size(energies)[1]
      energy = energies[i, 1]
      num = energies[i, 2]
      probs[indx_p, i] = num*boltzmann_factor_energy(energy, p)
      partition_func += probs[indx_p, i]
    end
    if partition_func == Inf
      probs[indx_p, :] = zeros(eltype(p_range[1]), length(probs[indx_p, :]))
    else
      probs[indx_p, :] = probs[indx_p, :]./partition_func
    end
  end
  return probs
end

# construct probability distribution based on individual probabilities
function thermal_distr(x, p, probs, p_range)
  if p == Inf
    p_indx = length(p_range)
  else
    p_indx = Int(round((p-p_range[1])/(p_range[2]-p_range[1])))+1
  end

  prob_distr = probs[p_indx, Int(x)]
  return prob_distr
end

# construct probability distribution of energy in ising model based on exact enumeration of all configurations
function ising_exact_distr(energies, p_range)
  probs = get_probabilities(p_range, energies)
  return (x, p) -> thermal_distr(x ,p, probs, p_range)
end

# probability for specific sample at fixed value of tuning parameter
function distr_approx_x_p(full_data, p_range, x, p)
  p_indx = Int(round((p-p_range[1])/(p_range[2]-p_range[1])))+1
  return full_data[p_indx, x]
end

# construct probability distribution based on Monte Carlo samples
function distr_approx(energies, unique_energies, numbers, p_range)
  full_data =zeros(eltype(p_range[1]),(length(p_range),length(unique_energies)))
  for i in 1:length(p_range)
    numtot = zero(eltype(p_range[1]))
    for j in 1:length(unique_energies)
      number_list = numbers[i]
      energy_list = energies[i]
      sample = unique_energies[j]
      indices = findall(x -> x==sample, energy_list)
      if length(indices) == 0
        full_data[i, j] = zero(eltype(sample))
      else
        full_data[i, j] = number_list[indices[1]]
        numtot += number_list[indices[1]]
      end
    end
    full_data[i, :] = full_data[i, :]/numtot
  end

  return (x, p) -> distr_approx_x_p(full_data, p_range, x, p), collect(1:length(unique_energies))
end


text:
using FractalTools
using Debugger

ifs = Sierpinski()
generator = ifs.generator
# generator = randalg_sequential_generator(ifs.ws, ifs.probs)

function indicator_measure_2D(x, x0, ϵ)
    return (1/3) * is_in_Ball_2D(x,x0,ϵ)
end

function deneme_2D(x,k)
    1 * ones(size(x,2))
end

elton_integral_2D(deneme_2D, indicator_measure_2D, generator,(2,), ([0.25,0.25],[0.25,0.25]), 1e-7, 1e7)

# elton_integral(deneme_2D, generator,(1,), 1e-8, 1e7)


text:
using Test, Random, AutoGrad, Knet
if isdefined(AutoGrad, :gradcheck); @eval begin
    using AutoGrad: gradcheck
end; end
if isdefined(AutoGrad, :gcheck); @eval begin
    using AutoGrad: gcheck, @gcheck
end; end


text:
immutable _Transform
    a_1::Cfloat
    a_2::Cfloat
    a_3::Cfloat
    a_4::Cfloat
    a_5::Cfloat
    a_6::Cfloat
    a_7::Cfloat
    a_8::Cfloat
    a_9::Cfloat
end

type Transform
    matrix::Array{Real, 1}
end

const transform_identity = Transform([1, 0, 0, 0, 1, 0, 0, 0, 1])

function to_array(t::_Transform)
    Transform([t.a_1, t.a_2, t.a_3, t.a_4, t.a_5, t.a_6, t.a_7, t.a_8, t.a_9])
end


text:
"""
# `MOON_WEBER_2011`

1-D layered model of the Moon, containing density.

## Reference

Weber, R.C.,, Lin, P.-Y.,, Garnero, E.J., Williams, Q., Lognonné, P., 2011.
Seismic Detection of the Lunar Core. Science, 331, 309–312.
https://doi.org/10.1126/science.1199375
"""
const MOON_WEBER_2011 = SteppedLayeredModel(
    r = 1737.1 .- reverse([  0.0,   1.0,   15.0,   40.0,  238.0,
                           488.0, 738.0, 1257.1, 1407.1, 1497.1]),
    vp      = reverse([1.0, 3.2, 5.5, 7.7, 7.8, 7.6, 8.5, 7.5, 4.1, 4.3]),
    vs      = reverse([0.5, 1.8, 3.2, 4.4, 4.4, 4.4, 4.5, 3.2, 0.0, 2.3]),
    density = reverse([2.6, 2.7, 2.8, 3.3, 3.4, 3.4, 3.4, 3.4, 5.1, 8.0]),
)


text:
# Visualize the Results
function plot_trajectory!(plt, model::AbstractGameModel, traj::Algames.Traj)
    plot!(plt, legend=false, aspect_ratio=:equal)
    N = length(traj)
    for i = 1:model.p
        xi = [Algames.state(traj[k])[model.pz[i][1]] for k=1:N]
        yi = [Algames.state(traj[k])[model.pz[i][2]] for k=1:N]
        plot!(xi, yi, label=false)
        scatter!(xi, yi)
    end
    return nothing
end

# Return state information
function trajectory(model::AbstractGameModel, traj::Algames.Traj)
    N = length(traj)
    x = [[Algames.state(traj[k])[model.pz[i][1]] for k=1:N] for i = 1:model.p]
    y = [[Algames.state(traj[k])[model.pz[i][2]] for k=1:N] for i = 1:model.p]
    s = [[Algames.state(traj[k])[model.pz[i][3]] for k=1:N] for i = 1:model.p]
    v = [[Algames.state(traj[k])[model.pz[i][4]] for k=1:N] for i = 1:model.p]
    return x, y, s, v
end

max_stars_repo_name: lassepe/DecNashPlanning

max_issues_repo_name: lassepe/DecNashPlanning

max_forks_repo_name: lassepe/DecNashPlanning


text:
module SimpleIPC

export
 ipc_connect,
 ipc_eval,
 ipc_listen,
 ipc_run,
 ipc_set

import Serialization
import Sockets

macro async_showerr(ex)
 quote
  t = @async while true
   try
    eval($(esc(ex)))
   catch err
    err isa EOFError && break
    bt = catch_backtrace()
    println("Async error: ")
    showerror(stderr, err, bt)
   end
  end
 end
end

function ipc_listen(addr::Sockets.InetAddr)
 server = Sockets.listen(addr)
 @async while true
  socket = Sockets.accept(server)
  Sockets.nagle(socket, false)
  @async_showerr begin
   b = read(socket, UInt8)
   if b == 0x01
    Main.eval(Serialization.deserialize(socket))
   elseif b == 0x02
    name = Serialization.deserialize(socket)
    val = Serialization.deserialize(socket)
    @eval(Main, $name = $val)
   elseif b == 0x03
    fun = Serialization.deserialize(socket)
    args = Serialization.deserialize(socket)
    @eval(Main, $fun($args...))
   end
  end
 end
 server
end
ipc_listen(host::Sockets.IPAddr, port::Integer) = ipc_listen(Sockets.InetAddr(host, port))
ipc_listen(host::AbstractString, port::Integer) = ipc_listen(Sockets.InetAddr(Sockets.getaddrinfo(host), port))
ipc_listen(port::Integer) = ipc_listen(Sockets.localhost, port)

function ipc_connect(addr::Sockets.InetAddr)
 socket = Sockets.connect(addr)
 Sockets.nagle(socket, false)
 socket
end
ipc_connect(host::Sockets.IPAddr, port::Integer) = ipc_connect(Sockets.InetAddr(host, port))
ipc_connect(host::AbstractString, port::Integer) = ipc_connect(Sockets.InetAddr(Sockets.getaddrinfo(host), port))
ipc_connect(port::Integer) = ipc_connect(Sockets.localhost, port)

function ipc_eval(h::Sockets.TCPSocket, x::Expr)
 write(h, 0x01)
 Serialization.serialize(h, x)
 flush(h)
end

function ipc_set(h::Sockets.TCPSocket, name::Symbol, val)
 write(h, 0x02)
 Serialization.serialize(h, name)
 Serialization.serialize(h, val)
 flush(h)
end

function ipc_run(h::Sockets.TCPSocket, fun::Union{Expr,Symbol}, args::Tuple)
 write(h, 0x03)
 Serialization.serialize(h, fun)
 Serialization.serialize(h, args)
 flush(h)
end

end # module


text:
function log_skip_gram(vm::VectorModel, w::T1, s::T2, v::T1) where
	{T1 <: Integer, T2 <: Integer}
	code = view(vm.code, :, v)
	path = view(vm.path, :, v)
	return ccall(_c_skip_gram, Float32,
		(Ptr{Float32}, Ptr{Float32},
			Int,
			Ptr{Int32}, Ptr{Int8}, Int),
		view(vm.In, :, s, w), vm.Out,
			M(vm),
			path, code, size(vm.code, 1))
end

export log_skip_gram


text:
module DensityConsistency

using ExtractMacro, LinearAlgebra, SparseArrays, Random
export DCState, density_consistency

include("types.jl")
include("closures.jl")
include("dc_main.jl")

end # end module

max_stars_repo_name: giovact/DensityConsistency

max_issues_repo_name: giovact/DensityConsistency

max_forks_repo_name: giovact/DensityConsistency


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module Materials

using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters

abstract type AbstractMaterial end
abstract type AbstractMaterialState end

@generated function Base.:+(state::T, dstate::T) where {T <: AbstractMaterialState}
   expr = [:(state.$p + dstate.$p) for p in fieldnames(T)]
   return :(T($(expr...)))
end

export AbstractMaterial, AbstractMaterialState

function integrate_material!(material::M) where {M<:AbstractMaterial}
    error("One needs to define how to integrate material $M!")
end

function update_material!(material::M) where {M <: AbstractMaterial}
    material.drivers += material.ddrivers
    material.parameters += material.dparameters
    material.variables = material.variables_new
    reset_material!(material)
    return nothing
end

function reset_material!(material::M) where {M <: AbstractMaterial}
    material.ddrivers = typeof(material.ddrivers)()
    material.dparameters = typeof(material.dparameters)()
    material.variables_new = typeof(material.variables_new)()
    return nothing
end

export integrate_material!, update_material!, reset_material!

function isotropic_elasticity_tensor(lambda, mu)
    delta(i,j) = i==j ? 1.0 : 0.0
    g(i,j,k,l) = lambda*delta(i,j)*delta(k,l) + mu*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k))
    jacobian = SymmetricTensor{4, 3, Float64}(g)
    return jacobian
end

include("idealplastic.jl")
export IdealPlastic, IdealPlasticDriverState, IdealPlasticParameterState, IdealPlasticVariableState

include("chaboche.jl")
export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState

# include("viscoplastic.jl")
# export ViscoPlastic

include("uniaxial_increment.jl")
export uniaxial_increment!

include("biaxial_increment.jl")
export biaxial_increment!

include("stress_driven_uniaxial_increment.jl")
export stress_driven_uniaxial_increment!

include("DSA.jl")
export DSA, DSADriverState, DSAParameterState, DSAVariableState

end


text:
"""
PrivateKey(𝑒) represents an S256Point determined by 𝑃 = 𝑒G,
where 𝑒 is an integer and G the scep256k1 generator point.
"""
struct PrivateKey
    𝑒::Integer
    𝑃::AbstractPoint
    PrivateKey(𝑒) = new(𝑒, 𝑒 * G)
end

"""
Returns a Signature for a given PrivateKey and data 𝑧
pksign(pk::PrivateKey, 𝑧::Integer) -> Signature
"""
function pksign(pk::PrivateKey, 𝑧::Integer)
    𝑘 = rand(big.(0:N))
    𝑟 = (𝑘 * G).𝑥.𝑛
    𝑘⁻¹ = powermod(𝑘, N - 2, N)
    𝑠 = mod((𝑧 + 𝑟 * pk.𝑒) * 𝑘⁻¹, N)
    if 𝑠 > N / 2
        𝑠 = N - 𝑠
    end
    return Signature(𝑟, 𝑠)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function tau_int_from_atc end

function bat_integrated_autocorr_len_impl(v::AbstractVector{<:Real}, algorithm::AutocorLenAlgorithm)
    atc = fft_autocor(v)
    tau_int_est = tau_int_from_atc(atc, algorithm)
    (result = tau_int_est,)
end

function bat_integrated_autocorr_len_impl(v::AbstractVectorOfSimilarVectors{<:Real}, algorithm::AutocorLenAlgorithm)
    atc = fft_autocor(v)
    flat_atc = flatview(atc)
    tau_int_est = map(axes(flat_atc, 1)) do i
        tau_int_from_atc(view(flat_atc, i, :), algorithm)
    end
    (result = tau_int_est,)
end

"""
    GeyerAutocorLen <: AutocorLenAlgorithm

Integrated autocorrelation length estimation based on Geyer’s initial monotone sequence criterion

See [C. J. Geyer, "Praktical Markov Chain Monte Carlo" (1992)](https://projecteuclid.org/download/pdf_1/euclid.ss/1177011137)
and [C. J. Geyer, "Introduction to Markov Chain Monte Carlo" (2011)](https://www.semanticscholar.org/paper/1-Introduction-to-Markov-Chain-Monte-Carlo-Geyer/21a92825dcec23c743e77451ff5b5ee6b1091651).

Constructor:

```julia
GeyerAutocorLen()
```

The same algorithm is used by
[STAN (v2.21)](https://mc-stan.org/docs/2_21/reference-manual/effective-sample-size-section.html#estimation-of-effective-sample-size)
and (MCMCChains.jl (v3.0, function `ess_rhat`))[https://github.com/TuringLang/MCMCChains.jl/blob/v4.0.0/src/ess.jl#L288].
"""
struct GeyerAutocorLen <: AutocorLenAlgorithm
end

export GeyerAutocorLen

function tau_int_from_atc(atc::AbstractVector{<:Real}, algorithm::GeyerAutocorLen)
    s = zero(eltype(atc))
    Γ_min = eltype(atc)(Inf)

    i = firstindex(atc)
    while i < lastindex(atc) - 1
        Γ = min(atc[i] + atc[i+1], Γ_min)
        if Γ >= 0
            s = s + Γ
            Γ_min = Γ
        else
            break
        end
        i = i + 2
    end

    -1 + 2 * s
end

"""
    SokalAutocorLen <: AutocorLenAlgorithm

Integrated autocorrelation length estimation based on the automated windowing
procedure descibed in
[A. D. Sokal, "Monte Carlo Methods in Statistical Mechanics" (1996)](https://pdfs.semanticscholar.org/0bfe/9e3db30605fe2d4d26e1a288a5e2997e7225.pdf)

Constructor:

```julia
SokalAutocorLen(;c::Integer = 5)
```

* `c`: Step size for window search.

Same procedure is used by the emcee Python package (v3.0).
"""
@with_kw struct SokalAutocorLen <: AutocorLenAlgorithm
    c::Int = 5
end

export SokalAutocorLen

function tau_int_from_atc(atc::AbstractVector{<:Real}, algorithm::SokalAutocorLen)
    c = algorithm.c
    idxs = eachindex(atc)
    idx1 = first(idxs)

    tau_int::eltype(atc) = -1
    for M in idxs
        tau_int += 2 * atc[M]
        if M - idx1 >= c * tau_int
            break
        end
    end

    tau_int
end

function bat_eff_sample_size_impl(v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}}, algorithm::AutocorLenAlgorithm)
    tau_int = bat_integrated_autocorr_len(v, algorithm).result
    n = length(eachindex(v))
    ess = min.(n, n./ tau_int)
    (result = ess,)
end

function bat_eff_sample_size_impl(smpls::DensitySampleVector, algorithm::AutocorLenAlgorithm)
    vs = varshape(smpls)
    unshaped_smpls = unshaped.(smpls)
    n = length(eachindex(unshaped_smpls))

    W = unshaped_smpls.weight
    w0 = first(W)

    unshaped_ess = if all(w -> w ≈ w0, W)
        bat_eff_sample_size_impl(unshaped_smpls.v, algorithm).result
    else
        # If weights not uniform, resample to get unweighted samples. Kish's
        # approximation of ESS for weighted samples is often not good enough.

        # Empirical resampling factor:
        resampling_factor = min(mean(W .^ 2) / mean(W)^2, 10)

        n_resample = round(Int, n * resampling_factor)

        # RNG seed for resampling should be the same for the same samples:
        rng_seed = trunc(UInt64, mean(W) * n)
        rng = Philox4x((0x0, rng_seed))::Philox4x{UInt64,10}

        unweighted_smpls = bat_sample(rng, unshaped_smpls, OrderedResampling(nsamples = n_resample)).result
        resampled_ess = bat_eff_sample_size_impl(unweighted_smpls.v, algorithm).result
        min.(n, resampled_ess)
    end

    result_vs = replace_const_shapes(s::ConstValueShape -> ConstValueShape(Fill(n, size(shape.value)...)), vs)
    ess = result_vs(unshaped_ess)

    (result = ess,)
end


text:
using PyCall
using Conda

if PyCall.conda
    Conda.add_channel("conda-forge")
    Conda.add("pydrive")
else
    try
        pyimport("pydrive")
        # See if it works already
    catch ee
        typeof(ee) <: PyCall.PyError || rethrow(ee)
        error("""
Python pydrive not installed
Please either:
 - Rebuild PyCall to use Conda, by running in the julia REPL:
    - `ENV["PYTHON"]=""; Pkg.build("PyCall"); Pkg.build("PyDrive")`
 - Or install the python binding yourself, eg by running pip
    - `pip install pydrive`
    - then rebuilding PyDrive.jl via `Pkg.build("PyDrive")` in the julia REPL
    - make sure you run the right pip, for the instance of python that PyCall is looking at.
""")
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# This is Agha's example 3.2.2, a recursive factorial
#
using Actors
import Actors: spawn, newLink

# implement the behaviors
function rec_factorial(f::Request)
    if f.x == 0
        send(f.from, Response(1))
    else
        c = spawn(Bhv(rec_customer, f.x, f.from))
        send(self(), Request(f.x - 1, c))
    end
end

function rec_customer(n::Integer, u::Link, k::Response)
    send(u, Response(n * k.y))
    stop()
end

# setup factorial actor and response link
F = spawn(Bhv(rec_factorial))
resp = newLink()

for i ∈ 0:5:50      # send and receive loop
    send(F, Request(big(i), resp))
    println(receive(resp).y)
end

for i ∈ 0:5:50      # send all requests in one loop
    send(F, Request(big(i), resp))
end
for i ∈ 0:10      # receive all results
    println(receive(resp).y)
end


text:
export Weighted
struct Weighted{T}
    val :: T
    logWeight :: Float64
end


text:
# TODO over polytope https://hal.archives-ouvertes.fr/hal-00320882/document and https://math.nist.gov/mcsd/Seminars/2012/2012-09-19-DeLoera-presentation.pdf

module AlgebraicIntegration

using LinearAlgebra
using DynamicPolynomials
import GSL: sf_gamma
include("closed_form.jl")

using PolyJuMP
using SumOfSquares
include("sum_of_squares.jl")

end # module


text:
#TermStructure volatility surface

abstract type VarianceSection end
struct SVISection <: VarianceSection
    a::Float64
    b::Float64
    rho::Float64
    s::Float64
    m::Float64
    tte::Float64
    f::Float64
end

function varianceByLogmoneyness(s::SVISection, y::Float64)
    sqrsy = s.s^2 + (y - s.m)^2
    return s.a + s.b * (s.rho * (y - s.m) + sqrt(sqrsy))
end

function varianceSlopeCurvature(s::SVISection, y::Float64)
    ym = (y - s.m)
    sqrsy = s.s * s.s + ym * ym
    sqrtsy = sqrt(sqrsy)
    variance = s.a + s.b * (s.rho * ym + sqrtsy)
    slope = s.b * (ym / sqrtsy + s.rho)
    curvature = s.b * (1 - ym^2 / sqrsy) / sqrtsy
    return variance, slope, curvature
end

function varianceByLogmoneyness(s::SVISection, y::Array{Float64})
    sqrsy = @. s.s^2 + (y - s.m)^2
    return @. s.a + s.b * (s.rho * (y - s.m) + sqrt(sqrsy))
end

function varianceSlopeCurvature(s::SVISection, y::Array{Float64})
    ym = @. (y - s.m)
    sqrsy = @. (s.s^2 + ym^2)
    sqrtsy = sqrt.(sqrsy)
    variance = @. s.a + s.b * (s.rho * ym + sqrtsy)
    slope = @. s.b * (ym / sqrtsy + s.rho)
    curvature = @. s.b * (1 - ym^2 / sqrsy) / sqrtsy
    return variance, slope, curvature
end
struct VarianceSurfaceBySection{T}
    sections::Vector{T}
    expiries::Vector{Float64}
    #    indexMap::Any
end

#abstract type LocalVarianceSurface end

function varianceByLogmoneyness(
    s::VarianceSurfaceBySection{T},
    y::Float64,
    t::Float64,
    indexTime::Int = 0,
) where {T}
    if t <= s.expiries[1]
        return varianceByLogmoneyness(s.sections[1], y)
    elseif t >= s.expiries[end]
        return varianceByLogmoneyness(s.sections[end], y)
    else
        if indexTime == 0
            indexTime = searchsortedlast(s.expiries, t)
        end
        t0 = s.expiries[indexTime]
        var0 = varianceByLogmoneyness(s.sections[indexTime], y)
        t1 = s.expiries[indexTime+1]
        var1 = varianceByLogmoneyness(s.sections[indexTime+1], y)
        #linear interpolation in total variance along same logmoneyness
        v = (var1 * t1 * (t - t0) + var0 * t0 * (t1 - t)) / (t1 - t0)
        return v / t
    end
end

function varianceSlopeCurvature(
    s::VarianceSurfaceBySection{T},
    y::Float64, #could be a vector?
    t::Float64,
    indexTime::Int = 0, #for perf, have a version with t, tIndex
) where {T}
    if t <= s.expiries[1]
        return varianceSlopeCurvature(s.sections[1], y)
    elseif t >= s.expiries[end]
        return varianceSlopeCurvature(s.sections[end], y)
    else
        if indexTime == 0
            indexTime = searchsortedlast(s.expiries, t)
        end
        t0 = s.expiries[indexTime]
        var0, slope0, curvature0 = varianceSlopeCurvature(s.sections[indexTime], y)
        t1 = s.expiries[indexTime+1]
        var1, slope1, curvature1 = varianceSlopeCurvature(s.sections[indexTime+1], y)
        #linear interpolation in total variance along same logmoneyness
        v = (var1 * t1 * (t - t0) + var0 * t0 * (t1 - t)) / (t1 - t0)
        s = (slope1 * t1 * (t - t0) + slope0 * t0 * (t1 - t)) / (t1 - t0)
        c = (curvature1 * t1 * (t - t0) + curvature0 * t0 * (t1 - t)) / (t1 - t0)
        return v / t, s / t, c / t
    end
end

#length(s::VarianceSurfaceBySection{T}) where {T} = 1

function varianceByLogmoneyness(
    s::VarianceSurfaceBySection{T},
    y::Array{Float64,1},
    t::Float64,
    indexTime::Int = 0,
) where {T}
    if t <= s.expiries[1]
        return varianceByLogmoneyness(s.sections[1], y)
    elseif t >= s.expiries[end]
        return varianceByLogmoneyness(s.sections[end], y)
    else
        if indexTime == 0
            indexTime = searchsortedlast(s.expiries, t)
        end
        t0 = s.expiries[indexTime]
        var0 = varianceByLogmoneyness(s.sections[indexTime], y)
        t1 = s.expiries[indexTime+1]
        var1 = varianceByLogmoneyness(s.sections[indexTime+1], y)
        #linear interpolation in total variance along same logmoneyness
        v = @. (var1 * t1 * (t - t0) + var0 * t0 * (t1 - t)) / (t1 - t0)
        @. v /= t
        return v
    end
end

function varianceSlopeCurvature(
    s::VarianceSurfaceBySection{T},
    y::Array{Float64,1}, #could be a vector?
    t::Float64,
    indexTime::Int = 0, #for perf, have a version with t, tIndex
) where {T}
    if t <= s.expiries[1]
        return varianceSlopeCurvature(s.sections[1], y)
    elseif t >= s.expiries[end]
        return varianceSlopeCurvature(s.sections[end], y)
    else
        if indexTime == 0
            indexTime = searchsortedlast(s.expiries, t)
        end
        t0 = s.expiries[indexTime]
        var0, slope0, curvature0 = varianceSlopeCurvature(s.sections[indexTime], y)
        t1 = s.expiries[indexTime+1]
        var1, slope1, curvature1 = varianceSlopeCurvature(s.sections[indexTime+1], y)
        #linear interpolation in total variance along same logmoneyness
        v = @. (var1 * t1 * (t - t0) + var0 * t0 * (t1 - t)) / (t1 - t0)
        s = @. (slope1 * t1 * (t - t0) + slope0 * t0 * (t1 - t)) / (t1 - t0)
        c = @. (curvature1 * t1 * (t - t0) + curvature0 * t0 * (t1 - t)) / (t1 - t0)
        @. v /= t
        @. s /= t
        @. c /= t
        return v, s, c
    end
end

function localVarianceByLogmoneyness(
    surface::VarianceSurfaceBySection{T}, #another implementation could be based on call surface
    y::Array{Float64,1},
    t0::Float64,
    t1::Float64,
) where {T}
    indexT0 = searchsortedlast(surface.expiries, t0)
    indexT1 = searchsortedlast(surface.expiries, t1)

    t = t0
    #out = map(y -> varianceSlopeCurvature(surface, y, t, indexT0), y)
    #w = getindex.(out, 1)
    #dw = getindex.(out, 2)
    #d2w = getindex.(out, 3)
    w, dw, d2w = varianceSlopeCurvature(surface, y, t, indexT0)
    denom = @. ifelse(
        w == 0,
        1,
        1 - y / w * dw + 0.25 * (-0.25 * t^2 - t / w + y^2 / w^2) * dw^2 + 0.5 * d2w * t,
    )
    #w1 = map(y -> varianceByLogmoneyness(surface, y, t1, indexT1), y)
    w1 = varianceByLogmoneyness(surface, y, t1, indexT1)
    #may be worth using mid point t= (t0+t1)/2 and w(t0) below.
    num = @. (w1 * t1 - w * t0) / (t1 - t0)
    lv = map(
        (num, denom) -> ifelse(num < 0, 0, ifelse(denom < 0, 10.0, num / denom)),
        num,
        denom,
    )

end


text:
module HeatTransfer2D

using Modia

# Desired:
#   using ModiaMath: plot
#
# In order that ModiaMath need not to be defined in the user environment, it is included via Modia:
using Modia.ModiaMath: plot

# Definition of model
const N      = 30 # 100 # Number of nodes in x-direction (same number of nodes in y-direction
const L      = 0.2   # Length in x-, and y-direction
const T0     = 290.0 # Initial temperature
const TNx    = 310.0 # Temperature along the y-line outside x-node N
const TNy    = 330.0 # Temperature along the x-line outside y-node N
const cp     = 910.0 # Material Heat Capacity
const lambda = 237.0 # Material thermal conductivity
const rho    = 271.0 # Material Density
const dL     = L / N # Element length in x- and y-direction
const c      = lambda/(cp*rho*dL*dL)
const N2     = N*N

geti(i,j)  = (j-1)*N + i

""" der(T) = heatTransfer2D(T)"""
function heatTransfer2D(T::Matrix{Float64})
   der_T::Matrix{Float64} = zeros(N,N)
   for i in 1:N
      for j in 1:N
         qx1 = if i > 1 ; T[i - 1, j] - T[i,j] else 0.0              end
         qx2 = if i < N ; T[i + 1, j] - T[i,j] else 2*(TNx - T[i,j]) end
         qy1 = if j > 1 ; T[i, j - 1] - T[i,j] else 0.0              end
         qy2 = if j < N ; T[i, j + 1] - T[i,j] else 2*(TNy - T[i,j]) end

         der_T[i,j] = c*(qx1 + qx2 + qy1 + qy2)
      end
   end
   return der_T
end

# Not used presently:
""" Structure of der( heatTransfer2D, T ) """
function jacobian_incidence(::typeof(heatTransfer2D), args...)
   I::Vector{Int} = fill(1,5*N2)
   J::Vector{Int} = fill(1,5*N2)

   r = 0
   for i in 1:N
      for j in 1:N
         i0 = geti(i,j)
         i1 = geti(max(i-1,1),j)
         i2 = geti(min(i+1,N),j)
         i3 = geti(i,max(j-1,1))
         i4 = geti(i,min(j+1,N))

         r += 1
         I[r] = i0
         J[r] = i0

         r+= 1
         I[r] = i0
         J[r] = i1

         r+= 1
         I[r] = i0
         J[r] = i2

         r+= 1
         I[r] = i0
         J[r] = i3

         r+= 1
         I[r] = i0
         J[r] = i4
      end
   end

   return sparse(I,J,1)
end

@model HeatTransfer begin
   T = Float(start=fill(T0,N,N), fixed=true) # Temperature of the nodes
@equations begin
   der(T) = heatTransfer2D(T)
  end
end

result = simulate(HeatTransfer, 30)
#plot(result, "T", figure=22)

res = Dict{AbstractString,Any}()
res["time"] = result["time"]
T = result["T"]
T = reshape(T, size(T,1), N, N)
v1 = "T[$(div(N,2)),1]"
v2 = "T[1, $(div(N,2))]"
v3 = "T[$N,1]"
v4 = "T[1, $N]"
res[v1] = T[:, div(N,2),1]
res[v2] = T[:, 1, div(N,2)]
res[v3] = T[:, N, 1]
res[v4] = T[:, 1, N]

plot(res, (v1, v2, v3, v4), figure=20, heading="HeatTransfer2D: N=$N, temperature [K]")

end


text:
function check_reshop(m)
    error("Solver $m is not of type ReSHOP")
end

# trickery here
# if we store emp directly, print blows up ...
function setemp_solver!(m::ReSHOP.Optimizer, emp)
#    ReSHOP.reshop_setemp!(m, function () return _solveEMP(emp) end)
    return true
end

function setemp!(m, emp)
    error("Model $(typeof(m)) not supported: only JuMP and Convex are")
end

function setemp!(m::JuMP.Model, emp)
    #setemp_solver!(m.moi_backend.optimizer.model, emp)
end

function _getrealmoibackend(o::ReSHOP.Optimizer)
    return o
end

function _getrealmoibackend(o::MOI.Bridges.AbstractBridgeOptimizer)
    return o.model
end

function _getrealmoibackend(o::MOI.Utilities.CachingOptimizer)
    return _getrealmoibackend(o.optimizer)
end

function _getrealmoibackend(o)
    error("Unknown Optimizer of type $(typeof(o))")
end

function getReSHOPModel(emp::EMPmaster)
    jmodel = emp.model_ds

    return _getrealmoibackend(jmodel.moi_backend)
end

function _setvalues(m, x::Vector) error("_setvalues: unsupported type $(typeof(m))") end

function _setvalues(m::JuMP.Model, x::Vector)
    @assert length(m.colVal) == length(x)
    m.colVal[:] = x[:]
end

function setvarnames(ctx, m)
end

function setvarnames(ctx, m::JuMP.Model)
    if isnothing(m.internalModel.inner.d)
        ReSHOP.ctx_setvarnames(ctx, m.colNames)
    end
end


text:
#=
    As yet unused:
=#

isreduced(x::FastRational{T,IsReduced}) where {T} = true
isreduced(x::FastRational{T,MayReduce}) where {T} = false

mayreduce(x::FastRational{T,IsReduced}) where {T} = false
mayreduce(x::FastRational{T,MayReduce}) where {T} = true

# denominator is strictly positive
# numerator   is negative, zero, or positve

iszero(x::FastRational{T,IsReduced}) where {T} = iszero(x.num)
isone(x::FastRational{T,IsReduced})  where {T} = x.num === x.den

isinteger(x::FastRational{T,IsReduced}) where {T} = x.den == 1
isinteger(x::FastRational{T,MayReduce}) where {T} = canonical(x.num,x.den)[2] == 1


text:
click_dot!(game, dot, gui::Nothing) = click_dot!(game, dot)

function click_dot!(game, dot)
    !can_select(game, dot) ? (return false) : nothing
    if dot.selected
        dot.selected = false
        remove_dot!(game, dot)
    else
        dot.selected = true
        add_dot!(game, dot)
    end
    return true
end

function can_select(game, dot)
    game_over!(game) ? (return false) : nothing
    dots = game.selected_dots
    # first dot
    if isempty(dots)
        return true
    end
    # unclick last last clicked dot
    if dot == dots[end]
        return true
    end
    # cannot select dot again
    if dot ∈ dots
        return false
    end
    # colors do not match
    if dots[end].color != dot.color
        return false
    end
    # not adjecent
    if !is_adjecent(dots[end], dot)
        return false
    end
    return true
end

function is_adjecent(dot1, dot2)
    if (dot1.row + 1 == dot2.row) && (dot1.col == dot2.col)
        return true
    end
    if (dot1.row - 1 == dot2.row) && (dot1.col == dot2.col)
        return true
    end
    if (dot1.row == dot2.row) && (dot1.col == dot2.col + 1)
        return true
    end
    if (dot1.row  == dot2.row) && (dot1.col == dot2.col - 1)
        return true
    end
    return false
end

function max_coords(dots)
    min_row = minimum(x->x.row, dots)
    max_row = maximum(x->x.row, dots)
    min_col = minimum(x->x.col, dots)
    max_col = maximum(x->x.col, dots)
    return min_row, max_row, min_col, max_col
end

is_rectangular(game) = is_rectangular(game.selected_dots, game.dots)

function is_rectangular(selected_dots, dots)
    min_row, max_row, min_col, max_col = max_coords(selected_dots)
    if ((max_row - min_row) == 0) || ((max_col - min_col) == 0)
        return false
    end
    return rectangle_sizes(dots, min_row, max_row, min_col, max_col)
end

function rectangle_sizes(dots, min_row, max_row, min_col, max_col)
    max_height = max_row - min_row + 1
    max_width = max_col - min_col + 1
    for width in 2:max_width
        for height in 2:max_height
            if rectangle_offset(dots, width, height, min_row, max_height, min_col, max_width)
                return true
            end
        end
    end
    return false
end

function rectangle_offset(dots, width, height, min_row, max_height, min_col, max_width)
    for c in 0:(max_width - width)
        for r in 0:(max_height - height)
            if trace_rectangle(
                    dots,
                    min_row + r,
                    height + min_row + r - 1,
                    min_col + c,
                    width + min_col + c - 1
                )
                return true
            end
        end
    end
    return false
end

function trace_rectangle(dots, min_row, max_row, min_col, max_col)
    for r in [min_row,max_row]
        for c in min_col:max_col
            !dots[r,c].selected ? (return false) : nothing
        end
    end
    for c in [min_col,max_col]
        for r in min_row:max_row
            !dots[r,c].selected ? (return false) : nothing
        end
    end
    return true
end

function rand_color!(x::Dot)
    x.color = rand(["red","blue","yellow","green","orange"])
    return nothing
end

populate(n_rows, n_cols) = [Dot(;row, col) for row in 1:n_rows, col in 1:n_cols]

game_over!(game) = game.round == 0

set_unselected!(game::Game) = set_unselected!(game.selected_dots)

function set_unselected!(dots)
    map(x->x.selected=false, dots)
end

function update_score!(game)
    game.score += length(game.selected_dots)
end

function update_round!(game)
    game.round -= 1
end

remove_dot!(game::Game, dot) = remove_dot!(game.selected_dots, dot)

function remove_dot!(selected_dots, dot)
    filter!(x->x != dot, selected_dots)
end

add_dot!(game::Game, dot) = add_dot!(game.selected_dots, dot)

function add_dot!(selected_dots, dot)
    push!(selected_dots, dot)
end

select_all_color!(game) = select_all_color!(game.dots, game.selected_dots)

function select_all_color!(dots, selected_dots)
    color = selected_dots[1].color
    color_dots = filter(x->x.color == color, dots)
    map(x->x.selected = true, color_dots)
    return color_dots
end

function clear_selected!(game)
    empty!(game.selected_dots)
end

sort_by_row!(game::Game) = sort_by_row!(game.selected_dots)

function sort_by_row!(dots)
    rows = map(x->x.row, dots)
    idx = sortperm(rows)
    dots .= dots[idx]
end

not_connected(game) = length(game.selected_dots) < 2

function shift_colors!(game)
    for dot in game.selected_dots
        shift_color!(dot, game)
    end
end

function shift_color!(dot, game)
    dots = game.dots
    row = dot.row
    col = dot.col
    while (row > 1)
        dots[row,col].color = dots[row-1,col].color
        row -= 1
    end
    rand_color!(dots[row,col])
end

click_submit!(game, gui::Nothing) = click_submit!(game)

function click_submit!(game)
    game_over!(game) ? (return false) : nothing
    not_connected(game) ? (return false) : nothing
    update_round!(game)
    if is_rectangular(game)
        game.selected_dots = select_all_color!(game)
    end
    update_score!(game)
    sort_by_row!(game)
    shift_colors!(game)
    set_unselected!(game)
    clear_selected!(game)
    return true
end


text:
#=====
Serialization
	Tools
=====#
module data
using DataFrames
struct ImageDataFrame
    filenames::Vector{String}
    labels::Vector
	resize::Union{Nothing,Tuple}

	function ImageDataFrame(filenames, labels; resize=nothing)
		@assert length(filenames) == length(labels)
		try
			@eval import Images
			@eval import ImageMagick
		catch
			@warn "Package Images or ImageMagick not installed"
		end
		new(filenames, labels, resize)
	end

end

Base.length(ds::ImageDataFrame) = length(ds.filenames)

function Base.getindex(ds::ImageDataFrame, idx)
	filename = ds.filenames[idx]
	# img = Images.load(filename) has multi-thread issue
    img = ImageMagick.load(filename)
	img = Images.channelview(img)
	img = permutedims(img, [2,3,1])
	img = convert(Array{Float32}, img)
	if ds.resize !== nothing
		img = Images.imresize(img, ds.resize...)
	end
	return (img, ds.labels[idx])
end

struct DFDataFrame
	df
    X::Vector{Symbol}
	Y::Vector{Symbol}

	DFDataFrame(df,X,Y) = new(df, makeArray(X), makeArray(Y))
end

Base.length(ds::DFDataFrame) = size(df,1)

function Base.getindex(ds::DFDataFrame, idx)
	df = ds.df
	(Vector(df[idx, ds.X]), Vector(df[idx, ds.Y]))
end
#==
TRANSFORMERS
==#
abstract type Transformer end

Base.length(t::Transformer) = length(t.ds)
mutable struct NoisingTransfomer
    ds::Union{DataFrame, Nothing}
	noiselevel
	axis

    NoisingTransfomer(noiselevel=0.01; axis=[:X]) = new(nothing,noiselevel,axis)
end

genNoise(arr::AbstractArray, level::Number) = randn(eltype(arr), size(arr)) .* level

function Base.getindex(t::NoisingTransfomer, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X += genNoise(X, t.noiselevel)
	end

	if :Y in t.axis
		Y += genNoise(Y, t.noiselevel)
	end

	return (X,Y)
end
mutable struct Normalizer <: Transformer
    ds::Union{DataFrame, Nothing}
	dims::Array{Int}
	means
	stds
	axis

    Normalizer(means=0, stds=1, dims=[1,2]; axis=[:X]) =
		new(nothing,dims, means, stds, axis)
end

function normalize(x, means, stds)
	(x .- means) ./ stds
end

function Base.getindex(t::Normalizer, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X += normalize(X, t.means, t.stds)
	end

	if :Y in t.axis
		Y += normalize(Y, t.means, t.stds)
	end

	return (X,Y)
end
mutable struct ImageCrop <: Transformer
    ds::Union{DataFrame, Nothing}
	shapeX
	shapeY

    ImageCrop(shapeX , shapeY=nothing) = new(nothing,shapeX, shapeY)
end

function imageCrop(arr::AbstractArray,shape)
	maxX = size(arr,1) - shape[1]
	maxY = size(arr,2) - shape[2]
	offX = rand(1:maxX)
	offY = rand(1:maxY)
	arr[offX:offX+shape[1]-1,offY:offY+shape[2]-1, :]
end

function Base.getindex(t::ImageCrop, idx)
	X,Y = t.ds[idx]
	if t.shapeX !== nothing
		X = imageCrop(X, t.shapeX)
	end

	if t.shapeY !== nothing
		Y = imageCrop(Y, t.shapeY)
	end

	return (X,Y)
end

function onehot(x, labels, dtype::Type)
	result = zeros(dtype, length(labels))
    result[findfirst(x .== labels)] = 1
    result
end

function onehot(X::AbstractArray, labels, dtype::Type)
	result = zeros(dtype, length(labels))
	for x in X
		result[findfirst(x .== labels)] = 1
	end
	result
end

mutable struct OneHotEncoder <: Transformer
    ds::Union{DataFrame, Nothing}
	labels
	axis
	dtype::Type

    OneHotEncoder(labels; axis=[:Y], dtype=getContext().dtype) =
		new(nothing, labels, axis, dtype)
end

function Base.getindex(t::OneHotEncoder, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X = onehot(X, t.labels, t.dtype)
	end

	if :Y in t.axis
		Y = onehot(Y, t.labels, t.dtype)
	end

	return (X,Y)
end
mutable struct Subset <: Transformer
    ds::Union{DataFrame, Nothing}
	idxs::AbstractArray{Int}

    Subset(idxs::AbstractArray{Int}) = new(nothing, idxs)
end

Base.length(t::Subset) = length(t.idxs)
Base.getindex(t::Subset, idx) = t.ds[t.idxs[idx]]

mutable struct Split <: Transformer
    ds::Union{DataFrame, Nothing}
	valid_perc::Float64
	shuffle::Bool

    Split(valid_perc=0.2; shuffle=true) =
		new(nothing, valid_perc, shuffle)
end

function (t::Split)(ds::DataFrame)
    t.ds = ds
	maxl = length(ds)

	idxs = t.shuffle ? Random.shuffle(1:maxl) : 1:maxl
	cut = round(Int, maxl*(1.0-t.valid_perc))

    a, b = (Subset(idxs[1:cut]), Subset(idxs[cut+1:end]))
	return (a(ds), b(ds))
end

#=
helper functions
=#
function update_mb!(arr::AbstractArray, elem::AbstractArray, idx::Int)
	@assert size(arr)[1:end-1] == size(elem) "$(size(arr)) $(size(elem))"
	idxs = Base.OneTo.(size(elem))
	arr[idxs..., idx] = elem
end

function update_mb!(t::Tuple, elems::Tuple, idx::Int)
	@assert length(t) == length(elems)
	for (arr,elem) in zip(t,elems)
		update_mb!(arr, elem, idx)
	end
end
create_mb(arr::AbstractArray, batchsize) = similar(arr, size(arr)..., batchsize)
create_mb(t::Tuple, batchsize)= Tuple(collect(create_mb(elem, batchsize) for elem in t))
mutable struct MiniBatch <: Transformer
    ds::Union{DataFrame, Nothing}
    batchsize::Int
    shuffle::Bool

    MiniBatch(batchsize=8; shuffle=true) =
        new(nothing, batchsize, shuffle)
end

Base.length(dl::MiniBatch) = length(dl.ds) ÷ dl.batchsize

function Base.iterate(dl::MiniBatch, state=undef)
    maxl = length(dl.ds)
    bs = dl.batchsize

    if state == undef
        idxs = dl.shuffle ? Random.shuffle(1:maxl) : 1:maxl
        state = (idxs,1)
    end
    idxs, count = state

    if count > (maxl-bs) return nothing end

	l = Threads.SpinLock()
	minibatch = nothing

    Threads.@threads for i in 1:bs

		idx = i + count - 1
		sample = dl.ds[idx]
		@assert sample isa Tuple "DataFrames should return Tuples, not $(typeof(sample))"

		if minibatch === nothing
			Threads.lock(l)
			if minibatch === nothing
				minibatch = create_mb(sample, bs)
			end
			Threads.unlock(l)
		end

		update_mb!(minibatch, sample, i)
    end
	Threads.unlock(l)
    return ((minibatch), (idxs, count + bs))
end
#-----------------------------------------------------------
end


text:
module InheritanceTest
using ComputedFieldTypes
using Test

abstract type Bar end

@computed struct Parametric{T} <: Bar
    x::Base.promote_op(+, T, Float64)
end

@computed struct NonParametric <: Bar
    x::Base.promote_op(+, Int, Float64)
end

@computed struct NonParametricNoParent
    x::Base.promote_op(+, Int, Float64)
end

@testset "inheritance" begin
    @test isa(@inferred(Parametric{Int}(1)), Parametric{Int, Float64})
    @test isa(Parametric{Int}(1), Bar)
    @test isa(@inferred(Parametric{Float64}(1.0)), Parametric{Float64, Float64})

    @test isa(@inferred(NonParametric(1)), NonParametric{Float64})
    @test isa(NonParametric(1), Bar)

    @test isa(@inferred(NonParametricNoParent(1)), NonParametricNoParent{Float64})
    @test !isa(NonParametricNoParent(1), Bar)
end

@computed struct ParametricRef{T}
    x::Base.RefValue{Some{T}}
    y::getfield(Base,:RefValue){Some{T}}
    ParametricRef(x::Some{T}) where {T} = new{T}(Ref(x), Ref(x))
end

@computed struct ParametricWhere{T}
    a::Pair{S, Some{T}} where S
    b::Pair{S, Some{T}} where S<:typeassert(T,DataType)
    c::Pair{S, Some{T}} where S>:(T::DataType)
    d::Pair{S, Some{T}} where typeassert(T,DataType)<:S<:typeassert(T,DataType)
    ## this is currently illegal, since we cannot determine when to compute identity(S):
    # x::Pair{identity(S), Some{T}} where S
    function ParametricWhere(x::Some{T}) where {T}
        p = something(x) => x
        return new{T}(p, p, p, p)
    end
end

@testset "nested expressions" begin
    @test isa(@inferred(ParametricRef(Some(1.0))), ParametricRef{Float64})
    T = typeof(ParametricRef(Some(1)))
    @test T === fulltype(ParametricRef{Int})
    @test isa(ParametricRef(Some(1)).x[], Some{Int})
    @test isa(ParametricRef(Some(1)).y[], Some{Int})

    @test isa(@inferred(ParametricWhere(Some(1.0))), ParametricWhere{Float64})
    T = typeof(ParametricWhere(Some(1)))
    @test T === fulltype(ParametricWhere{Int})
    @test fieldtype(T, :a) === Pair{S, Some{Int}} where S
    @test fieldtype(T, :b) === Pair{S, Some{Int}} where S<:Int
    @test fieldtype(T, :c) === Pair{S, Some{Int}} where S>:Int
    @test fieldtype(T, :d) === Pair{S, Some{Int}} where Int<:S<:Int
end

end

max_stars_repo_name: vtjnash/ComputedFieldTypes

max_issues_repo_name: vtjnash/ComputedFieldTypes

max_forks_repo_name: vtjnash/ComputedFieldTypes


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

"""
IN PROGRESS:
The mutable struct representing an FMU in the FMI 1.0 Standard.
Also contains the paths to the FMU and ZIP folder as well als all the FMI 1.0 function pointers
"""
mutable struct FMU1

    # Constructor
    fmu1() = new()
end


text:
using YaoBlocks
using LuxurySparse
using LinearAlgebra
using ..SymEngine
using ..SymEngine: BasicType, BasicOp, BasicTrigFunction

op_types = [:Mul, :Add, :Pow]
const BiVarOp = Union{[SymEngine.BasicType{Val{i}} for i in op_types]...}

export @vars

simag = SymFunction("Im")
sreal = SymFunction("Re")
sabs = SymFunction("abs")

Base.promote_rule(::Type{Bool}, ::Type{Basic}) = Basic
Base.conj(x::Basic) = Basic(conj(SymEngine.BasicType(x)))
Base.conj(x::BasicType) = real(x) - im * imag(x)
Base.conj(x::BiVarOp) = juliafunc(x)(conj.(get_args(x.x))...)
Base.conj(x::BasicTrigFunction) = juliafunc(x)(conj.(get_args(x.x)...)...)
# WARNING: symbols and constants are assumed real!
Base.imag(x::BasicType{Val{:Constant}}) = Basic(0)
Base.imag(x::BasicType{Val{:Symbol}}) = Basic(0)
Base.abs(x::Basic) = Basic(abs(SymEngine.BasicType(x)))
Base.abs(x::BasicType{Val{:Constant}}) = x
Base.abs(x::BasicType{Val{:Symbol}}) = x
function Base.imag(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(imag, +, args)
end

function Base.real(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(real, +, args)
end

function Base.abs(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(abs, +, args)
end

function Base.imag(x::BasicType{Val{:Mul}})
    args = (get_args(x.x)...,)
    get_mul_imag(args)
end

function Base.real(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    if imag(a) == 0 && imag(b) == 0
        return x.x
    else
        if imag(a) == 0
            return a^real(b) * cos(log(a) * imag(b))
        else
            return sreal(x.x)
        end
    end
end

function Base.imag(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    if imag(a) == 0 && imag(b) == 0
        return Basic(0)
    else
        if imag(a) == 0
            return a^real(b) * sin(log(a) * imag(b))
        else
            return simag(x.x)
        end
    end
end

function Base.abs(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    abs(a)^real(b)
end

function Base.real(x::BasicType{Val{:Mul}})
    args = (get_args(x.x)...,)
    get_mul_real(args)
end

function get_mul_imag(args::NTuple{N,Any}) where {N}
    imag(args[1]) * get_mul_real(args[2:end]) + real(args[1]) * get_mul_imag(args[2:end])
end
get_mul_imag(args::Tuple{Basic}) = imag(args[1])

function get_mul_real(args::NTuple{N,Any}) where {N}
    real(args[1]) * get_mul_real(args[2:end]) - imag(args[1]) * get_mul_imag(args[2:end])
end
get_mul_real(args::Tuple{Basic}) = real(args[1])

function Base.real(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return x.x
    else
        return sreal(x.x)
    end
end

function Base.abs(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return x.x
    else
        return sabs(x.x)
    end
end

function Base.imag(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return Basic(0)
    else
        return simag(x.x)
    end
end

@generated function juliafunc(x::BasicType{Val{T}}) where {T}
    SymEngine.map_fn(T, SymEngine.fn_map)
end

const SymReal = Union{Basic,SymEngine.BasicRealNumber}
YaoBlocks.RotationGate(block::GT, theta::T) where {N,T<:SymReal,GT<:AbstractBlock{N}} =
    RotationGate{N,T,GT}(block, theta)

YaoBlocks.phase(θ::SymReal) = PhaseGate(θ)
YaoBlocks.shift(θ::SymReal) = ShiftGate(θ)

YaoBlocks.mat(::Type{Basic}, ::HGate) = 1 / sqrt(Basic(2)) * Basic[1 1; 1 -1]
YaoBlocks.mat(::Type{Basic}, ::XGate) = Basic[0 1; 1 0]
YaoBlocks.mat(::Type{Basic}, ::YGate) = Basic[0 -1im; 1im 0]
YaoBlocks.mat(::Type{Basic}, ::ZGate) = Basic[1 0; 0 -1]
YaoBlocks.mat(::Type{Basic}, gate::ShiftGate) = Diagonal([1, exp(im * gate.theta)])
YaoBlocks.mat(::Type{Basic}, gate::PhaseGate) = exp(im * gate.theta) * IMatrix{2}()
function YaoBlocks.mat(::Type{Basic}, R::RotationGate{N}) where {N}
    I = IMatrix{1 << N}()
    return I * cos(R.theta / 2) - im * sin(R.theta / 2) * mat(Basic, R.block)
end
for GT in [:XGate, :YGate, :ZGate]
    @eval YaoBlocks.mat(::Type{Basic}, R::RotationGate{1,T,<:$GT}) where {T} =
        invoke(mat, Tuple{Type{Basic},RotationGate}, Basic, R)
end

for T in [:(RotationGate{N,<:SymReal}), :(PhaseGate{<:SymReal}), :(ShiftGate{<:SymReal})]
    @eval YaoBlocks.mat(gate::$T) = mat(Basic, gate)
end

YaoBlocks.PSwap{N}(locs::Tuple{Int,Int}, θ::SymReal) where {N} =
    YaoBlocks.PutBlock{N}(rot(ConstGate.SWAPGate(), θ), locs)

YaoBlocks.pswap(n::Int, i::Int, j::Int, α::SymReal) = PSwap{n}((i, j), α)
YaoBlocks.pswap(i::Int, j::Int, α::SymReal) = n -> pswap(n, i, j, α)

export subs, chiparams
chiparams(blk::RotationGate, param) = rot(blk.block, param)
chiparams(blk::ShiftGate, param) = shift(param)
chiparams(blk::PhaseGate, param) = phase(param)
chiparams(blk::TimeEvolution, param) = time_evolve(blk.H, param, tol = blk.tol)
chiparams(blk::AbstractBlock, params...) =
    niparams(blk) == length(params) == 0 ? blk :
    throw(NotImplementedError(:chiparams, (blk, params...)))

SymEngine.subs(c::AbstractBlock, args...; kwargs...) = subs(Basic, c, args...; kwargs...)
function SymEngine.subs(::Type{T}, c::AbstractBlock, args...; kwargs...) where {T}
    c = chiparams(c, map(x -> T(subs(x, args...; kwargs...)), getiparams(c))...)
    chsubblocks(c, [subs(T, blk, args..., kwargs...) for blk in subblocks(c)])
end

# dumpload
YaoBlocks.tokenize_param(param::Basic) = QuoteNode(Symbol(param))
YaoBlocks.parse_param(x::QuoteNode) = :(Basic($x))


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

using Markdown
using InteractiveUtils

# ╔═╡ 794a5356-6e05-44b4-ad9d-375a1fa3613b
begin
	using PlutoUI
	using DataFrames

	using Test
end

# ╔═╡ b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
md"""
## Introduction DataFrames / 9 - Reshaping DataFrames

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/09_reshaping.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ d6993a5b-8158-4ce3-accb-8e8c514229b3
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ de97ded0-3952-4188-9a43-a8f383e1ae65
md"""
### Wide to long
"""

# ╔═╡ 3ecd5008-03f0-4727-abba-345b8ac08e1d
x = DataFrame(id=[1, 2, 3, 4], id2=[1, 1, 2, 2], M1=[11, 12, 13, 14], M2=[111, 112, 113, 114])

# ╔═╡ b294ecbe-2f49-4e5d-bda9-78096463da07
stack(x, [:M1, :M2], :id)  ## first pass measure variables and then id-variable

# ╔═╡ 53eb03ad-2c1d-44e6-9f4d-e56c121d0845
md"""
add `view=true` keyword argument to make a view; in that case columns of the resulting dataframe share memory with columns of the source data frame, so the operation is *potentially unsafe*
"""

# ╔═╡ 59d52da1-01b9-4df2-b18a-db87f834eaac
## optionally you can rename columns
stack(x, ["M1", "M2"], "id", variable_name="key", value_name="observed")

# ╔═╡ 8b34aaef-4a4e-4751-ac4f-4765910244aa
md"""
if second argument is omitted in stack, all other columns are assumed to be the id-variables
"""

# ╔═╡ 196f3066-f8a8-40c4-8db8-719108b86e38
stack(x, Not([:id, :id2]))

# ╔═╡ 8dae04b5-bb93-4f1a-9272-84ac47935394
stack(x, Not([1, 2]))  ## we can use index instead of symbol

# ╔═╡ f8eeacd5-5085-4652-9aef-d31a128fe682
nx = DataFrame(id=[1, 1, 1], id2=['a', 'b', 'c'], a1=rand(3), a2=rand(3))

# ╔═╡ f2ad2b0f-ff3f-49b0-8c4a-b28359ebabe7
md"""
if stack is not passed any measure variables by default numeric variables are selected as measures
"""

# ╔═╡ 25116592-668c-4ce9-b3d6-c04faab3ea21
stack(nx)

# ╔═╡ 0b94505f-3e6d-427f-939f-72c2bbf0aab2
md"""
Here all columns are treated as measures:
"""

# ╔═╡ db1ec4bd-60c5-4dac-a8a7-63731e4199d1
stack(DataFrame(rand(3, 2), :auto))

# ╔═╡ 1a78b75c-1c77-4c81-81d1-0769f7e2542f
begin
	df = DataFrame(rand(3, 2), :auto)
	df.key = [1, 1, 1]
	mdf = stack(df)                       ## duplicates in key are silently accepted
end

# ╔═╡ ca9927e5-3ff2-4c73-8206-4d7495668312
md"""
### Long to wide
"""

# ╔═╡ 806fff39-f475-4c30-8964-a7858c91bc81
x₁ = DataFrame(id=[1, 1, 1], id2=['a', 'b', 'c'], a1=rand(3), a2=rand(3))

# ╔═╡ 24e09155-38ec-4afb-82ef-17ab32d149fc
y₁ = stack(x₁)

# ╔═╡ c5ae24cf-c09d-45c7-bf45-3161250af616
unstack(y₁, :id2, :variable, :value)  ## standard unstack with a specified key

# ╔═╡ a16e6b3d-2a20-4f04-b34a-9428b8edc853
unstack(y₁, :variable, :value)        ## all other columns are treated as keys

# ╔═╡ 16886bba-11e7-4d71-8df6-6f958c4ee1a4
unstack(y₁)    ## all columns other than named :variable and :value are treated as keys

# ╔═╡ c986ec52-c3b8-4f8b-af87-df04c70c46a4
unstack(y₁, renamecols=n->string("unstacked_", n))  ## we can rename the unstacked columns

# ╔═╡ 3dad50c9-e60b-468c-be15-eab978c80128
df₁ = stack(DataFrame(rand(3, 2), :auto))

# ╔═╡ a9ebc088-20ab-44f9-88a9-aac9964796d7
@test_throws ArgumentError unstack(df₁, :variable, :value)
## unable to unstack when no key column is present

# ╔═╡ c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e
html"""

"""

# ╔═╡ Cell order:
# ╟─b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
# ╠═794a5356-6e05-44b4-ad9d-375a1fa3613b
# ╟─d6993a5b-8158-4ce3-accb-8e8c514229b3
# ╟─de97ded0-3952-4188-9a43-a8f383e1ae65
# ╠═3ecd5008-03f0-4727-abba-345b8ac08e1d
# ╠═b294ecbe-2f49-4e5d-bda9-78096463da07
# ╟─53eb03ad-2c1d-44e6-9f4d-e56c121d0845
# ╠═59d52da1-01b9-4df2-b18a-db87f834eaac
# ╟─8b34aaef-4a4e-4751-ac4f-4765910244aa
# ╠═196f3066-f8a8-40c4-8db8-719108b86e38
# ╠═8dae04b5-bb93-4f1a-9272-84ac47935394
# ╠═f8eeacd5-5085-4652-9aef-d31a128fe682
# ╟─f2ad2b0f-ff3f-49b0-8c4a-b28359ebabe7
# ╠═25116592-668c-4ce9-b3d6-c04faab3ea21
# ╟─0b94505f-3e6d-427f-939f-72c2bbf0aab2
# ╠═db1ec4bd-60c5-4dac-a8a7-63731e4199d1
# ╠═1a78b75c-1c77-4c81-81d1-0769f7e2542f
# ╟─ca9927e5-3ff2-4c73-8206-4d7495668312
# ╠═806fff39-f475-4c30-8964-a7858c91bc81
# ╠═24e09155-38ec-4afb-82ef-17ab32d149fc
# ╠═c5ae24cf-c09d-45c7-bf45-3161250af616
# ╠═a16e6b3d-2a20-4f04-b34a-9428b8edc853
# ╠═16886bba-11e7-4d71-8df6-6f958c4ee1a4
# ╠═c986ec52-c3b8-4f8b-af87-df04c70c46a4
# ╠═3dad50c9-e60b-468c-be15-eab978c80128
# ╠═a9ebc088-20ab-44f9-88a9-aac9964796d7
# ╠═c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
import Base.var, Base.std

# mcvar and mcse stand for Monte Carlo variance and Monte Carlo error respectively
export mcvar, var, mcse, std

# Variance of MCMCChain
mcvar_iid(c::MCMCChain, pars::Ranges=1:ncol(c.samples)) =
  Float64[var(c.samples[:, pars[i]]) for i = 1:pars.len]/nrow(c.samples)

mcvar_iid(c::MCMCChain, par::Real) = mcvar_iid(c, par:par)

# Standard deviation of MCMCChain
msce_iid(c::MCMCChain, pars::Ranges=1:ncol(c.samples)) = sqrt(mcvar_iid(c, pars))

msce_iid(c::MCMCChain, par::Real) = msce_iid(c, par:par)

# Function for estimating Monte Carlo error using batch means, see for instance
# Flegal J.M, Jones, G.L. Batch Means and Spectral Variance Estimators in Markov chain Monte Carlo. Annals of
# Statistics, 2010, 38 (2), pp 1034-1070
function mcvar_bm(x::AbstractVector; batchlen::Int=100)
  nbatches = div(length(x), batchlen)
  @assert nbatches > 1 "Choose batch size such that the number of batches is greather than one"
  nbsamples = nbatches*batchlen
  batchmeans = Float64[mean(x[((j-1)*batchlen+1):(j*batchlen)]) for j = 1:nbatches]
  return batchlen*var(batchmeans)/nbsamples
end

mcvar_bm(x::AbstractMatrix, pars::Ranges=1:size(x, 2); batchlen::Int=100) =
  Float64[mcvar_bm(x[:, pars[i]]; batchlen=batchlen) for i = 1:pars.len]

mcvar_bm(x::AbstractMatrix, par::Real; batchlen::Int=100) = mcvar_bm(x, par:par; batchlen=batchlen)

mcvar_bm(c::MCMCChain, pars::Ranges=1:ncol(c.samples); batchlen::Int=100) =
  Float64[mcvar_bm(c.samples[:, pars[i]]; batchlen=batchlen) for i = 1:pars.len]

mcvar_bm(c::MCMCChain, par::Real; batchlen::Int=100) = mcvar_bm(c, par:par; batchlen=batchlen)

# Monte Carlo standard error using batch means
mcse_bm(c::MCMCChain, pars::Ranges=1:ncol(c.samples); batchlen::Int=100) = sqrt(mcvar_bm(c, pars; batchlen=batchlen))

mcse_bm(c::MCMCChain, par::Real; batchlen::Int=100) = mcse_bm(c, par:par; batchlen=batchlen)

# Function for estimating the variance of a single MCMC chain using the initial monotone sequence estimator (IMSE), see
# Geyer C.J. Practical Markov Chain Monte Carlo. Statistical Science, 1992, 7 (4), pp 473-483
function mcvar_imse(x::AbstractVector; maxlag::Int=length(x)-1)
  k = convert(Int, floor((maxlag-1)/2))
  m = k+1

  # Preallocate memory
  g = Array(Float64, k+1)

  # Calculate empirical autocovariance
  acv = acf(x, 0:maxlag, correlation=false)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Create the monotone sequence of g
  if m > 1
    for j = 2:m
      if g[j] > g[j-1]
        g[j] = g[j-1]
      end
    end
  end

  # Calculate the initial monotone sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(x)
end

mcvar_imse(x::AbstractMatrix, pars::Ranges=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcvar_imse(x[:, pars[i]]; maxlag=maxlag) for i = 1:pars.len]

mcvar_imse(x::AbstractMatrix, par::Real; maxlag::Int=size(x, 1)-1) = mcvar_imse(x, par:par; maxlag=maxlag)

mcvar_imse(c::MCMCChain, pars::Ranges=1:ncol(c.samples); maxlag::Int=nrow(c.samples)-1) =
  Float64[mcvar_imse(c.samples[:, pars[i]]; maxlag=maxlag) for i = 1:pars.len]

mcvar_imse(c::MCMCChain, par::Real; maxlag::Int=nrow(c.samples)-1) = mcvar_imse(c, par:par; maxlag=maxlag)

# Standard deviation using Geyer's IMSE
mcse_imse(c::MCMCChain, pars::Ranges=1:ncol(c.samples); maxlag::Int=nrow(c.samples)-1) =
  sqrt(var_imse(c, pars; maxlag=maxlag))

mcse_imse(c::MCMCChain, par::Real; maxlag::Int=nrow(c.samples)-1) = mcse_imse(c, par:par; maxlag=maxlag)

# Function for estimating the variance of a single MCMC chain using the initial positive sequence estimator (IPSE), see
# Geyer C.J. Practical Markov Chain Monte Carlo. Statistical Science, 1992, 7 (4), pp 473-483
function mcvar_ipse(x::AbstractVector; maxlag::Int=length(x)-1)
  k = convert(Int, floor((maxlag-1)/2))

  # Preallocate memory
  g = Array(Float64, k+1)
  m = k+1

  # Calculate empirical autocovariance
  acv = acf(x, 0:maxlag, correlation=false)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Calculate the initial monotone sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(x)
end

mcvar_ipse(x::AbstractMatrix, pars::Ranges=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcvar_ipse(x[:, pars[i]]; maxlag=maxlag) for i = 1:pars.len]

mcvar_ipse(x::AbstractMatrix, par::Real; maxlag::Int=size(x, 1)-1) = mcvar_ipse(x, par:par; maxlag=maxlag)

mcvar_ipse(c::MCMCChain, pars::Ranges=1:ncol(c.samples); maxlag::Int=nrow(c.samples)-1) =
  Float64[mcvar_ipse(c.samples[:, pars[i]]; maxlag=maxlag) for i = 1:pars.len]

mcvar_ipse(c::MCMCChain, par::Real; maxlag::Int=nrow(c.samples)-1) = mcvar_ipse(c, par:par; maxlag=maxlag)

# Standard deviation using Geyer's IMSE
mcse_ipse(c::MCMCChain, pars::Ranges=1:ncol(c.samples); maxlag::Int=nrow(c.samples)-1) =
  sqrt(mcvar_ipse(c, pars; maxlag=maxlag))

mcse_ipse(c::MCMCChain, par::Real; maxlag::Int=nrow(c.samples)-1) = mcse_ipse(c, par:par; maxlag=maxlag)

# Wrapper function for computing MCMC variance using various approaches
vtypes = (:bm, :iid, :imse, :ipse)

function var(c::MCMCChain, pars::Ranges=1:ncol(c.samples); vtype::Symbol=:imse, args...)
  @assert in(vtype, vtypes) "Unknown variance type $vtype"

  if vtype == :bm
    return mcvar_bm(c, pars; args...)
  elseif vtype == :iid
    return mcvar_iid(c, pars)
  elseif vtype == :imse
    return mcvar_imse(c, pars; args...)
  elseif vtype == :ipse
    return mcvar_ipse(c, pars; args...)
  end
end

var(c::MCMCChain, par::Real; vtype::Symbol=:imse, args...) = var(c, par:par; vtype=vtype, args...)

# Wrapper function for computing Monte Carlo (standard) error using various approaches
function std(c::MCMCChain, pars::Ranges=1:ncol(c.samples); vtype::Symbol=:imse, args...)
  @assert in(vtype, vtypes) "Unknown standard error type $vtype"

  if vtype == :bm
    return mcse_bm(c, pars; args...)
  elseif vtype == :iid
    return mcse_iid(c, pars)
  elseif vtype == :imse
    return mcse_imse(c, pars; args...)
  elseif vtype == :ipse
    return mcse_ipse(c, pars; args...)
  end
end