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

using Markdown
using InteractiveUtils

# ╔═╡ 6ca83e83-423a-4621-ae14-dc8da5913394
using Plots, PlutoTest

# ╔═╡ df80bc4c-770a-4acd-a4f3-337c6b720c71
plotly()

# ╔═╡ 20fb22d0-b0e9-11eb-1561-8b6505f5cb45
md"""
# [Day 18: Like a GIF For Your Yard](https://adventofcode.com/2015/day/18)

After the [million lights incident](https://adventofcode.com/2015/day/6), the fire code has gotten stricter: now, at most ten thousand lights are allowed. You arrange them in a 100x100 grid.

Never one to let you down, Santa again mails you instructions on the ideal lighting configuration. With so few lights, he says, you'll have to resort to animation.

Start by setting your lights to the included initial configuration (your puzzle input). A `#` means "on", and a `.` means "off".

Then, animate your grid in steps, where each step decides the next configuration based on the current one. Each light's next state (either on or off) depends on its current state and the current states of the eight lights adjacent to it (including diagonals). Lights on the edge of the grid might have fewer than eight neighbors; the missing ones always count as "off".

For example, in a simplified 6x6 grid, the light marked `A` has the neighbors numbered `1` through `8`, and the light marked `B`, which is on an edge, only has the neighbors marked `1` through `5`:

```
1B5...
234...
......
..123.
..8A4.
..765.

```

The state a light should have next is based on its current state (on or off) plus the number of neighbors that are on:

-   A light which is on stays on when `2` or `3` neighbors are on, and turns off otherwise.
-   A light which is off turns on if exactly `3` neighbors are on, and stays off otherwise.

All of the lights update simultaneously; they all consider the same current state before moving to the next.

Here's a few steps from an example configuration of another 6x6 grid:

```
Initial state:
.#.#.#
...##.
#....#
..#...
#.#..#
####..

After 1 step:
..##..
..##.#
...##.
......
#.....
#.##..

After 2 steps:
..###.
......
..###.
......
.#....
.#....

After 3 steps:
...#..
......
...#..
..##..
......
......

After 4 steps:
......
......
..##..
..##..
......
......

```

After `4` steps, this example has four lights on.

In your grid of 100x100 lights, given your initial configuration, **how many lights are on after 100 steps?**
"""

# ╔═╡ c73f4666-7302-4743-bb13-07ca4aea6440
puzzle_input = "#...##......#......##.##..#...##......##.#.#.###.#.#..#..#......####..#......###.#.#....#..##..###..
####..#.#...#....#.#####.##.##.#..#.......#....#.##...###.###..#.#.#........#..#.#.##...##..#.####.#
...#..##...#.#.###.#.###..#.##.####.###...#...........#.###..##.#.##.#.###...#.#..###....#.###.#..#.
.#...##...####.#..#.....#..#...#.#.##...#...##..#.#.###....#..###.....##..#.###..###.....##..###...#
..##.#####....##..#.#..##.##..######...#..###.######.....#..##...#.#..##..##..#..#..#..##.#.#.#.#...
.###.###.###...##...##..###..##.###.#.....##..##.#.#########...##..##.#..##.#..##..####..#.#.#.#####
#.#####..###.###.##.##.#...#.#.#.#..#.###...#..##.###.#...####.#..#.#.....###..#..####..#.#.#...##..
....#...##.....#....####.##.#.###..#.#.##..#.#...##.###.###..#.##..#.#.##..##..#.##.###..#.#.###.###
##.##...#.##...#.#..#.#..#...###...###.#..#..#.#####..###.#......#.....###.#####.#.#..#.#.#.##..#.#.
#.#..#.....#.....##.#..##...###..##...##...###.#.###.#..#.#.###...##..##..#.###...#.#######.#...#.#.
#.#.....####.#..#.##...#.##....#####.###.#.....#####....###..#........##..####...#...#.###....#..###
##.#.##..#.#.##.#.....##.#.....###.####.#..######.....####.#.#..##.#.##...#..#.#.....#.####.#.......
#..#..#.#..#.######.##..##.####.....##.#.##.#.######..#.#....#.#...#.#..#..#.#.###.#..#.#.#..#...###
####..####.#.#.###.....#.#.#.##..#.##.##.##.#..##..##.#.##.....#.#..#.####.....###.#..#.####.#.#..##
###.##..##.#.##..#..##...#.#####.##.#....##.####.#.##....#..###.#.#.##...#.....#.#.#.#.#..##.#.#..#.
......#..####...##.##...#.##.##...##..#..##.###..#...#..##...#.#....###.####...#.##.###.#.##.####.##
..#...#####.#.#..#.##....#..#...#..####.....###...##.###....#..#.###...#........#.#.##..#..#.#.....#
#######.#.#.###.###..######.##..#####.##.###.###....####.#..##.##...###.#..############.#.##....##.#
#.#...##.###.#.###..#.#.#.#.#.#..##..####.#..##.....#.##..#.##...##.#..##..#.#.#....##....##.#..#.#.
..#.#.####.....###..#######.#.#.#.#...##.#####.....##...##...##.###..######.###..#...####.#..###.###
.#.##....#.#.##..##.#.##.##..######...#.....#..#.#.#.#.....#.#..##.#.#.......#######....#.......#...
..###.##.##..##....#.###...#.....##..##......###...##..###.##...##.###.#.#.#.###.###.#.#...###..#...
.##.#.#...#...##.#.#...#..#..#.#...##.#.##...##..#....#.#..##.#..#.#..#.#.....#..#.#...#######.#.##.
...####....#.###.#..###..##...##..#.#.#.###...#..##.##.##..##.#...#..#.##.....#.#........#..#.#.####
.....##..###...#....#.#.#.#...###.###...#.#...#.#.####....#..####...###..#..######..##.##..###.#####
#####.##..#....###.###....##.....#.#..#....#.#####.##.#.####.#.##...#..###...###..##...#.###.#####..
###.##..........########.######....####.###.#..##...#.##.####.#.....##..#####..###...#####.....#.#.#
##..#####.##.#.#####.#.##.##..#.##....########.#####.#...#.###.##...#.###.#.#..#....##.#..#...#.#.#.
.##.#....#..#...#..#####..#..##.#......#..#....########...#..#...#.....####.#...##...#.###.#.#..##.#
.##.##.#.##.#.##...#.#.#..##.##.###.#..##..#...###.##.###.#####.#.###..#..###.#...#.###.#...#..#.#.#
.#..#..#.#..#..###..#....###.####.##.#.###.#.##.###.#.##.###.###...###...###.#...####...#.##.##.#.#.
###..##...###...#..##.#..#.#...##....###.##.##..#####....###..#..#....#..###.###.#...#.##...#.#.#..#
#....#.......##.....#.##...#..#.###.#.##..##..#.##..#.###..##.##...#####.#..#####..#####..#####....#
.####.####....###..###.#.##.####.##.#...####.#.###.#.....#...####..#####.###..#.#.###.##.##...##..#.
####..##...##.########...##..###..#..###.##.#.#.#........#.#####.#...#.###.####.#..####..#.#.#....##
###.#..#...###.#..#..#.###...##..###.##.#.#...#..#...####..##....#.#..#..##.#.#...#####.###.#..#.#.#
...##....#.###.#.#..##...##.###.#..#..#......#...#.#..####.#.##..######.####.#...#..#..#..##.#.#.##.
##.####.#...#..#.#.##..##.#.#.###..##...####......#..######.#......#.##.#....##...###.#.#..#......##
#.....#...#######.##.#..#.#...###.#..#.####....#.#.##.#.##...###..#...#.###.##..#.###..#.##...#####.
#####.##...#..#.#.#.......#.##..#####..#####...###..##.#.#..###.#.#####.####..#.#..##...#.##...#.###
.##.#..#######.###.#.####.....##...#.##.#.#..#...##....####......######.#..######.....##########.##.
##...#.#..#.##.###.#.#.#.##.###.##..##.##.##...#.#..###.#######..#.....#####..#....######.#..##..###
.#.#.###.....#..##..#.#..##..#.###...###.#..##...#...#.#####.#.#####..###.#..#...##..#.#..#..####...
.#......##..#.....####.###....##.###.....###.##........#.###.##..#..#.#######.#.######..##..###.....
..##.#.#..#.##...#.###.###...######..#..#.#..#....###.#.#....#..........#...##.##.##.#..##..#.#####.
###.###.#..#.##..##.#..#..##.....##.....#..#######.#..#.#.#.####.###..###.#.#..#.##.##.####.###.####
#.#.#..#....########.#..#..#...##..#.##..#.#..##..####...##.....#.##.#.#...########..#.###.#..#.#.##
.##.....#...#.#...##.##....###...##..#.####...#..#.#..#..#.##..#.###.##.####.##..####.....##.#.....#
....####.#.##.#.##.#..##.#.######.##.####..#...####.#..###.#.#..#..##.#.#.....##.#####.#.####...#.#.
#..#####.#####.....##....######..##....#..#.#.###.#####.....##.##.####.#...##...#.##.#.#####.##.#...
##.####..###.#....#...#.#.#.#.###.#####.#.####..####...####......##..#..#..#.#.##...########....#...
.###.#.#.#.#..####.##.#..######..#.#.###.....#.#......#.#.#.#..####.##...##.#####.#.##..##..#..#.#..
.....###...#...#.####.###.#.#.#.#.....#....#.####.###.##.##.##.#######......#.####......#....##.....
##..#..#.#.##..#...#..##.##.##..###.#....##.##....####.#.##.###....#.##.#.#.##...##.###...#..#..####
...#.#..##..##.#...##.##...#.#......#.#.##..###....####.##...#.#.###.#..#..#.####..##..##..#####.###
.##.##..##########.##...#.##.####.#.#######.##.#.##.##..#...##....########.###..##.##.##.#..##.#.#.#
#####.#....#.##..#.....#......##.##..#.##.###..##.......###..##.#.###.##.###....####.#..#.###..#.#.#
.#...#..#.##....##....#...####....#...#..#...####...########.###.#..##.#.#.##..###..#.#.###.....##.#
##..##.....###......#..###.##.####.##.####.#.#....#..#...#..#.#..#.###.#...#...#..##.##...#..#######
.....##..###..##...#####.#.#.....###.#.#..####...#.#.#..#..####..##.#..###.####.#....##..###....#..#
#.#.##.#....#.#####.#....##...#...##...##....#.#.......#....#..#...###.###.#.####..####....#.##.#.#.
..##...##..###.#.#.##.#..#....#.#.....##.###.#.###.###.....#...#.#..#######.#####..#.###...##......#
#......###..#....#.#..#.###.##.#...##..###.####.#.#....#.##..#.###..##.#..#####..##.###.....#..###..
##.#.##..##.###.#..##.....#.##.....###....##.####.######.#...#..###....#.#...#.##.....###....#..#.#.
.##.#.#.#.##..#.#.#..##..#.###.####....#..###.######..####.#.....###.##..#...###.#..######.##.#.##..
...##.####.#..##.#####.##.#...##..#..#...#.#.#.#####...#....#..###...#..#....#.#.##.#.######.#..####
..#.#.#.#...#.######.#.....#..#.#..###....#.#.########...#....#.#.##..#...##...#.#..#.#.###....##...
#####..#..##..#..##..#..#.#.##.#....#####.####.##.#.###..##..##....#.....#.#####.#...#.#####.##.#.#.
#.#..#####...####.###.###.....####.###.....##...##...#..#..#######.#.##....##..####.....##...#..#..#
#.#.###.#.#..##..#....#.#...#.#.##.##..#.##.....##...#.#..##.......##.#.###..#####.#.##....#.##.....
...#.......#....#.#.####.#.###.###..#....#..##.#..####........#.##..#...#.#...###.#..#.#.#...#...#..
...##.#####.##.#.###.##.##.#.##..##.#.#.#.#.#.##.#..##...##.#.#..#..##.##.#####.#.###...#####..#..#.
#######.#..#..#....##.#.#..####.#..#..###...#..#.......###.#.#.####....#.###...#.#.###.#.#.#.#..###.
..##.##.#.##.###....###.##.#.###.#...#....#.####..###..###.#.#..#...##.#.#.#..##.###..###.#.##...###
######..######..##..##.#.#.##.##.#..##..#.#.#.##..#.#...#...#.#.#..######.#..#.#.######..#......##.#
#.#####.....#.......#########..###.##...#...##.#.#..#...#####...#...#..#.###.#..#.#...###.#.#.#...#.
#....##....###...##.##.#...##.........##.#.#..#.#.##.#.######.#####..#..###.###.#...#.#.##.######...
#.#...###.#.###.##.#.######.#######.###.##..#.#.#...######.##.####.##..#.#.#.#......##..##.........#
..###..##....#.....##...#.#.###.#.#.....##.#...###.####.#...#...##..##.#.#.####..###...######....#.#
..###.#.##.####.#..#.##....##..#####....#..##.##.#..#######...#.####...##.#.#.##.........#....#....#
.##.#...#.####..#.#...#.##..######.##..##.#.###.##..###.###....##..#.##.##..##.#...###.##.##.###....
#...###.###.#..#....#.......#..#.....###..#.###.##.##....#.####.#.####.##..##..#..#.....#....##.#.#.
.##.#..#..#.##.......#.####.#######.....#.##.##.#.....#.#..#....######.#..###.##.##.....#.####..##.#
###..#.###.#..####.....##....#..####....#.##.##..#...######.#########...#.#....##...###.#..#.##...#.
#..###..##..#.#.##.###.#.#.##...###.#...##.##..#.###....###..#.#...#.###..######.#..#.###..#..#..#.#
.#........##.#.###..###.#.#.##.....##.##.#.#...##..#.##....###..#.#.#.#.##....#.##..#.#...###...#...
####.####..#....#.#.#..#..##.......##.####...###.##..#.#.##.#..##..######.......##.#.##..#...#.....#
..#..#..###..##.##..######.#..###..###.#.##..##.#..#####.#.#.#.##..#.##..##.##......####.#..........
...##.##..###.#...###....#.#.#.#.....#.##.....##...#...#......####...##.##....##.#..#.####.#..###.#.
..#.....####.#.###.#####..#..###..#..#.#...#####...###.###....#.###..#...#..#..#.#..#.##..##.#.#....
..##.#####...###.###.........#....##.####.##..#.#..#.#...#...##.##.##..#.#.##.########......#####...
...###.#.#..#...#.###.###.......##.###.#..#.##########...#..#.#.#.##.#.###...######..#.#...###.##...
.#.#.#######.#..##.##..##...#...####...#..#####.#..##...###.#.#...#.##...#......#..##.####..#.....##
.##.##.#.#......#######..###.....##.#.##..###......#....####...#.###.#.##.#........#..#....##.....##
#...#.###.#.##...##.####....#...#.###..#.#.....#.#....#.#.#.##...#.#..#####.#.#..#..#..#....#...####
.....##...###......#####..##.##.##...##.#.#####..##...#.#.#.#.###...###.##.####..#.#..#.#..#.####.##
#..#..##.#.##.#.##.#.#.#..###....###.##.#.##.#...#.#..#...#....###.#..#.#.######.#...####..#..##.#.#
#..#.#..#...###.#..##.#...#...##.#......#...#..#..####..##.....#.###...#.#..#.#....#.#####.##.###...
###....#.#..#.#..###..#.##......#...#..#..##.#..###..##..#..#.####..#...########..##.#.##.#.#.#...#.
.#.#.##.##.###..#...#.#....#..#.##..#.#.#.#.##.##.#####...#........####..###..####.#####..#.##.#.##."

# ╔═╡ de9c2f12-f7bb-4e0e-badc-1fffb75d569a
example_input = """
.#.#.#
...##.
#....#
..#...
#.#..#
####..
"""

# ╔═╡ 210f6c99-e328-44c4-b7bf-aa6714787310
function evolve(state)
	nextstate = []

	rows = split(strip(state), "\n")

	width = length(rows[1])
	height = length(rows)

	for (row_idx, row) in enumerate(rows)
		nextrow = ""

		for (col_idx, col) in enumerate(row)
			nsum = 0
			if col_idx > 1 # left
				nsum += row[col_idx-1] == '#' ? 1 : 0
			end
			if col_idx < width # right
				nsum += row[col_idx+1] == '#' ? 1 : 0
			end
			if row_idx > 1 # up
				nsum += rows[row_idx-1][col_idx] == '#' ? 1 : 0

				if col_idx > 1 # up left
					nsum += rows[row_idx-1][col_idx-1] == '#' ? 1 : 0
				end
				if col_idx < width # up right
					nsum += rows[row_idx-1][col_idx+1] == '#' ? 1 : 0
				end
			end
			if row_idx < height # down
				nsum += rows[row_idx+1][col_idx] == '#' ? 1 : 0

				if col_idx > 1 # down left
					nsum += rows[row_idx+1][col_idx-1] == '#' ? 1 : 0
				end
				if col_idx < width # down right
					nsum += rows[row_idx+1][col_idx+1] == '#' ? 1 : 0
				end
			end

			#m[height-row_idx+1, col_idx] = col == '#' ? 1 : 0
			# A light which is on stays on when 2 or 3 neighbors are on, and turns off otherwise.

			if col == '#'
				if nsum == 2 || nsum == 3
					nextrow *= '#'
				else
					nextrow *= '.'
				end
			else
			# A light which is off turns on if exactly 3 neighbors are on, and stays off otherwise.
				if nsum == 3
					nextrow *= '#'
				else
					nextrow *= '.'
				end
			end
		end
		push!(nextstate, nextrow)
	end
	join(nextstate, "\n")
end

# ╔═╡ add036cf-c9cf-4042-94e5-c786eecfa2fa
function repeatevolve(state, n)
	for i in 1:n
		state = evolve(state)
	end
	state
end

# ╔═╡ 8e0ff6c5-a830-4b14-b638-9ce8b0db85ad
lights(state) = count(f->f == '#', state)

# ╔═╡ 1bc1d59d-cf9f-4bb4-be4e-786b5864926d
@assert lights(repeatevolve(example_input, 4)) == 4

# ╔═╡ 9c5b6723-ae66-427c-919a-370ace69c63d
function grid2matrix(state)
	rows = split(strip(state), "\n")

	width = length(rows[1])
	height = length(rows)
	m = falses(height, width)

	for (row_idx, row) in enumerate(rows)
		for (col_idx, col) in enumerate(row)
			m[height-row_idx+1, col_idx] = col == '#' ? 1 : 0
		end
	end
	m
end

# ╔═╡ 3408aee3-3ee1-4fd9-926a-23f87ec96aa8
md"""
```
Initial state:
.#.#.#
...##.
#....#
..#...
#.#..#
####..
```
"""

# ╔═╡ 5d81ca32-cd4f-4164-a021-67c3f54b9370
heatmap(grid2matrix(repeatevolve(example_input, 0)))

# ╔═╡ d7ce4f93-e416-440d-9778-22898990edc7
md"""
```
After 1 step:
..##..
..##.#
...##.
......
#.....
#.##..
```"""

# ╔═╡ 642279fe-e0b8-4dfd-a015-f9f3e94c6e3a
heatmap(grid2matrix(repeatevolve(example_input, 1)))

# ╔═╡ 1afe4e36-f7df-4f67-a8ee-f639cd03d964
md"""
```
After 2 steps:
..###.
......
..###.
......
.#....
.#....
```"""

# ╔═╡ 6e4148ac-ba55-4dd4-bf3b-ce66bb28d2ee
heatmap(grid2matrix(repeatevolve(example_input, 2)))

# ╔═╡ 6e37530a-aba4-4520-a81e-28a47cf471f9
md"""
```
After 3 steps:
...#..
......
...#..
..##..
......
......
```"""

# ╔═╡ 65534039-fa37-486c-88dd-2c9cdae7f435
heatmap(grid2matrix(repeatevolve(example_input, 3)))

# ╔═╡ d7c592dc-fa4b-4e3d-81a0-cea70369e63c
md"""
```
After 4 steps:
......
......
..##..
..##..
......
......
```"""

# ╔═╡ 9c132f3d-4d90-46b9-b1e8-c055b0722fb7
heatmap(grid2matrix(repeatevolve(example_input, 4)))

# ╔═╡ 905eccfc-b5ec-4ffd-b055-33b16aa72f70
@gif for i=1:99
   heatmap(grid2matrix(repeatevolve(example_input, floor(i/20))))
end every 1

# ╔═╡ 0349d3fa-2eb2-43a4-8cc5-e995472eb3c9
heatmap(grid2matrix(repeatevolve(puzzle_input, 0)))

# ╔═╡ c9234ee3-9206-40fb-b398-571222598774
@gif for i=1:100
   heatmap(grid2matrix(repeatevolve(puzzle_input, i)))
end every 1

# ╔═╡ fd116e3b-2eda-4c68-bca4-152aaf1d05c0
heatmap(grid2matrix(repeatevolve(puzzle_input, 100)))

# ╔═╡ 214c1e6e-7035-4cdb-a380-d86919173535
part1 = lights(repeatevolve(puzzle_input, 100))

# ╔═╡ c44422ed-30e3-4ea7-a81a-291ec6e7084a
@test part1 == 814

# ╔═╡ c3ef9fe2-7a92-45dd-91fe-e31e132b62ec
md"Your puzzle answer was `814`."

# ╔═╡ 75f28fc6-29a1-4733-bd5b-88ac67c844b4
md"""
# Part Two

You flip the instructions over; Santa goes on to point out that this is all just an implementation of [Conway's Game of Life](https://en.wikipedia.org/wiki/Conway's_Game_of_Life). At least, it was, until you notice that something's wrong with the grid of lights you bought: four lights, one in each corner, are stuck on and can't be turned off. The example above will actually run like this:

```
Initial state:
##.#.#
...##.
#....#
..#...
#.#..#
####.#

After 1 step:
#.##.#
####.#
...##.
......
#...#.
#.####

After 2 steps:
#..#.#
#....#
.#.##.
...##.
.#..##
##.###

After 3 steps:
#...##
####.#
..##.#
......
##....
####.#

After 4 steps:
#.####
#....#
...#..
.##...
#.....
#.#..#

After 5 steps:
##.###
.##..#
.##...
.##...
#.#...
##...#

```

After `5` steps, this example now has `17` lights on.

In your grid of 100x100 lights, given your initial configuration, but with the four corners always in the on state, **how many lights are on after 100 steps?**
"""

# ╔═╡ 72976f12-9139-4d89-8329-ec57b9535410
function fixrowcorners(row)
	chars = collect(row)
	chars[1] = '#'
	chars[length(chars)] = '#'
	join(chars)
end

# ╔═╡ 385a7cab-ea88-4ad6-b614-6d446549eef0
function fixcorners(state)
	rows = split(strip(state), "\n")

	rows[1] = fixrowcorners(rows[1])
	rows[length(rows)] = fixrowcorners(rows[length(rows)])

	join(rows, "\n")
end

# ╔═╡ bac289d5-d7f5-4fa9-a512-72d24d44d1a2
function repeatevolve2(state, n)
	for i in 1:n
		state = evolve(fixcorners(state))
	end
	fixcorners(state)
end

# ╔═╡ 6a36690d-682f-43c7-8a65-34a6b3b5da1a
@assert lights(repeatevolve2(strip(example_input), 5)) == 17

# ╔═╡ f87a84a5-a5cc-421f-9065-a846ca5b9319
md"""

```
Initial state:
##.#.#
...##.
#....#
..#...
#.#..#
####.#
```

"""

# ╔═╡ 259aabd1-703c-45d2-b93e-3f7782b09e58
heatmap(grid2matrix(repeatevolve2(strip(example_input), 0)))

# ╔═╡ 84f4f18a-665a-4667-8733-edcd479d9028
md"""

```

After 1 step:
#.##.#
####.#
...##.
......
#...#.
#.####

```

"""

# ╔═╡ c8752e7b-a4d1-4cb1-bfb3-b5c131d6dec7
heatmap(grid2matrix(repeatevolve2(strip(example_input), 1)))

# ╔═╡ f85bce95-8a94-4658-91fa-aa96dc1e86ca
md"""

```

After 2 steps:
#..#.#
#....#
.#.##.
...##.
.#..##
##.###

```

"""

# ╔═╡ 22a54b94-9e5a-4d5f-924b-fc4a2d585224
heatmap(grid2matrix(repeatevolve2(strip(example_input), 2)))

# ╔═╡ eff70454-34fc-435c-9d8e-d0da4cf440a1
md"""

```

After 3 steps:
#...##
####.#
..##.#
......
##....
####.#

```

"""

# ╔═╡ 8c17a7ae-f5a0-439a-b8ac-7ab106b88988
heatmap(grid2matrix(repeatevolve2(strip(example_input), 3)))

# ╔═╡ cbcdafd8-a2ab-406d-b0a7-cc154f2beb8c
md"""

```

After 4 steps:
#.####
#....#
...#..
.##...
#.....
#.#..#

```

"""

# ╔═╡ e13011e3-5441-419c-b01b-1810dafc5166
heatmap(grid2matrix(repeatevolve2(strip(example_input), 4)))

# ╔═╡ b5342e62-14db-4597-82a9-19d1662cf02c
md"""

```

After 5 steps:
##.###
.##..#
.##...
.##...
#.#...
##...#

```

"""

# ╔═╡ 6ceea610-1318-4458-98cf-aba5a40b3ead
heatmap(grid2matrix(repeatevolve2(strip(example_input), 5)))

# ╔═╡ a6fc1934-51d3-40e7-9098-282f18edc10f
@gif for i=1:100
   heatmap(grid2matrix(repeatevolve2(puzzle_input, i)))
end every 1

# ╔═╡ 6cc2df3c-7532-4b5b-b306-24bc6b06000f
heatmap(grid2matrix(repeatevolve2(puzzle_input, 100)))

# ╔═╡ 038dca7c-0369-41f4-8902-56ba5081b0fc
part2 = lights(repeatevolve2(puzzle_input, 100))

# ╔═╡ 73ce4e85-0913-4ae8-97dd-dca6d09f2ac8
@test part2 == 924

# ╔═╡ 2e0eb85c-f975-4760-bf4e-e4319d5f1cce
md"Your puzzle answer was `924`."

# ╔═╡ Cell order:
# ╠═6ca83e83-423a-4621-ae14-dc8da5913394
# ╠═df80bc4c-770a-4acd-a4f3-337c6b720c71
# ╟─20fb22d0-b0e9-11eb-1561-8b6505f5cb45
# ╟─c73f4666-7302-4743-bb13-07ca4aea6440
# ╠═de9c2f12-f7bb-4e0e-badc-1fffb75d569a
# ╠═210f6c99-e328-44c4-b7bf-aa6714787310
# ╠═add036cf-c9cf-4042-94e5-c786eecfa2fa
# ╠═8e0ff6c5-a830-4b14-b638-9ce8b0db85ad
# ╠═1bc1d59d-cf9f-4bb4-be4e-786b5864926d
# ╠═9c5b6723-ae66-427c-919a-370ace69c63d
# ╟─3408aee3-3ee1-4fd9-926a-23f87ec96aa8
# ╠═5d81ca32-cd4f-4164-a021-67c3f54b9370
# ╟─d7ce4f93-e416-440d-9778-22898990edc7
# ╠═642279fe-e0b8-4dfd-a015-f9f3e94c6e3a
# ╟─1afe4e36-f7df-4f67-a8ee-f639cd03d964
# ╠═6e4148ac-ba55-4dd4-bf3b-ce66bb28d2ee
# ╟─6e37530a-aba4-4520-a81e-28a47cf471f9
# ╠═65534039-fa37-486c-88dd-2c9cdae7f435
# ╟─d7c592dc-fa4b-4e3d-81a0-cea70369e63c
# ╠═9c132f3d-4d90-46b9-b1e8-c055b0722fb7
# ╠═905eccfc-b5ec-4ffd-b055-33b16aa72f70
# ╠═0349d3fa-2eb2-43a4-8cc5-e995472eb3c9
# ╠═c9234ee3-9206-40fb-b398-571222598774
# ╠═fd116e3b-2eda-4c68-bca4-152aaf1d05c0
# ╠═214c1e6e-7035-4cdb-a380-d86919173535
# ╠═c44422ed-30e3-4ea7-a81a-291ec6e7084a
# ╟─c3ef9fe2-7a92-45dd-91fe-e31e132b62ec
# ╟─75f28fc6-29a1-4733-bd5b-88ac67c844b4
# ╠═385a7cab-ea88-4ad6-b614-6d446549eef0
# ╠═72976f12-9139-4d89-8329-ec57b9535410
# ╠═bac289d5-d7f5-4fa9-a512-72d24d44d1a2
# ╠═6a36690d-682f-43c7-8a65-34a6b3b5da1a
# ╟─f87a84a5-a5cc-421f-9065-a846ca5b9319
# ╠═259aabd1-703c-45d2-b93e-3f7782b09e58
# ╟─84f4f18a-665a-4667-8733-edcd479d9028
# ╠═c8752e7b-a4d1-4cb1-bfb3-b5c131d6dec7
# ╟─f85bce95-8a94-4658-91fa-aa96dc1e86ca
# ╠═22a54b94-9e5a-4d5f-924b-fc4a2d585224
# ╟─eff70454-34fc-435c-9d8e-d0da4cf440a1
# ╠═8c17a7ae-f5a0-439a-b8ac-7ab106b88988
# ╟─cbcdafd8-a2ab-406d-b0a7-cc154f2beb8c
# ╠═e13011e3-5441-419c-b01b-1810dafc5166
# ╟─b5342e62-14db-4597-82a9-19d1662cf02c
# ╠═6ceea610-1318-4458-98cf-aba5a40b3ead
# ╠═a6fc1934-51d3-40e7-9098-282f18edc10f
# ╠═6cc2df3c-7532-4b5b-b306-24bc6b06000f
# ╠═038dca7c-0369-41f4-8902-56ba5081b0fc
# ╠═73ce4e85-0913-4ae8-97dd-dca6d09f2ac8
# ╟─2e0eb85c-f975-4760-bf4e-e4319d5f1cce


text:
import Base: -, +, zero,
             <, >, <=, >=

"""
    zero(::Type{Point2D})
    zero(p::Point2D)
Gives origin for a point of a particular type
"""
zero(::Type{Point{T}}) where T <: Number = Point{T}(zero(T), zero(T))
zero(p::Point{T})      where T <: Number = zero(typeof(p))

-(p::Point, q::Point) = Vec2D(p.x - q.x, p.y - q.y)
-(p::Point, v::Vec2D) = Point(p.x - v.dx, p.y - v.dy)
+(p::Point, v::Vec2D) = Point(p.x + v.dx, p.y + v.dy)
+(v::Vec2D, p::Point) = p + v

function <(p::Point, q::Point)
    p.x < q.x || (p.x == q.x && p.y < q.y)
end

function >(p::Point, q::Point)
    p.x > q.x || (p.x == q.x && p.y > q.y)
end

<=(p::Point, q::Point) = p < q || p == q
>=(p::Point, q::Point) = p > q || p == q

ext: jl


text:
include("utils.jl")

struct AcrobotDrawParams <: AbstractDrawParams end

function Ctx(env::AcrobotEnv, mode::Symbol=:no_render)
    if mode == :human_pane
        draw_parameters = AcrobotDrawParams()
        viewer = CairoRGBSurface(500, 500)  # screen_size = 500 x 500

        CairoCtx(draw_parameters, viewer)
    elseif mode == :human_window
        draw_parameters = AcrobotDrawParams()
        viewer = CairoRGBSurface(draw_parameters.screen_width, draw_parameters.screen_height)

        canvas = @GtkCanvas()
        canvas.backcc = CairoContext(viewer)
        win = GtkWindow(canvas, "Acrobot",
                draw_parameters.screen_width, draw_parameters.screen_height; resizable=false)
        show(canvas)
        visible(win, false)
        signal_connect(win, "delete-event") do widget, event
            ccall((:gtk_widget_hide_on_delete, Gtk.libgtk), Bool, (Ptr{GObject},), win)
        end

        GtkCtx(draw_parameters, canvas, win)
    elseif mode == :rgb
        draw_parameters = AcrobotDrawParams()
        viewer = CairoRGBSurface(draw_parameters.screen_width, draw_parameters.screen_height)

        RGBCtx(draw_parameters, viewer)
    elseif mode == :no_render
        NoCtx()
    else
        error("Unrecognized mode in Ctx(): $(mode)")
    end
end

function drawcanvas!(env::AcrobotEnv, viewer::CairoContext, parameters::AcrobotDrawParams)
    # Set background
    set_source_rgb(viewer, 1, 1, 1)
    rectangle(viewer, 0, 0, 500, 500) # screen_size = 500 x 500
    fill(viewer)

    # Draw goal
    set_source_rgb(viewer, 0, 0, 0)
    move_to(viewer, 0, 136)
    line_to(viewer, 500, 136)
    stroke(viewer)
    close_path(viewer)

    # First translation and rotation
    translate(viewer, 250, 250)  # Bring the origin to the center of the screen
    theta1 = -env.state[1].data
    rotate(viewer, theta1)

    # Draw the first link
    set_source_rgb(viewer, 0, 0.8, 0.8)
    set_line_width(viewer, 19)
    move_to(viewer, 0, 0)
    line_to(viewer, 0, 113)  # Link length translated into the view coordinates
    stroke(viewer)
    set_line_width(viewer, 1)  # Reset line width to avoid unexpected surprises

    # Draw the first joint
    set_source_rgb(viewer, 0.8, 0.8, 0)
    circle(viewer, 0, 0, 10)
    fill(viewer)

    # Second translation and rotation
    translate(viewer, 0, 113)
    theta2 = -env.state[2].data
    rotate(viewer, theta2)

    # Draw the second link
    set_source_rgb(viewer, 0, 0.8, 0.8)
    set_line_width(viewer, 19)
    move_to(viewer, 0, 0)
    line_to(viewer, 0, 113)
    stroke(viewer)
    close_path(viewer)
    set_line_width(viewer, 1)

    # Draw the second joint
    set_source_rgb(viewer, 0.8, 0.8, 0)
    circle(viewer, 0, 0, 10)
    fill(viewer)

    # Undo all translations and rotations
    rotate(viewer, -theta2)
    translate(viewer, 0, -113)
    rotate(viewer, -theta1)
    translate(viewer, -250, -250)

    #=
    # Calibration circle
    set_source_rgb(viewer, 0, 0, 0)
    calibration_radius = 10
    circle(viewer, calibration_radius, calibration_radius, calibration_radius)
    set_line_width(viewer, 1)
    stroke(viewer)
    move_to(viewer, calibration_radius, calibration_radius)
    line_to(viewer, 2calibration_radius, calibration_radius)
    close_path(viewer)
    set_line_width(viewer, 1)
    stroke(viewer)
    set_line_width(viewer, 1)
    =#
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function runcode(app)

print("run code...\n");

cd(app.sourcepath * "C++/Main");
run(Gencode.string2cmd("./cpu" * app.codename * " " * app.appname * " out"));
cd(app.currentdir);

# if numpde==1
#     mystr = app.appname * " 1 datain/ dataout/out";
# else
#     mystr = app.appname * " " * string(numpde);
#     for i = 1:numpde
#         mystr = mystr * " datain" * string(i) * "/ " * "dataout" * string(i) * "/out";
#     end
# end

# mpirun = app.mpirun;
# if app.platform == "cpu"
#     if app.mpiprocs==1
#         runstr = "./app/serial" * mystr;
#     else
#         runstr = mpirun * " -np " * string(app.mpiprocs) * " ./app/mpi" * mystr;
#     end
#     run(string2cmd(runstr));
# elseif app.platform == "gpu"
#     if app.mpiprocs==1
#         runstr = "./app/gpu" * mystr;
#     else
#         runstr = mpirun * " -np " * string(app.mpiprocs) * " ./app/gpumpi" * mystr;
#     end
#     run(string2cmd(runstr));
# else
#     error("app.platform must be either cpu or gpu");
# end

# return runstr

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Base.Test
using MeshCat
using GeometryTypes
using CoordinateTransformations
using Colors
using MeshIO, FileIO
using Polyhedra
using CDDLib

@testset "MeshCat" begin
    include("server_client.jl")
    include("visualizer.jl")
    include("notebook.jl")

    @testset "deprecations" begin
        # Deprecated in v0.0.2
        @testset "cylinder" begin
            l = 1.0
            r = 2.0
            c = HyperCylinder(l, r)
            @test c.origin == Point(0., 0, 0)
            @test c.extremity == Point(0., 0, l)
            @test radius(c) == r
        end
    end
end


text:
@defcomp impactvectorbornediseases begin
    regions = Index()

    dengue = Variable(index=[time,regions])
    schisto = Variable(index=[time,regions])
    malaria = Variable(index=[time,regions])

    dfbs = Parameter(index=[regions])
    dfch = Parameter(index=[regions])
    smbs = Parameter(index=[regions])
    smch = Parameter(index=[regions])
    malbs = Parameter(index=[regions])
    malch = Parameter(index=[regions])

    gdp90 = Parameter(index=[regions])
    pop90 = Parameter(index=[regions])

    income = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])
    temp = Parameter(index=[time,regions])

    dfnl = Parameter(default = 1)
    vbel = Parameter(default = -2.65)
    smnl = Parameter(default = 1)
    malnl = Parameter(default = 1)

    function run_timestep(p, v, d, t)

        for r in d.regions
            ypc = 1000.0 * p.income[t, r] / p.population[t, r]
            ypc90 = p.gdp90[r] / p.pop90[r] * 1000.0

            v.dengue[t, r] = p.dfbs[r] * p.population[t, r] * p.dfch[r] * p.temp[t, r]^p.dfnl * (ypc / ypc90)^p.vbel

            v.schisto[t, r] = p.smbs[r] * p.population[t, r] * p.smch[r] * p.temp[t, r]^p.smnl * (ypc / ypc90)^p.vbel

            if v.schisto[t, r] < -p.smbs[r] * p.population[t, r] * (ypc / ypc90)^p.vbel
                v.schisto[t, r] = -p.smbs[r] * p.population[t, r] * (ypc / ypc90)^p.vbel
            end

            v.malaria[t, r] = p.malbs[r] * p.population[t, r] * p.malch[r] * p.temp[t, r]^p.malnl * (ypc / ypc90)^p.vbel
        end
    end
end


text:
@show PROGRAM_FILE
@show @__DIR__
@show abspath(PROGRAM_FILE) == @__FILE__

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
filter_with(df::DataFrame, set::Any; kwargs...)

# Arguments
- `df::DataFrame` to filter.
- `set::Dict` or `set::NamedTuple`: Values to keep in the DataFrame.

# Keywords
- `extrapolate::Bool = false`: Add missing regions/years to the DataFrame?
    If `extrapolate` is set to true, the following `kwargs` become relevant:
    - When extrapolating over years,
        - `backward::Bool = true`: Do we extrapolate backward in time?
        - `forward::Bool = true`: Do we extrapolate forward in time?
        Currently, "extrapolating" means copying the closest
    - When extrapolating across regions,
        - `r::Pair = "md" => "dc"`: `Pair` indicating a region (`r.first`) to extrapolate to
            another region (`r.second`). A suggested regional extrapolation: MD data will be
            used to approximate DC data in the event that it is missing.
        - `overwrite::Bool = false`: If data in the target region `r.second` is already present,
            should it be overwritten?

# Returns
- `df::DataFrame` with only the desired keys.

# Examples

```jldoctest filter_with
julia> df = read_file(joinpath(SLIDE_DIR,"docs","src","assets","data","filter_use.csv"))
14×4 DataFrame
│ Row │ yr    │ i      │ j      │ value   │
│     │ Int64 │ String │ String │ Float64 │
├─────┼───────┼────────┼────────┼─────────┤
│ 1   │ 2015  │ agr    │ agr    │ 69.42   │
│ 2   │ 2015  │ agr    │ fbp    │ 277.179 │
│ 3   │ 2015  │ fbp    │ agr    │ 49.132  │
│ 4   │ 2015  │ fbp    │ fbp    │ 210.998 │
│ 5   │ 2015  │ uti    │ agr    │ 4.846   │
│ 6   │ 2015  │ uti    │ fbp    │ 10.102  │
│ 7   │ 2015  │ uti    │ uti    │ 35.093  │
│ 8   │ 2016  │ agr    │ agr    │ 60.197  │
│ 9   │ 2016  │ agr    │ fbp    │ 264.173 │
│ 10  │ 2016  │ fbp    │ agr    │ 47.739  │
│ 11  │ 2016  │ fbp    │ fbp    │ 205.21  │
│ 12  │ 2016  │ uti    │ agr    │ 4.548   │
│ 13  │ 2016  │ uti    │ fbp    │ 9.152   │
│ 14  │ 2016  │ uti    │ uti    │ 27.47   │

julia> df = filter_with(df, (i = ["agr","fbp"], j = ["agr","fbp"]))
8×4 DataFrame
│ Row │ yr    │ i      │ j      │ value   │
│     │ Int64 │ String │ String │ Float64 │
├─────┼───────┼────────┼────────┼─────────┤
│ 1   │ 2015  │ agr    │ agr    │ 69.42   │
│ 2   │ 2015  │ agr    │ fbp    │ 277.179 │
│ 3   │ 2015  │ fbp    │ agr    │ 49.132  │
│ 4   │ 2015  │ fbp    │ fbp    │ 210.998 │
│ 5   │ 2016  │ agr    │ agr    │ 60.197  │
│ 6   │ 2016  │ agr    │ fbp    │ 264.173 │
│ 7   │ 2016  │ fbp    │ agr    │ 47.739  │
│ 8   │ 2016  │ fbp    │ fbp    │ 205.21  │

julia> filter_with(df, (yr = 2016,); drop = true)
4×3 DataFrame
│ Row │ i      │ j      │ value   │
│     │ String │ String │ Float64 │
├─────┼────────┼────────┼─────────┤
│ 1   │ agr    │ agr    │ 60.197  │
│ 2   │ agr    │ fbp    │ 264.173 │
│ 3   │ fbp    │ agr    │ 47.739  │
│ 4   │ fbp    │ fbp    │ 205.21  │
```
"""
function filter_with(
    df::DataFrame,
    set::Any;
    drop = false,
    extrapolate::Bool = false,
    forward::Bool = true,
    backward::Bool = true,
    r::Pair = "md" => "dc",
    overwrite::Bool = false
)
    cols = propertynames(df)

    # Find keys that reference both column names in the input DataFrame df and
    # values in the set Dictionary. Then, created a DataFrame containing all permutations.
    df_set = _intersect_with(set, df; intersect_values=true)

    if df_set === nothing
        @warn("Returning unfiltered DataFrame.")
        return df
    end

    # Drop values that are not in the current set.
    idx_set = propertynames(df_set)
    df = innerjoin(df, df_set, on = idx_set)

    if extrapolate
        :yr in idx_set && (df = extrapolate_year(df, set; forward = forward, backward = backward))
        :r in idx_set  && (df = extrapolate_region(df, r; overwrite = overwrite))
    end

    # If one of the filtered DataFrame columns contains only one unique value, drop it.
    # If drop specifies columns to drop, drop only these IFF they contain one unique value.
    # NEVER drop units and never drop columns containing multiple unique values.
    if drop !== false
        idx_drop = setdiff(_find_constant(df[:,idx_set]), propertynames_with(df,:units))
        drop !== true && intersect!(idx_drop, ensurearray(drop))
        setdiff!(cols, idx_drop)
    end

    return sort(df[:,cols])
end

"""
    extrapolate_year(df::DataFrame, yr::Array{Int64,1}; kwargs...)
    extrapolate_year(df::DataFrame, set::Any; kwargs...)

# Arguments
- `df::DataFrame` that might be in need of extrapolation.
- `yr::Array{Int64,1}`: List of years overwhich extrapolation is possible (depending on the kwargs)
- `set::Dict` or `set::NamedTuple` containing list of years, identified by the key `:yr`.

# Keywords
- `backward::Bool = true`: Do we extrapolate backward in time?
- `forward::Bool = true`: Do we extrapolate forward in time?

# Returns
- `df::DataFrame` extrapolated in time.

# Example
Continuing with the DataFrame from [`SLiDE.filter_with`](@ref),

```jldoctest extrapolate_year; setup = :(df = filter_with(read_file(joinpath(SLIDE_DIR,"docs","src","assets","data","filter_use.csv")), (i = ["agr","fbp"], j = ["agr","fbp"])))
julia> df
8×4 DataFrame
│ Row │ yr    │ i      │ j      │ value   │
│     │ Int64 │ String │ String │ Float64 │
├─────┼───────┼────────┼────────┼─────────┤
│ 1   │ 2015  │ agr    │ agr    │ 69.42   │
│ 2   │ 2015  │ agr    │ fbp    │ 277.179 │
│ 3   │ 2015  │ fbp    │ agr    │ 49.132  │
│ 4   │ 2015  │ fbp    │ fbp    │ 210.998 │
│ 5   │ 2016  │ agr    │ agr    │ 60.197  │
│ 6   │ 2016  │ agr    │ fbp    │ 264.173 │
│ 7   │ 2016  │ fbp    │ agr    │ 47.739  │
│ 8   │ 2016  │ fbp    │ fbp    │ 205.21  │

julia> extrapolate_year(df, Dict(:yr => 2014:2017))
16×4 DataFrame
│ Row │ yr    │ i      │ j      │ value   │
│     │ Int64 │ String │ String │ Float64 │
├─────┼───────┼────────┼────────┼─────────┤
│ 1   │ 2014  │ agr    │ agr    │ 69.42   │
│ 2   │ 2014  │ agr    │ fbp    │ 277.179 │
│ 3   │ 2014  │ fbp    │ agr    │ 49.132  │
│ 4   │ 2014  │ fbp    │ fbp    │ 210.998 │
│ 5   │ 2015  │ agr    │ agr    │ 69.42   │
│ 6   │ 2015  │ agr    │ fbp    │ 277.179 │
│ 7   │ 2015  │ fbp    │ agr    │ 49.132  │
│ 8   │ 2015  │ fbp    │ fbp    │ 210.998 │
│ 9   │ 2016  │ agr    │ agr    │ 60.197  │
│ 10  │ 2016  │ agr    │ fbp    │ 264.173 │
│ 11  │ 2016  │ fbp    │ agr    │ 47.739  │
│ 12  │ 2016  │ fbp    │ fbp    │ 205.21  │
│ 13  │ 2017  │ agr    │ agr    │ 60.197  │
│ 14  │ 2017  │ agr    │ fbp    │ 264.173 │
│ 15  │ 2017  │ fbp    │ agr    │ 47.739  │
│ 16  │ 2017  │ fbp    │ fbp    │ 205.21  │
```
"""
function extrapolate_year(
    df::DataFrame,
    yr::Array{Int64,1};
    backward::Bool = true,
    forward::Bool = true
)
    df = copy(df)
    yr_diff = setdiff(yr, unique(df[:,:yr]))
    length(yr_diff) == 0 && (return df)

    cols = setdiff(propertynames(df), [:yr])
    cols_ans = propertynames(df)

    df_ext = []

    if backward
        yr_min = minimum(df[:,:yr])
        df_min = filter_with(df, (yr = yr_min,))[:,cols]

        yr_back = yr_diff[yr_diff .< yr_min]
        df_back = crossjoin(DataFrame(yr = yr_back), df_min)[:,cols_ans]

        push!(df_ext, df_back)
    end

    if forward
        yr_max = maximum(df[:,:yr])
        df_max = filter_with(df, (yr = yr_max,))[:,cols]

        yr_forward = yr_diff[yr_diff .> yr_max]
        df_forward = crossjoin(DataFrame(yr = yr_forward), df_max)[:,cols_ans]

        push!(df_ext, df_forward)
    end
    return sort([df_ext...; df])
end

function extrapolate_year(
    df::DataFrame,
    set;
    backward::Bool = true,
    forward::Bool = true
)
    extrapolate_year(df, set[:yr]; forward = forward, backward = backward)
end

function extrapolate_year(
    df::DataFrame,
    yr::UnitRange{Int64};
    backward::Bool = true,
    forward::Bool = true
)
    extrapolate_year(df, ensurearray(yr); forward = forward, backward = backward)
end

"""
    extrapolate_region(df::DataFrame; kwargs...)
    extrapolate_region(df::DataFrame, r::Pair; kwargs...)

Fills in missing data in the input DataFrame `df` by filling it with existing information in
`df`. Here, "extrapolate" makes a direct copy of the data.

# Arguments
- `df::DataFrame` that might be in need of extrapolation.
- `r::Pair = "md" => "dc"`: `Pair` indicating a region (`r.first`) to extrapolate to another
    region (`r.second`). A suggested regional extrapolation: MD data will be used to
    approximate DC data in the event that it is missing. To fill multiple regions with data,
    use "md" => ["dc","va"].

# Keyword Argument:
- `overwrite::Bool = false`: If data in the target region `r.second` is already present,
    should it be overwritten?

# Returns
- `df::DataFrame` extrapolated in region.

# Example

```jldoctest extrapolate_region
julia> df = read_file(joinpath(SLIDE_DIR,"docs","src","assets","data","filter_utd.csv"))
8×5 DataFrame
│ Row │ yr    │ r      │ s      │ t       │ value     │
│     │ Int64 │ String │ String │ String  │ Float64   │
├─────┼───────┼────────┼────────┼─────────┼───────────┤
│ 1   │ 2015  │ md     │ agr    │ exports │ 0.0390152 │
│ 2   │ 2015  │ md     │ agr    │ imports │ 0.778159  │
│ 3   │ 2015  │ va     │ agr    │ exports │ 1.11601   │
│ 4   │ 2015  │ va     │ agr    │ imports │ 0.88253   │
│ 5   │ 2016  │ md     │ agr    │ exports │ 0.0330508 │
│ 6   │ 2016  │ md     │ agr    │ imports │ 0.762089  │
│ 7   │ 2016  │ va     │ agr    │ exports │ 1.16253   │
│ 8   │ 2016  │ va     │ agr    │ imports │ 0.86741   │

julia> extrapolate_region(df)
12×5 DataFrame
│ Row │ r      │ yr    │ s      │ t       │ value     │
│     │ String │ Int64 │ String │ String  │ Float64   │
├─────┼────────┼───────┼────────┼─────────┼───────────┤
│ 1   │ dc     │ 2015  │ agr    │ exports │ 0.0390152 │
│ 2   │ dc     │ 2015  │ agr    │ imports │ 0.778159  │
│ 3   │ dc     │ 2016  │ agr    │ exports │ 0.0330508 │
│ 4   │ dc     │ 2016  │ agr    │ imports │ 0.762089  │
│ 5   │ md     │ 2015  │ agr    │ exports │ 0.0390152 │
│ 6   │ md     │ 2015  │ agr    │ imports │ 0.778159  │
│ 7   │ md     │ 2016  │ agr    │ exports │ 0.0330508 │
│ 8   │ md     │ 2016  │ agr    │ imports │ 0.762089  │
│ 9   │ va     │ 2015  │ agr    │ exports │ 1.11601   │
│ 10  │ va     │ 2015  │ agr    │ imports │ 0.88253   │
│ 11  │ va     │ 2016  │ agr    │ exports │ 1.16253   │
│ 12  │ va     │ 2016  │ agr    │ imports │ 0.86741   │
```

If we instead want to copy VA data into DC, specify:

```jldoctest extrapolate_region
julia> extrapolate_region(df, "va" => "dc")
12×5 DataFrame
│ Row │ r      │ yr    │ s      │ t       │ value     │
│     │ String │ Int64 │ String │ String  │ Float64   │
├─────┼────────┼───────┼────────┼─────────┼───────────┤
│ 1   │ dc     │ 2015  │ agr    │ exports │ 1.11601   │
│ 2   │ dc     │ 2015  │ agr    │ imports │ 0.88253   │
│ 3   │ dc     │ 2016  │ agr    │ exports │ 1.16253   │
│ 4   │ dc     │ 2016  │ agr    │ imports │ 0.86741   │
│ 5   │ md     │ 2015  │ agr    │ exports │ 0.0390152 │
│ 6   │ md     │ 2015  │ agr    │ imports │ 0.778159  │
│ 7   │ md     │ 2016  │ agr    │ exports │ 0.0330508 │
│ 8   │ md     │ 2016  │ agr    │ imports │ 0.762089  │
│ 9   │ va     │ 2015  │ agr    │ exports │ 1.11601   │
│ 10  │ va     │ 2015  │ agr    │ imports │ 0.88253   │
│ 11  │ va     │ 2016  │ agr    │ exports │ 1.16253   │
│ 12  │ va     │ 2016  │ agr    │ imports │ 0.86741   │
```
"""
function extrapolate_region(df::DataFrame, r::Pair = "md" => "dc"; overwrite = false)
    df = copy(df)
    if !overwrite
        r = r.first => setdiff(ensurearray(r.second), unique(df[:,:r]))
        length(r.second) == 0 && (return df)
    else
        df = edit_with(df, Drop.(:r, r.second, "=="))
    end

    cols = setdiff(propertynames(df), [:r])
    df_close = crossjoin(DataFrame(r = r.second), filter_with(df, (r = r.first,))[:,cols])

    return sort([df_close; df])
end

"""
"""
function _intersect_with(
    x::Union{Dict,NamedTuple},
    df::DataFrame;
    intersect_values::Bool=false
)
    T = typeof(x)
    idx = intersect(findindex(df), collect(keys(x)))

    if isempty(idx)
        @error("Cannot overlap input $T and DataFrame. No overlapping keys/columns.")
        return nothing
    end

    val = [intersect_values ? intersect(unique(df[:,k]), ensurearray(x[k])) : x[k] for k in idx]

    ii = .!isempty.(val)

    if any(.!ii)
        idx_drop = idx[.!ii]
        @warn("No overlapping values in $idx_drop. These will be removed.")
        idx = idx[ii]
        val = val[ii]

        if isempty(idx)
            @error("Cannot overlap input $T and DataFrame. No overlapping values.")
            return nothing
        end
    end

    return DataFrame(permute(NamedTuple{Tuple(idx,)}(val,)))
end

"""
"""
function _find_constant(df::DataFrame)
    idx = findindex(df)
    return idx[length.(unique.(eachcol(df[:,idx]))) .== 1]
end

"""
    dropzero!(df::DataFrame, x::Float64)
Removes values `x` in columns of type AbstractFloat from `df`.
"""
function dropvalue!(df::DataFrame, x::Float64)
    cols = find_oftype(df, typeof(x));
    if isnan(x); [filter!(row -> .!isnan.(row[col]), df) for col in cols]
    else;        [filter!(row -> row[col] .!== x, df) for col in cols]
    end
    return df
end

dropvalue(df::DataFrame, x::Float64) = dropvalue!(copy(df), x)

"""
    dropzero!(df::DataFrame)
Removes zero values in columns of type AbstractFloat from `df`
"""
dropzero!(df::DataFrame) = dropvalue!(dropvalue!(df, 0.0), -0.0)
dropzero(df::DataFrame) = dropzero!(copy(df))

"""
    dropnan!(df::DataFrame)
Removes `NaN` values in columns of type AbstractFloat from `df`
"""
dropnan!(df::DataFrame) = dropvalue!(df, NaN)
dropnan(df::DataFrame) = dropnan!(copy(df))

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
iterations = 1000

pop = [6000 1000 1000 0 0 0 0
2000 5000 1000 0 10 10 0
500 3000 12000 0 30 0 0
0 0 10 0 0 2 0
2 10 20 10 3 20 0
5 2 0 0 0 10000 2000
3 1 0 0 1000 1000 5000].*1000

road = [11 122122
12 11132122
13 1223
21 1112223331
22 111221233233
23 13223233
31 2132
32 31223323
33 3221321223
35 25
25 35
43 5354
53 4354
52 61
61 5271
75 76
76 666777
66 766777
67 667677
77 666776]

aeroplane = [11 33
33 1177
77 33]

port = [33 43
33 35
43 33
35 33
53 52
52 53
25 26
26 25
26 46
46 26
34 56
56 34
56 66
66 56]

L = vec([.6 .5 .5 0 0 0 0
.5 .6 .6 0 .4 .3 0
.5 .6 .7 0 .3 0 0
0 0 .5 0 0 .3 0
.4 .5 .3 .2 .4 .5 0
.5 .3 0 0 0 .9 .8
.6 .4 0 0 .8 .8 .9]')'

beta = 0.7

B = vec([10000 1000 1000 0 0 0 0
1500 8000 1000 0 200 200 0
500 2000 20000 0 200 0 0
0 0 200 0 0 20 0
40 200 400 200 30 300 0
80 40 0 0 0 10000 1200
50 20 0 0 600 900 6000]')'

pop_a = vec(pop')'
pop_all = pop_a
for uk in range(1,size(pop)[1]^2-1)
  pop_all = vcat(pop_all,pop_a)
end

for jkl in range(1,2*size(pop)[1]^2)
  B = hcat(B,0)
end

Lst = copy(L)
for jkl in range(1,2*size(pop)[1]^2)
  Lst = hcat(Lst,0)
end

D = B.*(1.-Lst)

#----------------------------------------#
# Roads

roadmax = size(road)[1]

roads = eye(size(pop)[1]^2)
ports = eye(size(pop)[1]^2)
aeroplanes = eye(size(pop)[1]^2)

roadstart = 1
while roadstart <= roadmax
  place = road[roadstart,1]
  pdg = digits(place)
  pdg2 = pdg[end:-1:1,end:-1:1]
  down = (pdg2[1]-1)*7+pdg2[2]

  other = road[roadstart,2]
  digs = digits(other)[end:-1:1,end:-1:1]
  len = size(digs)[1]
  for xyz in range(1,convert(Int64,len/2))
    down2 = (digs[(xyz*2)-1]-1)*size(pop)[1]+digs[(xyz*2)]
    roads[down,down2]=1
  end
  roadstart += 1
end
currentroad = (1-eye(size(pop)[1]^2)).*roads.*pop_all
totalroad = sum(currentroad,2)
roadx = 1.6 /(1.6+0.4+0.1)
ratioroad = currentroad ./ totalroad

roadcorrect = copy(ratioroad)
for m=1:size(roadcorrect,1)
	for l=1:size(roadcorrect,2)
  		 isnan(roadcorrect[m,l]) ? roadcorrect[m,l]=0 :	 roadcorrect[m,l]=roadcorrect[m,l]
	end
end

#----------------------------------------#
# Aeroplane

aeromax = size(aeroplane)[1]
aerostart = 1
while aerostart <= aeromax
  place = aeroplane[aerostart,1]
  pdg = digits(place)
  pdg2 = pdg[end:-1:1,end:-1:1]
  down = (pdg2[1]-1)*7+pdg2[2]

  other = aeroplane[aerostart,2]
  digs = digits(other)[end:-1:1,end:-1:1]
  len = size(digs)[1]
  for xyz in range(1,convert(Int64,len/2))
    down2 = (digs[(xyz*2)-1]-1)*size(pop)[1]+digs[(xyz*2)]
    aeroplanes[down,down2]=1
  end
  aerostart += 1
end
currentaero = (1-eye(size(pop)[1]^2)).*aeroplanes.*pop_all
totalaero = sum(currentaero,2)
aerox = 0.1 /(1.6+0.4+0.1)
ratioaero = currentaero ./ totalaero

aerocorrect = copy(ratioaero)
for m=1:size(aerocorrect,1)
	for l=1:size(aerocorrect,2)
  		 isnan(aerocorrect[m,l]) ? aerocorrect[m,l]=0 :	 aerocorrect[m,l]=aerocorrect[m,l]
	end
end

#----------------------------------------#
# Ports

portmax = size(port)[1]

portstart = 1
while portstart <= portmax
  place = port[portstart,1]
  pdg = digits(place)
  pdg2 = pdg[end:-1:1,end:-1:1]
  down = (pdg2[1]-1)*7+pdg2[2]

  other = port[portstart,2]
  digs = digits(other)[end:-1:1,end:-1:1]
  len = size(digs)[1]
  for xyz in range(1,convert(Int64,len/2))
    down2 = (digs[(xyz*2)-1]-1)*size(pop)[1]+digs[(xyz*2)]
    ports[down,down2]=1
  end
  portstart += 1
end
currentport = (1-eye(size(pop)[1]^2)).*ports.*pop_all
totalport = sum(currentport,2)
portx = 0.4 /(1.6+0.4+0.1)
ratioport = currentport ./ totalport

portcorrect = copy(ratioport)
for m=1:size(portcorrect,1)
	for l=1:size(portcorrect,2)
  		 isnan(portcorrect[m,l]) ? portcorrect[m,l]=0 :	 portcorrect[m,l]=portcorrect[m,l]
	end
end

#----------------------------------------#
## ADDITIONS
rr = roadcorrect.*roadx.*0.015
pp = portcorrect.*portx.*0.015
aa = aerocorrect.*aerox.*0.015
needleft = 1 .- sum(rr+pp+aa,2)

eyeget = eye(size(pop)[1]^2) .* needleft

ssss = eyeget + aa + rr + pp

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

large = size(ssss)[1]

Ipart = zeros(1,size(pop)[1]^2)
Ipart[1,13]=1
Rpart = zeros(1,size(pop)[1]^2)

SIR = hcat(pop_a,Ipart,Rpart)

table = ["s1"]
for x in range(2,size(pop)[1]^2-1)
  h = string("s",x)
  table = hcat(table,h)
end
for x in range(1,size(pop)[1]^2)
  h = string("i",x)
  table = hcat(table,h)
end
for x in range(1,size(pop)[1]^2)
  h = string("r",x)
  table = hcat(table,h)
end

old_SIR = zeros(1,large*3)

mosquito = zeros(1,large)
old_ir = zeros(large,1)

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

for x in range(1,iterations)

  temp = rand(10:44,1,large)
  prec = rand(100:400,1,large)

#----------------------------------------#
  # Temperature first
  sd = 7;  mean = 27;  sc = 13.5
  h = (1/(sd*sqrt(2*pi)))*sc
  f(t) = (exp((-(t-mean).^2)/(2*sd^2)))*h

  T=f(temp)
  total = sum(T)
  ss = size(T)[2]
  temperature = eye(ss)

  for uu in 1:ss
    Q = copy(vec(T))
    Q[uu] = 0
    jj = 0

    for x in Q
      if x == 0
        d = T[uu]
      else
        y = total-T[uu]
        b = x/y
        c = 1-T[uu]
        d = c*b
      end
      jj += 1

      temperature[uu,jj] = d
    end
  end
#----------------------------------------#

  # Precipitation next
  sd = 50;  mean = 250;  sc = 100
  h = (1/(sd*sqrt(2*pi)))*sc
  f(t) = (exp((-(t-mean).^2)/(2*sd^2)))*h

  T=f(prec)
  total = sum(T)
  ss = size(T)[2]
  precipitation = eye(ss)

  for uu in 1:ss
    Q = copy(vec(T))
    Q[uu] = 0
    jj = 0

    for x in Q
      if x == 0
        d = T[uu]
      else
        y = total-T[uu]
        b = x/y
        c = 1-T[uu]
        d = c*b
      end
      jj += 1

      precipitation[uu,jj] = d
    end
  end
#----------------------------------------#

  # Living standards next
  T = 1 - L
  total = sum(T)
  ss = size(T)[2]
  living = eye(ss)

  for uu in 1:ss
    Q = copy(vec(T))
    Q[uu] = 0
    jj = 0

    for x in Q
      if x == 0
        d = T[uu]
      else
        y = total-T[uu]
        b = x/y
        c = 1-T[uu]
        d = c*b
      end
      jj += 1

      living[uu,jj] = d
    end
  end
#----------------------------------------#

  # Now probability for mosquito

  mosq_p = temperature*precipitation*living

  mosq_b = rand(1000:2000,1,large)
  mosq_d = rand(500:600,1,large)
  w_scalar = 0.7

  mosquito = (mosquito + (mosq_b*0.7 - mosq_d))*mosq_p

  # Now scaling for mosquito

  addup = sum(mosquito)
  elements = size(mosquito)[2]
  average = addup / elements

  addon_all = 0
  for ippa in mosquito
    addon = (ippa - average)^2
    addon_all = addon + addon_all
  end

  sd = (1/elements*addon_all)^0.5
  z = (mosquito .- average) ./ sd

  a = 111.*atan(27.*z/294)
  b = -358.*z./23
  c = b+a
  d = exp(c)
  attack_rate = 1./(1.+d)
#----------------------------------------#

  # S -> I matrix

  u = SIR[large+1:2*large] ./ (SIR[1:large]+SIR[large+1:2*large]+SIR[2*large+1:3*large])
  # pq = (SIR[large+1:2*large] .- old_SIR[1+large:2*large])
  # qrd = (SIR[1:large] .- old_SIR[1:large])

  ratio = 0.5
  IR = attack_rate' .* u .* ratio
  # constant = 0.1
  # IRRR = (2*exp(1/2*constant*x))/(1+exp(1/2*constant*x))
  # IRR = IRRR .* (SIR[10:18]./(SIR[1:9].+SIR[10:18].+SIR[19:27]))
  # IR = IRRR .* IRR
  kk = 1
  infection = zeros(elements)
  for xp in isnan(IR)*1
    if xp == 0
      infection[kk] = IR[kk]
    else
      infection[kk] = 0.0
    end
    kk += 1
  end

  old_SIR = copy(SIR)
  old_ir = copy(IR)

  s_i = infection' .* eye(elements)
  global(s_i)

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

  # I -> R matrix
  opp = 20
  i_r1 = L ./opp
  i_r = eye(elements).*i_r1

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

  # S -> S matrix
  nrow = size(ssss)[1]
  h = sum(eye(nrow,nrow)-s_i,2)

  s_s = ssss.*h

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

  # I -> I matrix
  leftover = (1-sum(i_r,2))
  scaled = leftover.*s_s
  other = (scaled.*(1-eye(large))).*0.5
  othertotal = sum(other,2)
  correct = (scaled.+othertotal).*eye(large).+other
  i_i = (leftover-sum(correct,2).+correct).*eye(large).+correct.*(1-eye(large))

#----------------------------------------#
  # R -> R matrix

  gg = (leftover.-s_s).*eye(nrow,nrow).*2/nrow
  hh = 1 .- sum(s_s+gg,2)
  qp = (1.-eye(nrow,nrow)).*s_s
  vv = qp ./ sum(qp,2) .*hh

  ooo = 1
  for xpd in isnan(vv)*1
    if xpd == 1
      vv[ooo] = 0
    end
    ooo += 1
  end
  r_r = gg+s_s+vv

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

  z = zeros(nrow,nrow)
  P1 = hcat(s_s,s_i,z)
  P2 = hcat(z,i_i,i_r)
  P3 = hcat(z,z,r_r)

  p = vcat(P1,P2,P3)

  SIR = (SIR+(B*beta-D))*p
  table = vcat(table,SIR)

end

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

table = table'
starting = ["Date"]
for x in range(1,iterations)
  starting = hcat(starting,x)
end

table = vcat(starting,table)
writecsv("Data.csv",table, header = false)

max_stars_repo_path: 49 square.jl

max_issues_repo_path: 49 square.jl

max_forks_repo_path: 49 square.jl


text:
using CircuitscapeUI

function change_dir_if_bundle()
    cd(joinpath(split(string(Base.julia_cmd()))[1][2:end], ".."))
end

# Blink
Core.eval(CircuitscapeUI.Blink.AtomShell, :(_electron = "../Libraries/Julia.app/Contents/MacOS/Julia"))
Core.eval(CircuitscapeUI.Blink.AtomShell, :(mainjs = "main.js"))
Core.eval(CircuitscapeUI.Blink, :(buzz = "main.html"))
Core.eval(CircuitscapeUI.Blink, :(resources = Dict("spinner.css" => "res/spinner.css",
                               "blink.js" => "res/blink.js",
                               "blink.css" => "res/blink.css",
                               "reset.css" => "res/reset.css")))
Core.eval(CircuitscapeUI.Blink, :(webiosetup = "webio_setup.js"))

# Mbedtls
Core.eval(CircuitscapeUI.MbedTLS, :(libmbedcrypto = "../Libraries/libmbedcrypto"))
Core.eval(CircuitscapeUI.MbedTLS, :(libmbedtls = "../Libraries/libmbedtls"))
Core.eval(CircuitscapeUI.MbedTLS, :(libmbedx509 = "../Libraries/libmbedx509"))

# Tachyons
Core.eval(CircuitscapeUI.Tachyons, :(path = "tachyons.min.css"))

# WebIO
Core.eval(CircuitscapeUI.WebIO, :(bundlepath = "bundle.js"))

# InteractBase
Core.eval(CircuitscapeUI.InteractBase, :(font_awesome = "all.js"))
Core.eval(CircuitscapeUI.InteractBase, :(style_css = "style.css"))

# InteractBulma
Core.eval(CircuitscapeUI.InteractBulma, :(font_awesome = "all.js"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_min_css = "bulma.min.css"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_slider_min_css = "bulma-slider.min.css"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_switch_min_css = "bulma-switch.min.css"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_checkradio_min_css = "bulma-checkradio.min.css"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_tooltip_min_css = "bulma-tooltip.min.css"))
Core.eval(CircuitscapeUI.InteractBulma, :(bulma_accordion_min_css = "bulma-accordion.min.css"))

# Knockout
Core.eval(CircuitscapeUI.Knockout, :(knockout_js = "knockout.js"))
Core.eval(CircuitscapeUI.Knockout, :(knockout_punches_js = "knockout_punches.js"))

# CircuitscapeUI
Core.eval(CircuitscapeUI, :(TESTPATH = "test/"))
Core.eval(CircuitscapeUI, :(logo = "cs_logo.ico"))

Base.@ccallable function julia_main(args::Vector{String})::Cint
    # ApplicationBuilder.App.change_dir_if_bundle()
    # change_dir_if_bundle()
    @show pwd()
    @show CircuitscapeUI.Blink.AtomShell._electron
    @show CircuitscapeUI.WebIO.bundlepath
    @show CircuitscapeUI.InteractBulma.bulma_tooltip_min_css
    @show CircuitscapeUI.InteractBulma.bulma_accordion_min_css
    w = CircuitscapeUI.run_ui()
    @show CircuitscapeUI.AssetRegistry.registry
    while CircuitscapeUI.Blink.active(w)
        sleep(0.01)
    end
    return 0
end


text:
abstract type GeometricInterpretor end

struct GNNInterpretation{Interpretor <: GeometricInterpretor}
    model :: Chain
    fg :: FeaturedGraph
    X :: Matrix{Real}
    y :: Array{Int}
    interpretor :: Interpretor
end

struct GNNInterpretationResults
    edge_mask :: FeaturedGraph
    node_feat_mask :: Union{Matrix{Real}, Nothing}
end

function explain_node(method :: GNNInterpretation{GeometricInterpretor},
    node :: Int) :: GNNInterpretationResults end

function visualize_explanation(edge_mask :: FeaturedGraph)
    swg = graph(edge_mask)
    graphplot(swg |> SimpleDiGraph, arrow_shift = 0.75, arrow_size=20, edge_width=2.5, edge_color=:darkgray, node_color=[3, 5, 3]) # add weights, node colors, node labels
end


text:
using DataGeneratorTranslators
using Base.Test

@testset "type" begin

	@testset "generators" begin
		include("test_generators.jl")
	end

end


text:
module QPWalkingControl

export
    PushRecoveryController,
    ICPController,
    ICPTrajectoryGenerator,
    ICPWalkingStateMachine,
    PushApplier,
    Widget, # from InteractBase
    PushRecoveryVisualizer,
    add_sole_frames!,
    make_foot_polygons,
    icp,
    horizontal_projection,
    transfer_weight!

# consider moving to RBD.jl
export
    in_frame,
    @framechecked,
    unwrap

# consider moving to RigidBodySim.jl
export
    JointDamping

using LinearAlgebra
using RigidBodyDynamics
using RigidBodyDynamics.Contact
using RigidBodyDynamics.PDControl
using StaticArrays
using Rotations
using QPControl
using QPControl.Trajectories
using Parametron
using PlanarConvexHulls
using StaticUnivariatePolynomials
import MathOptInterface

using Base.Iterators: cycle
using RigidBodyDynamics.Graphs: TreePath, target
using RigidBodyDynamics: frame_definition
using StaticUnivariatePolynomials: derivative
using QPControl.Trajectories: exponential_integral, breaks, subfunctions

const SUP = StaticUnivariatePolynomials
const MOI = MathOptInterface
const RBD = RigidBodyDynamics

const SPoint3D{T} = Point3D{SVector{3, T}}

include("frames.jl")
include("util.jl")

include("setup.jl")
include("planarcontactmode.jl")
include("controller.jl")

include("damping.jl")
include("pushapplier.jl")

include("foot_trajectories/basic_foot_trajectory.jl")
include("foot_trajectories/foot_trajectory_gen.jl")

include("icp/icp.jl")
include("icp/control.jl")
include("icp/trajectory_gen.jl")
include("icp/trajectory_init.jl")
include("icp/statemachine.jl")

include("visualization.jl")

using .Visualization

end # module


text:
module Points

export Point, get_x, get_y, get_z, get_coords, unit_vectors2, unit_vectors3
export UnitVector, Biased, Unbiased

###
### Point
###

struct Point{N, T}
    coords::NTuple{N, T}
end

get_coords(p::Point) = p.coords
Point(coords::T...) where T = Point{length(coords), T}(coords)
Point(coords::Tuple) = Point(coords...)

@generated function Point{N, T}() where {N, T}
    isa(N, Integer) || throw(TypeError(:Point, "", Integer, typeof(N)))
    N > 0 || throw(DomainError(N, "N must be a positive Integer"))
    return Point((zero(T) for i in 1:N)...,)
end

Point{N}() where {N} = Point{N, Float64}()
Point() = Point{1}()

Base.getindex(p::Point, inds...) = getindex(get_coords(p), inds...)
Base.length(p::Point) = length(get_coords(p))
Base.eltype(p::Point) = eltype(get_coords(p))

get_x(p::Point) = p[1]
get_y(p::Point) = p[2]
get_z(p::Point) = p[3]

import Base: +, -, *, ==
+(p::Point{1, <:Any}, x::Number) = Point(x + get_x(p))
+(p1::Point{N}, p2::Point{N}) where N = Point((get_coords(p1) .+ get_coords(p2))...)
*(p1::Point, n::Number) = Point((get_coords(p1) .* n))
*(n::Number, p::Point) = p * n
-(p::Point) = Point(broadcast(-, get_coords(p))...)
-(p1::Point{N}, p2::Point{N}) where N = p1 + (-p2)
==(p1::Point{N}, p2::Point{N}) where N = get_coords(p1) == get_coords(p2)
norm_squared(p::Point) = sum(x -> x^2, p.coords)

Base.zero(p::Point) = Base.zero(typeof(p))
Base.zero(::Type{Point{N, T}}) where {T, N} = Point{N, T}()
Base.iszero(p::Point{N, T}) where {N, T} = p == Point{N, T}()
Base.show(io::IO, p::Point{N, T}) where {T, N}  = print(io, "Point{$N, $T}", get_coords(p))
Base.show(io::IO, p::Point{1, T}) where {T}  = print(io, "Point{1, $T}(", get_coords(p)[1], ")")

const unit_vectors1 = (Point(1), Point(-1))

const unit_vectors2 = (Point(0, 1), Point(1, 0), Point(0, -1), Point(-1, 0))

const unit_vectors3 = (Point(0, 0, 1), Point(0, 1, 0), Point(1, 0, 0),
                       Point(0, 0, -1), Point(0, -1, 0), Point(-1, 0, 0))

struct UnitVector{P} end

struct Biased
    bias::Float64
    cutoff::Float64
end
Biased(bias=0.0) = Biased(bias, 1/2 + bias)

struct Unbiased end

Base.rand(ub::Unbiased, t::Type{UnitVector{Point{N, T}}}) where {N, T} = rand(t)
Base.rand(::Type{UnitVector{Point{N, T}}}) where {N, T} = rand(UnitVector{Point{N}})
Base.rand(::Type{UnitVector{Point{1}}}) = rand(Bool) ? Point(1) : Point(-1)
Base.rand(::Type{UnitVector{Point{2}}}) = unit_vectors2[rand(1:4)]
Base.rand(::Type{UnitVector{Point{3}}}) = unit_vectors3[rand(1:6)]
Base.rand(b::Biased, ::Type{UnitVector{Point{N, T}}}) where {N, T} = rand(b, UnitVector{Point{N}})
Base.rand(b::Biased, ::Type{UnitVector{Point{1}}}) = rand() > b.cutoff ? Point(1) : Point(-1)

end # module Point


text:
"""
Example using NLPJuMP with Ipopt
"""

using JuMP
push!(LOAD_PATH, "../src/")
using NLPSaUT

NLPSaUT.test_ipopt()

# """
# User-defined fitness function

# # Returns:
# 	- `Vector`: vector stores [f, h, g] where f is objective, h is equality, and g is inequality constraints
# """
# function f_fitness(x...)
# 	# objective
#     f = x[1]^2 - x[2]

#     # equality constraints
#     h = zeros(1,)
#     h = x[1]^3 + x[2] - 1

#     # inequality constraints
#     g = zeros(2,)
#     g[1] = x[1]^2 + x[2]^2 - 1
#     g[2] = x[2] - 2

#     fitness = vcat(f,h,g)[:]
#     return fitness
# end

# # ------------------------------------------------------------------------- #
# ## Solve model
# # problem dimensions
# nx = 2                   # number of decision vectors
# nh = 1                   # number of equality constraints
# ng = 2                   # number of inequality constraints
# lx = -10*ones(nx,)
# ux =  10*ones(nx,)
# x0 = [1.2, 0.9]

# # get model
# order = 2
# diff_f = "forward"
# x, model = build_model(f_fitness, nx, nh, ng, lx, ux, x0)#, order, diff_f)
# set_optimizer_attribute(model, "tol", 1e-6)
# set_optimizer_attribute(model, "print_level", 5)
# println(model)

# # run optimizer
# optimize!(model)

# # print results
# println(termination_status(model))
# println("Decision vector: ")
# println(value.(x))
# println("Objective: ")
# println(objective_value(model))

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
@defcomp LGforcing begin

    c_LGconcentration=Parameter(index=[time],unit="ppbv")
    c0_LGconcbaseyr=Parameter(unit="ppbv", default=0.11)
    f_LGforcing=Variable(index=[time],unit="W/m2")
    f0_LGforcingbase=Parameter(unit="W/m2", default=0.022)
    fslope_LGforcingslope=Parameter(unit="W/m2", default=0.2)

    function run_timestep(p, v, d, t)

        #eq.13 in Hope 2006
        v.f_LGforcing[t]=p.f0_LGforcingbase+p.fslope_LGforcingslope*(p.c_LGconcentration[t]-p.c0_LGconcbaseyr)

    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
# Optimization and solution query

function optimize(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(optimize, Cint, (Ptr{Cvoid},), model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function computeIIS(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(computeIIS, Cint, (Ptr{Cvoid},), model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

#################################################
#
#  solution status and optimization info
#
#################################################

const GRB_LOADED          = 1
const GRB_OPTIMAL         = 2
const GRB_INFEASIBLE      = 3
const GRB_INF_OR_UNBD     = 4
const GRB_BOUNDED         = 5
const GRB_CUTOFF          = 6
const GRB_ITERATION_LIMIT = 7
const GRB_NODE_LIMIT      = 8
const GRB_TIME_LIMIT      = 9
const GRB_SOLUTION_LIMIT  = 10
const GRB_INTERRUPTED     = 11
const GRB_NUMERIC         = 12
const GRB_SUBOPTIMAL      = 13
const GRB_INPROGRESS      = 14
const GRB_USER_OBJ_LIMIT  = 15

const status_symbols = [
    :loaded,
    :optimal,
    :infeasible,
    :inf_or_unbd,
    :unbounded,
    :cutoff,
    :iteration_limit,
    :node_limit,
    :time_limit,
    :solution_limit,
    :interrupted,
    :numeric,
    :suboptimal,
    :inprogress,
    :user_obj_limit
]

get_status_code(model::Model) = get_intattr(model, "Status")
get_status(model::Model) = status_symbols[get_status_code(model)]::Symbol

@grb_dbl_attr get_runtime      "Runtime"
@grb_dbl_attr get_objval       "ObjVal"
@grb_dbl_attr get_objbound     "ObjBound"

@grb_int_attr get_sol_count    "SolCount"
@grb_int_attr get_barrier_iter_count "BarIterCount"
#@grb_dbl_attr get_node_count   "NodeCount"

get_iter_count(model::Model) = convert(Int, get_dblattr(model, "IterCount"))
get_node_count(model::Model) = convert(Int, get_dblattr(model, "NodeCount"))

mutable struct OptimInfo
    status::Symbol
    runtime::Float64

    sol_count::Int
    iter_count::Int
    barrier_iter_count::Int
    node_count::Int
end

function get_optiminfo(model::Model)
    OptimInfo(
        get_status(model),
        get_runtime(model),

        get_sol_count(model),
        get_iter_count(model),
        get_barrier_iter_count(model),
        get_node_count(model)
    )
end

function show(io::IO, s::OptimInfo)
    println(io, "Gurobi Optimization Info")
    println(io, "    status   = $(s.status)")
    println(io, "    runtime  = $(s.runtime)")
    println(io, "    # solutions = $(s.sol_count)")
    println(io, "    # iters     = $(s.iter_count)")
    println(io, "    # bar iters = $(s.barrier_iter_count)")
    println(io, "    # nodes     = $(s.node_count)")
end

#################################################
#
#  solution query
#
#################################################

get_solution(model::Model) = get_dblattrarray(model, "X", 1, num_vars(model))

const varmap = Dict(
    -3 => :Superbasic,
    -2 => :NonbasicAtUpper,
    -1 => :NonbasicAtLower,
     0 => :Basic
)

const conmap = Dict(
     0 => :Basic,
    -1 => :Nonbasic)

const basicmap_rev = Dict(
    :Superbasic => Cint(-3),
    :NonbasicAtUpper => Cint(-2),
    :NonbasicAtLower => Cint(-1),
    :Basic => Cint(0),
    :Nonbasic => Cint(-1)
)

function get_basis(model::Model)
    cval = Array{Cint}(undef, num_vars(model))
    cbasis = Array{Symbol}(undef, num_vars(model))
    get_intattrarray!(cval, model, "VBasis", 1)
    for it in 1:length(cval)
        cbasis[it] = varmap[cval[it]]
    end

    rval = Array{Cint}(undef, num_constrs(model))
    rbasis = Array{Symbol}(undef, num_constrs(model))
    get_intattrarray!(rval, model, "CBasis", 1)
    rsense = Array{Cchar}(undef, num_constrs(model))
    get_charattrarray!(rsense, model, "Sense", 1)
    for it in 1:length(rval)
        rbasis[it] = conmap[rval[it]]
        if rbasis[it] == :Nonbasic
            if rsense[it] == convert(Cchar,'<')
                rbasis[it] = :NonbasicAtUpper
            else
                rbasis[it] = :NonbasicAtLower
            end
        end
    end
    return cbasis, rbasis
end

"""
    loadbasis(model::Model, x::Vector)

Load basis to a problem in form of primal solution.

    loadbasis(model::Model, rval::Vector{Symbol}, cval::Vector{Symbol})

Load basis to a problem in terms of basicness description Variables columns (in
`cval`) can be: `:Basic`, `:NonbasicAtLower`, `:NonbasicAtUpper` and
`:Superbasic` Constraints rows (in `rval`) can be: `:Basic`, `:Nonbasic`
"""
function loadbasis(model::Model, x::Vector)

    ncols = num_vars(model)
    nrows = num_constrs(model)

    length(x) != ncols && error("solution candidate size is different from the number of columns")

    cvals = Array{Cint}(undef, ncols)
    rvals = Array{Cint}(undef, nrows)

    # obtain situation of columns

    lb = get_dblattrarray(model, "LB", 1, ncols)
    ub = get_dblattrarray(model, "UB", 1, ncols)

    for i in 1:ncols
        if isapprox(x[i],lb[i])
            cvals[i] = basicmap_rev[:NonbasicAtLower]
        elseif isapprox(x[i],ub[i])
            cvals[i] = basicmap_rev[:NonbasicAtUpper]
        else
            cvals[i] = basicmap_rev[:Basic]
        end
    end

    # obtain situation of rows where: y = Ax

    A = get_constrmatrix(model) #A is sparse

    y = A*x

    senses = get_charattrarray(model, "Sense", 1, nrows)
    rhs    = get_dblattrarray(model, "RHS", 1, nrows)

    for j in 1:nrows
        if senses[j] == '=' && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:Nonbasic]#AtLower
        elseif senses[j] == '>' && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:Nonbasic]#AtLower
        elseif senses[j] == '<' && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:Nonbasic]#AtUpper
        else
            rvals[j] = basicmap_rev[:Basic]
        end
    end

    loadbasis(model, rvals, cvals)

    return nothing
end
function loadbasis(model::Model, rval::Vector{Symbol}, cval::Vector{Symbol})

    nrval = map(x->conmap[x], rval)
    ncval = map(x->varmap[x], cval)

    loadbasis(model, nrval, ncval)

    return nothing
end
function loadbasis(model::Model, rval::Vector{Cint}, cval::Vector{Cint})

    ncols = num_vars(model)
    nrows = num_constrs(model)

    set_intattrarray!(model, "VBasis", 1, num_vars(model), cval)
    set_intattrarray!(model, "CBasis", 1, num_constrs(model), rval) # r = row; c = constr

    return nothing
end


text:
obs(env::PendulumEnv, ::Nothing) = obs(env, [0.0, 0.0])
obs(env::PendulumEnv, (θ, θ̄)) = Flux.data(θ)

function Ctx(env::PendulumEnv)
    path = (p) -> normpath("$(@__DIR__)/$p")
    s = Scope(imports=path.(["../../assets/utils.js", "../../assets/pendulum/pendulum.js", "../../assets/pendulum/pendulum.css"]))
    o = Observable(s, "obs", obs(env, env.state))
    onimport(s, @js () -> begin
        window.pick = (e) -> document.querySelector(e)
        window.container = window.pick(".wio-scope") || window.pick(".webio-scope")
        window.p = __init__(container, $o[]);
        p.draw();
    end)

    onjs(o, @js function (x)
        p.set_theta(x)
        p.draw()
    end)

    Ctx(s, o)
end

function render(env::PendulumEnv, ctx::Ctx)
    ctx.o[] = obs(env, env.state)
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to compute statistics for estimating sectoral trade elasticities (method of Caliendo and Parro, 2015)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    CP_statistic(X::Matrix, τ_X::Matrix, N::Integer, S::Integer)

The function computes the LHS and RHS statistics of equation (23) in Caliendo and Parro (2015).

# Arguments
- `X::Matrix`: NS×N, origin country-industry destination country intermediate/final demand matrix X.
- `τ_X::Matrix`: NS×N, origin country-industry destination country intermediate/final demand tariff matrix τ_X.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `lhs_temp::Vectpr{Float64}`: (Σᴺ⁻²n(n+1)/2)×1, LHS of eq. (23) in logs for intermediate/final demand per industry.
- `rhs_temp::Vectpr{Float64}`: (Σᴺ⁻²n(n+1)/2)×1, RHS of eq. (23) in logs for intermediate/final demand per industry.

# Examples
```julia-repl
julia> lhs_Z_temp, rhs_Z_temp, lhs_F_temp, rhs_F_temp = create_reg_data(X_temp, F_temp, τ_X_temp, τ_F_temp)
```
"""

function CP_statistic(X::Matrix, τ_X::Matrix, N::Integer, S::Integer)

    # initialize loop
    LHS = Float64[]
    RHS = Float64[]

    for i in 1:N-2 # i -> j
        for j in i+1:N-1 # j -> n
            for n in j+1:N # n -> i
                push!( LHS, ( X[i,j]*X[j,n]*X[n,i] ) / ( X[j,i]*X[n,j]*X[i,n] ) )
                push!( RHS, ( τ_X[i,j]*τ_X[j,n]*τ_X[n,i] ) / ( τ_X[j,i]*τ_X[n,j]*τ_X[i,n] ) )
            end
        end
    end

    # create series in logs
    lhs_temp = log.(LHS)
    rhs_temp = log.(RHS)

    return lhs_temp, rhs_temp
end

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

"""
    elasticity_data(M::Matrix, τ::Matrix, demand::String, N::Integer, S::Integer)

A summary function to collect all industry LHS and RHS statistics of equation (23) in Caliendo and Parro (2015).
    Notice that intermediate and final demand are treated seperately.

# Arguments
- `M::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `τ::Matrix`: NS×NS, origin country-industry destination country-industry intermediate input tariff matrix τ_Z.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `lhs::Vectpr{Float64}`: S*(Σᴺ⁻²n(n+1)/2)×1, LHS of eq. (23) in logs for intermediate inputs and all industries.
- `rhs::Vectpr{Float64}`: S*(Σᴺ⁻²n(n+1)/2)×1, RHS of eq. (23) in logs for intermediate inputs and all industries.

# Examples
```julia-repl
julia> lhs, rhs = elasticity_data(Z, F, τ_Z, τ_F, N, S)
```
"""

function elasticity_data(M::Matrix, τ::Matrix, demand::String, N::Integer, S::Integer)

    if demand == "intermediate"
        # reduce intermediate demand to origin country-industry destination country
        X = [sum(M[i,j:j+S-1]) for i in 1:N*S, j in 1:S:N*S] # NS×N

    elseif demand == "total" # wrong
        X = [sum(M[i,j:j+S-1]) for i in 1:N*S, j in 1:S:N*S] # NS×N
        X = X .+ F
    elseif demand == "final" # wrong, dont understand how to reshape final demand matrix
        X = M
    else
        println(" × Failure! Specify demand schedule properly! \n")
    end

    τ_X = τ

    # initialize loop
    lhs = Float64[]
    rhs = Float64[]

    for j in 1:S

        # same sector in rows
        X_temp = X[j:S:(N-1)*S+j, :] # N×N
        τ_X_temp = τ_X[j:S:(N-1)*S+j, :] # N×N

        lhs_temp, rhs_temp = CP_statistic(X_temp, τ_X_temp, N, S)

        lhs = [lhs; lhs_temp]
        rhs = [rhs; rhs_temp]

    end

    # zero trade causes inf => create NaN instead which we let regression filter out later
    lhs = ifelse.(isinf.(lhs), NaN, lhs)
    rhs = ifelse.(isinf.(rhs), NaN, rhs)

    # check if computed length equals number of (theoretical) combinations
    combinations = Int64(sum([n*(n+1)/2 for n in 1:N-2])) # number of possible combinations from formula
    check = ([length(lhs) length(rhs)] .== S*combinations) == [1 1]
    if check == true
        println(" ✓ The length of the vector coincides with the theoretical length of: $combinations \n")
    else
        println(" × Failure! The length of the vector does not! coincide with the theoretical length of: $combinations \n")
    end

    return lhs, rhs

end


text:
export BayesNetDAI, BayesNetSampler

type BayesNetDAI
    verts::Vector{Var}
    fg::FactorGraph
    ia::InfAlg
    dirty::Bool
    memo_entropy::Float64
end

length(net::BayesNetDAI) = length(net.verts)

function BayesNetDAI(n::Int)
    verts = [Var(x,2) for x=1:n] # FIXME hardcoded binary
    fg = FactorGraph([Factor(x) for x in verts]...)
    ia = InfAlg(fg)
    BayesNetDAI(verts, fg, ia, true, -Inf)
end

function random_net_tree(n::Int=5, p=0.1)
    net = BayesNetDAI(n)
    if n==1
        move_params!(net,1)
        return net
    end

    for i=2:n
      add_edge!(net, i, rand(1:i-1))
    end

    for i=1:n
      move_params!(net, i, false)
    end
    net
end

function random_net_cascade(n::Int=5, p=0.1)
    net = BayesNetDAI(n)
    if n==1
        move_params!(net,1)
        return net
    end
    hubs = 1 #rand(1:3)
    n3 = fld(n,3)
    tert1 = 1:n3
    tert2 = n3+1:n3*2
    tert3 = n3*2+1:n

    for i=hubs+1:n3*2
      add_edge!(net, rand(1:i-1), i, false)
    end

    for i=n3*2+1:n
      if i<=n3*2+2 || rand() < p
        add_edge!(net, i, rand(1:hubs), false)
      else
        add_edge!(net, i, rand(n3*2+1:i-1), false)
      end
    end

    for i=1:n
      move_params!(net,i, false)
    end
    net
end

function random_net_erdos(n::Int=5, p=0.2)
        net = BayesNetDAI(n)
        if n==1
            move_params!(net,1)
            return net
        end
        for i=1:n
            for j=1:n
              if i != j && rand() < p
                add_edge!(net, j, i, false)
              end
            end
            move_params!(net,i, false)
        end
    net
end

function random_net_weights(n::Int=5, weights=WeightVec([0.5,0.3,0.1]))
    count = -1
    cyclic = true
    local net
    while cyclic
        net = BayesNetDAI(n)
        if n==1
            move_params!(net,1)
            return net
        end
        for i=1:n
            for j=1:(sample(weights)-1)
                out = i
                while out == i
                    out = rand(1:n)
                end
                add_edge!(net, out, i)
            end
            move_params!(net,i)
        end
        g = net2graph(net)
        cyclic = test_cyclic_by_dfs(g)
        count += 1
    end
    net
end

function net2graph(net)
    N = length(net.verts)
    g = simple_inclist(N)
    for i=1:N
        fac = net.fg[i]
        allvars = vars(fac)
        pars = allvars - net.verts[i]
        for n in labels(pars)
            add_edge!(g, n, i)
        end
    end
    g
end

function net2mat(net)
    g = net2graph(net)
    adjacency_matrix(g)
end

function draw_data(net, num)
    g = net2graph(net)
    tsort = topological_sort_by_dfs(g)
    data = zeros(Int, num_vertices(g), num)
    for i=1:num, j in tsort
        fac = net.fg[j]
        allvars = vars(fac)
        pars = allvars - net.verts[j]
        if length(pars) > 0
            pstates = Int[]
            for k in labels(pars)
                push!(pstates, data[findfirst(tsort,k),i])
            end
            pstate = calcLinearState(pars, pstates)
            i1,i2 = conditionalStateBoth(net.verts[j], pars, pstate)
            w = WeightVec([fac[i1], fac[i2]])
            data[j,i] = sample(1:2, w)
        else
            data[j,i] = sample(1:2, WeightVec(p(fac)))
        end
    end
    data
end

function clearBackups!(bnd::BayesNetDAI)
  clearBackups!(bnd.fg)
  bnd.dirty = true
end

function add_edge!(bnd::BayesNetDAI, u::Int, v::Int, backup=true)
    un,vn = bnd.verts[[u,v]]
    vnfac = bnd.fg[v]
    #assert(!(un in vars(vnfac)))
    un in vars(vnfac) && return
    bnd.dirty=true

    if backup
      setBackedFactor!(bnd.fg, v, embed(vnfac, un))
    else
      bnd.fg[v] = embed(vnfac, un) #Factor(vars(vnfac)+un)
    end
end

function del_edge!(bnd::BayesNetDAI, u::Int, v::Int, backup=true)
    un,vn = bnd.verts[[u,v]]
    bnd.dirty=true

    fac = bnd.fg[v]
    allvars = vars(fac)
    newvars = allvars - un
    newpars = newvars - vn
    newfac = Factor(newvars)

    for parstate in 1:nrStates(newpars)
        # FIXME: hardcoded binary
        i1,i2 = conditionalState2(allvars, vn, un, newpars, 1, parstate)
        newval = (fac[i1] + fac[i2])/2
        c1,c2 = conditionalStateBoth(vn, newpars, parstate)
        newfac[c1] = newval
        newfac[c2] = 1 - newval
    end
    if backup
      setBackedFactor!(bnd.fg, v, newfac)
    else
      bnd.fg[v] = newfac
    end
end

function set_factor!(bnd::BayesNetDAI, i::Int, vals::Vector{Float64}, backup=true)
    bnd.dirty = true
    fcopy = bnd.fg[i]
    assert(length(vals) == nrStates(fcopy))
    oldvals = p(fcopy)
    oldvals[:] = vals
    if backup
      setBackedFactor!(bnd.fg, i, fcopy)
    else
      bnd.fg[i] = fcopy
    end
end

function move_params!(bnd::BayesNetDAI, node, backup=true)
    bnd.dirty = true

    fac = bnd.fg[node]
    curr = bnd.verts[node]
    pars = vars(fac) - curr

    for p in 1:nrStates(pars) # FIXME: hardcoded binary
        newval = rand()
        c1,c2 = conditionalStateBoth(curr, pars, p)
        fac[c1] = newval
        fac[c2] = 1 - newval
    end
    if backup
      setBackedFactor!(bnd.fg, node, fac)
    else
      bnd.fg[node] = fac
    end
    return 0.0 # FIXME hardcoded binary: new_count * lgamma(self.arity)
end

function adjust_factor!(bnd::BayesNetDAI, node, addlist, dellist, backup=true)
    for n in dellist
        del_edge!(bnd, n, node, backup)
    end
    for n in addlist
        add_edge!(bnd, n, node, backup)
    end
    return 0.0 # FIXME hardcoded binary: new_count * lgamma(self.arity)
end

function show(io::IO, bnd::BayesNetDAI)
    show(io, bnd.verts)
    show(io, bnd.fg)
    #print(io, " : [")
    #varvec = vars(vars((n.fac)))
    #for v in varvec
        #print(io, "$(label(v)),")
    #end
    #print(io, "\b] : [")
    #for i in 1:nrStates(n.fac)
        #@printf io "%.2f," n.fac[i]
    #end
    #print(io, "\b]")
end

type BayesNetSampler <: Sampler
    bnd::BayesNetDAI
    data::Matrix{Int}
    changelist::Vector{Int}
    x::Vector{Int} # node ordering/mapping
    oldx::Vector{Int}
    fvalue::Vector{Float64}
    oldfvalue::Vector{Float64}
    logqfactor::Float64
    mat::Matrix{Bool}
    oldmat::Matrix{Bool}
    limparent::Int
    prior::Function
end

function BayesNetSampler(net::BayesNetDAI, data::Matrix{Int}, prior=basic_prior::Function)
    assert(all(0 .< data .<= 2)) # FIXME binary hardcode
    n = length(net)
    BayesNetSampler(net,
    data,
    [1:n;],
    randperm(n),
    randperm(n),
    zeros(Float64, n),
    zeros(Float64, n),
    0.0,
    zeros(Bool,n,n),
    zeros(Bool,n,n),
    4,
    prior)
end

function BayesNetSampler(n::Int, data::Matrix{Int}, prior=basic_prior::Function)
    BayesNetSampler(BayesNetDAI(n), data, prior)
end

function show(io::IO, bns::BayesNetSampler)
    show(io, bns.mat[bns.x,bns.x])
end

function energy(bns::BayesNetSampler)

    eprior = bns.prior(bns)

    ##### energy from data (likelihood) ####
    counts = {}
    buf = Int[]
    for node in bns.changelist
      fac = bns.bnd.fg[node]
      facvs = vars(fac)
      inds = labels(fac)
      count = zeros(Int,nrStates(fac))
      for row=1:size(bns.data,2)
          for ind in inds
              push!(buf, bns.data[ind,row])
          end
          index = calcLinearState(facvs, buf)
          empty!(buf)
          count[index] += 1
      end
      push!(counts, count)
    end

    for i=1:length(counts) # Calculate likelihood
        accum = 0.0
        node = bns.changelist[i]
        v = bns.bnd.verts[node]
        fac = bns.bnd.fg[node]
        pars = vars(fac) - v
        numparstates = nrStates(pars)
        for parstate in 1:numparstates
            i1,i2 = conditionalStateBoth(v, pars, parstate) # Hardcoded binary
            accum += log(fac[i1]) * counts[i][i1]
            accum += log(fac[i2]) * counts[i][i2]
        end
        bns.fvalue[node] = -accum
    end
    return sum(bns.fvalue) - bns.logqfactor + eprior # TODO - or + logqfactor here?
end

function propose!(bns::BayesNetSampler)
    bns.bnd.dirty = true
    clearBackups!(bns.bnd.fg)
    bns.oldx[:] = bns.x
    bns.oldfvalue[:] = bns.fvalue
    bns.oldmat[:] = bns.mat

    bns.logqfactor = 0.0
    resize!(bns.changelist, 0)

    bns.logqfactor = 0
    addlist = (Int=>Vector{Int})[]
    dellist = (Int=>Vector{Int})[]
    x = bns.x # convenience labels
    mat = bns.mat
    numnodes = length(bns.bnd.verts)
    scheme = rand(1:3)

    if scheme==1 # Temporal order change
        k = rand(1:numnodes-1)

        #For any parents of node k not shared with k+1 (and vice versa):
        for i=1:k-1
            if mat[i,k] && !mat[i,k+1]
                push!(get!(dellist, x[k], []), x[i])
                push!(get!(addlist, x[k+1], []), x[i])
            elseif !mat[i,k] && mat[i,k+1]
                push!(get!(dellist, x[k+1], []), x[i])
                push!(get!(addlist, x[k], []), x[i])
            end
        end
        push!(bns.changelist, x[k], x[k+1])

        #For any children of node k not shared with k+1 (and vice versa):
        for j=k+2:numnodes
            if mat[k,j] && !mat[k+1,j]
                push!(get!(dellist, x[j], []), x[k])
                push!(get!(addlist, x[j], []), x[k+1])
                push!(bns.changelist, x[j])
            elseif !mat[k,j] && mat[k+1,j]
                push!(get!(dellist, x[j], []), x[k+1])
                push!(get!(addlist, x[j], []), x[k])
                push!(bns.changelist, x[j])
            end
        end

        if mat[k,k+1]
            push!(get!(dellist, x[k+1], []), x[k])
            push!(get!(addlist, x[k], []), x[k+1])
        end

        bns.x[k], bns.x[k+1] = bns.x[k+1], bns.x[k]

    elseif scheme==2 # skeletal change
        overparlimit = true
        local i, j, edgedel
        while overparlimit
            i = rand(1:numnodes)
            j = rand(1:numnodes)
            while i==j
                j = rand(1:numnodes)
            end
            if i > j
                i,j = j,i
            end

            edgedel = mat[i,j]
            mat[i,j] = !edgedel
            overparlimit = sum(mat[:,j]) > bns.limparent
            if overparlimit
                mat[i,j] = edgedel # reset for another loop
            end
        end
        push!(bns.changelist, x[j])

        if edgedel
            push!(get!(dellist, x[j], []), x[i])
        else
            push!(get!(addlist, x[j], []), x[i])
        end

    elseif scheme==3 # Null (parameter only) move

        k = rand(1:numnodes)

        bns.logqfactor = move_params!(bns.bnd, bns.x[k])
        push!(bns.changelist, x[k])
    end

    if scheme != 3
        for node in union(Set(keys(addlist)), Set(keys(dellist)))
            bns.logqfactor += adjust_factor!(bns.bnd, node,
                get(addlist, node,[]),
                get(dellist, node,[]))
        end
    end
    return scheme
end

function reject!(bns::BayesNetSampler)
    bns.bnd.dirty = true
    bns.x[:] = bns.oldx
    bns.fvalue[:] = bns.oldfvalue
    bns.mat[:] = bns.oldmat
    restoreFactors!(bns.bnd.fg)
end

function save!(bns::BayesNetSampler)
    bns.oldx[:] = bns.x
    bns.oldfvalue[:] = bns.fvalue
    bns.oldmat[:] = bns.mat
    clearBackups!(bns.bnd.fg)
end

function record(bns::BayesNetSampler)
  return deepcopy(bns.bnd)
end

function check_bns(bns::BayesNetSampler)
    # Sanity check over topological node ordering
    n = size(bns.mat,1)
    for i=1:n, j=1:i
        assert(!bns.mat[i,j])
    end

    # Sanity checks over bayesnetdai
    check_bnd(bns.bnd)
    #@show bns.mat[bns.x,bns.x]
    #@show net2mat(bns.bnd)
    ord = sortperm(bns.x)
    Base.Test.@test all(bns.mat[ord,ord] .== net2mat(bns.bnd))
end

function check_bnd(bnd::BayesNetDAI)
    nedges = 0
    for i=1:numFactors(bnd.fg)
        v = bnd.verts[i]
        Base.Test.@test v in bnd.fg[i]
        arity = states(v)
        parvars = vars(bnd.fg[i]) - v
        nedges += length(parvars) + 1
        numparstates = nrStates(parvars)
        for parstate=1:numparstates
            tot = 0.0
            for state=1:arity
                index = conditionalState(v, parvars, state, parstate)
                amt = bnd.fg[i][index]
                Base.Test.@test 0.0 <= amt <= 1.0
                tot += amt
            end
            Base.Test.@test_approx_eq tot 1.0
        end
    end
    Base.Test.@test nedges == numEdges(bnd.fg)
end

function controllability(bnd::BayesNetDAI, target_vars, target_states::Vector{Int}, control_var)
  # Measure the amount of shift in probabliity we obtain by fixing control_var
  # to either 0 or 1 in the prob. of target_vars being in target_states
  @assert typeof(target_states) == Vector{Int}
  clearBackups!(bnd)
  conbnd = deepcopy(bnd) # Don't want to mutate incoming bnd

  target_vars = DAI.VarSet(conbnd.verts[target_vars]...)
  cvar = conbnd.verts[control_var]
  cvarset = DAI.VarSet(cvar)

  tarstate = calcLinearState(target_vars, target_states)

  # @show conbnd.fg
  conbnd.fg[control_var] = Factor(cvar)

  set_factor!(conbnd, control_var, [0.0, 1.0], false)
  setupIA!(conbnd)
  clearBackups!(conbnd)
  # @show conbnd.fg
  jointfac = marginal(conbnd.ia, target_vars, true)
  # normalize!(jointfac)

  ssd1 = jointfac[tarstate]
  # @show "ssd1", jointfac

  set_factor!(conbnd, control_var, [1.0, 0.0], false)
  setupIA!(conbnd)
  clearBackups!(conbnd)
  jointfac = marginal(conbnd.ia, target_vars, true)
  # normalize!(jointfac)

  ssd2 = jointfac[tarstate]
  # @show "ssd2", jointfac

  setupIA!(bnd)
  clearBackups!(bnd)
  jointfac = marginal(bnd.ia, target_vars)

  origssd = jointfac[tarstate]
  # @show "origssd", jointfac

  # @show ssd1, ssd2, origssd
  return max(ssd1,ssd2) - origssd
end

function ting_cod(bns::BayesNetSampler, target::Int, predictors::Vector{Int})
  data = bns.data
  # From Ting IEEE http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6190164
  # Eq 11
  N = size(data, 2)
  N1 = sum(data[target, :]-1)
  N0 = N - N1
  k = length(predictors)

  epsilon = 0.0
  for i=1:2^k
    states = digits(i-1, 2, k)
    states += 1
    U = V = 0
    for rowind=1:N
      if data[predictors, rowind] == states
        U += 1-(data[target, rowind]-1)
        V += data[target, rowind]-1
      end
    end
    epsilon += min(U, V)
  end
  if min(epsilon, N0, N1) == 0
    return 1
  end
  epsilon /= N
  return 1 - epsilon/min(N0/N, N1/N)
end

function cod(bnd::BayesNetDAI, target::Int, predictors::Vector{Int})
  clearBackups!(bnd)
  setupIA!(bnd)
  targv = bnd.verts[target]
  targvs = VarSet(targv)
  predvs = VarSet(bnd.verts[predictors]...)

  margfac = belief(bnd.ia, targvs)
  joint = marginal(bnd.ia, targvs+predvs)

  # normalize!(margfac)
  # normalize!(joint)

  numparstates = nrStates(predvs)
  epsdot = 0.0
  for parstate=1:numparstates
    i1,i2 = conditionalStateBoth(targv, predvs, parstate) # Hardcoded binary
    epsdot += min(joint[i1], joint[i2])
  end
  epszero = min(margfac[1], margfac[2])

  # if epszero == epsdot == 0.0
  #   return 0.0
  if (1.0 - epsdot/epszero) == -Inf || epszero == 0.0 || epsdot == 0.0
    println("#######################")
    @show bnd.fg
    @show joint
    @show margfac
    @show epsdot, epszero, 1.0 - epsdot/epszero
    check_bnd(bnd)
    @show numEdges(bnd.fg)
    # bnd.ia = InfAlg(copy(bnd.fg), "BP[inference=SUMPROD,verbose=1,updates=SEQMAX,logdomain=0,tol=1e-9,maxiter=10000,damping=0.0]")
    bnd.ia = InfAlg(copy(bnd.fg))
    init!(bnd.ia)
    run!(bnd.ia)
    margfac = belief(bnd.ia, targvs)
    joint = marginal(bnd.ia, targvs+predvs)
    @show joint
    @show margfac
    error("")
    return 0.0
  else
    return 1.0 - epsdot/epszero
  end
end

function setupIA!(bnd::BayesNetDAI)
  if bnd.dirty
    bnd.ia = InfAlg(bnd.fg)
    init!(bnd.ia)
    run!(bnd.ia)
  end
end

function kld(bnd::BayesNetDAI, other::BayesNetDAI)
    bnd.memo_entropy = entropy(bnd)
    # Now entropy and iaree are correct
    accum = 0.0
    for i=1:numFactors(bnd.fg)
      v = other.verts[i]
      arity = states(v)
      factor = other.fg[i]
      numparstates = div(nrStates(factor),arity)
      # Assumption: same ordering between bnd and other
      jointvars = vars(factor)
      parvars = jointvars - v
      joint = marginal(bnd.ia, jointvars)
      for parstate=1:numparstates
        i1,i2 = conditionalStateBoth(v, parvars, parstate) # Hardcoded binary
        accum += joint[i1] * log2(factor[i1])
        accum += joint[i2] * log2(factor[i2])
      end
    end
    return -bnd.memo_entropy - accum
end

function factor_klds(bnd1::MCBN.BayesNetDAI, bnd2::MCBN.BayesNetDAI, node::Int)
    bnd2.memo_entropy = MCBN.entropy(bnd2) # Prep the InfAlg
    accum = 0.0
    v = bnd1.verts[node]
    arity = states(v)
    factor = bnd1.fg[node]
    jointvars = vars(factor)
    parvars = jointvars - v
    numparstates = nrStates(parvars)
    joint = marginal(bnd2.ia, jointvars)
    for parstate=1:numparstates
      i1,i2 = conditionalStateBoth(v, parvars, parstate) # Hardcoded binary
      joint1 = joint[i1]
      joint2 = joint[i2]
      jointsum = joint1+joint2
      accum += factor[i1] * log2(factor[i1]/joint1*jointsum)
      accum += factor[i2] * log2(factor[i2]/joint2*jointsum)
    end
    return accum
end

function entropy(bnd::BayesNetDAI)
    if !bnd.dirty
        return bnd.memo_entropy
    end
    setupIA!(bnd)
    accum = 0.0
    for i in 1:numFactors(bnd.fg)
        v = bnd.verts[i]
        arity = states(v)
        pars = vars(bnd.fg[i]) - v
        numparstates = nrStates(pars)
        margfac = marginal(bnd.ia, pars)
        factor = bnd.fg[i]
        for parstate in 1:numparstates
            accumsum = 0.0
            i1,i2 = conditionalStateBoth(v, pars, parstate) # Hardcoded binary
            temp = factor[i1]
            accumsum -= temp == 0.0 ? 0.0 : temp * log2(temp)
            temp = factor[i2]
            accumsum -= temp == 0.0 ? 0.0 : temp * log2(temp)
            accum += margfac[parstate] * accumsum
        end
    end
    bnd.memo_entropy = accum
    bnd.dirty = false
    return accum
end

function naive_entropy(bnd::BayesNetDAI)
    numstates = 1
    accum = 0.0
    tot = 0.0
    allvars = VarSet(bnd.verts...)
    for n in bnd.verts
        numstates *= states(n)
    end
    assert(numstates == nrStates(allvars))

    for i=1:numstates
        accumprod = 1.0
        statemap = calcState(allvars, i)
        for j in 1:numFactors(bnd.fg)
            reducedstate = calcLinearState(vars(bnd.fg[j]), statemap, allvars)
            accumprod *= bnd.fg[j][reducedstate]
        end
        tot += accumprod
        accum += accumprod == 0 ? 0.0 : accumprod * log2(accumprod)
    end
    assert((tot-1.0)<=eps())
    return -accum
end

function basic_prior(bns::BayesNetSampler)
    0.0
end

#need to handle prior information: structure and cpd
#init cpds
#init structures


text:
using StructuresKit

#Verify Julia solution against Mathematica solution from
#Example 2 from Plaut, R.H., Moen, C.D.(2020). "Lateral-Torsional Deformations of C-Section and Z-Section Beams with Continuous Bracing".  Proceedings of the Structural Stability Research Council Annual Conference, Atlanta, Georgia.

#Z-section, no slope, simple span
#kϕ=300 N*mm/rad/mm, kx=0, gravity load
#This solution was calculated with Mathematica
#  max ϕ    q
plaut_solution=[0.0 0.0
0.0221 0.2
0.0476 0.4
0.0771 0.6
0.1116 0.8
0.1313 0.9
0.1527 1
0.1763 1.1
0.2022 1.2
0.2310 1.3
0.2630 1.4
0.2989 1.5
0.3393 1.6
0.3615 1.65
0.3756 1.68
0.3853 1.7
0.4003 1.73]

#*********** Julia solution

#inputs

#Ix Iy Ixy J Cw
section_properties = [(3.230E6,449530,-865760, 397.09, 3.4104E9)]

#E  ν
material_properties = [(200,0.30)]

#ax ay
load_location = [(27.826,101.6)]

#kx kϕ
spring_stiffness = [(0.0,300/1000)]

#ay_kx
spring_location = [(101.6)]

#member information
#L dL section_properties material_properties load_location spring_stiffness spring_location
member_definitions = [(7620, 7620/12 ,1, 1, 1, 1, 1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 7620]

#initialize twist vector
ϕmax=zeros(length(plaut_solution[:,1]))

#calculate max twist with PlautBeam
for i=1:length(plaut_solution[:,1])

    load = (0.0, plaut_solution[i,2]/1000)  #kN/mm

    z, u, v, ϕ, beam_properties = ThinWalledBeam.solve(member_definitions, section_properties, material_properties, spring_stiffness, spring_location, supports, load, load_location, end_boundary_conditions)

    ϕmax[i] = maximum(ϕ)

end

#calculate percent error at every load compared to Mathematica solution
error= abs.((ϕmax .- plaut_solution[:,1]))./ plaut_solution[:,1]

max_error=maximum(error[2:end])

#accept 1% error from numerical solution
@test max_error <= 0.01

#update with this format!!

using StructuresKit

#Ix Iy Ixy J Cw
section_properties = [(3.230E6,449530,-865760, 397.09, 3.4104E9)]

#E  ν
material_properties = [(200,0.30)]

#kx kϕ
spring_stiffness = [(0.0,300/1000)]

#ay_kx
spring_location = [(101.6)]

#qx qy
loads = [(0.0001, 0.00002)]

#ax ay
load_locations = [(27.826,101.6)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 7620]

#member information
#L(1) dL(2) section_properties(3) material_properties(4) spring_stiffness(5) spring_location(6) load(7) load_location(8)
member_definitions = [(7620, 7620/12, 1, 1, 1, 1, 1, 1)]

#Define model inputs.
model = ThinWalledBeam.user_interface(member_definitions, section_properties, material_properties, spring_stiffness, spring_location, loads, load_locations, end_boundary_conditions, supports)

#Calculate model stiffness and external force.
model = ThinWalledBeam.define(model)

#Solve model.
model = ThinWalledBeam.solve(model)


text:
using WebSockets
using Compat; import Compat.String
using Base.Test

import WebSockets: generate_websocket_key,
                   write_fragment,
                   read_frame,
                   is_websocket_handshake

import HttpCommon: Request

#is_control_frame is one line, checking one bit.
#get_websocket_key grabs a header.
#is_websocket_handshake grabs a header.

#generate_websocket_key makes a call to a library.
@test generate_websocket_key("dGhlIHNhbXBsZSBub25jZQ==") == "s3pPLMBiTxaQ9kYGzzhZRbK+xOo="

# Test writing

function xor_payload(maskkey, data)
    out = Array(UInt8, length(data))
  for i in 1:length(data)
    d = data[i]
    d = d $ maskkey[mod(i - 1, 4) + 1]
    out[i] = d
  end
  out
end

const io = IOBuffer()

# Length less than 126
for len = [8, 125], op = (rand(UInt8) & 0b1111), fin=[true, false]

    test_str = randstring(len)
    write_fragment(io, fin, test_str, op)

    frame = takebuf_array(io)

    @test bits(frame[1]) == (fin ? "1" : "0") * "000" * bits(op)[end-3:end]
    @test frame[2] == @compat UInt8(len)
    @test String(frame[3:end]) == test_str

    # Check to see if reading message without a mask fails
    in_buf = IOBuffer(frame)
    @test_throws ErrorException read_frame(in_buf)
    close(in_buf)

    # add a mask
    maskkey = rand(UInt8, 4)
    data = vcat(
        frame[1],
        frame[2] | 0b1000_0000,
        maskkey,
        xor_payload(maskkey, frame[3:end])
    )
    frame_back = read_frame(IOBuffer(data))

    @test frame_back.is_last == fin
    @test frame_back.rsv1 == false
    @test frame_back.rsv2 == false
    @test frame_back.rsv3 == false
    @test frame_back.opcode == op
    @test frame_back.is_masked == true
    @test frame_back.payload_len == len
    @test all(map(==, frame_back.maskkey, maskkey))
    @test test_str == String(frame_back.data)
end

# Length 126 or more
for len = 126:129, op = 0b1111, fin=[true, false]

    test_str = randstring(len)
    write_fragment(io, fin, test_str, op)

    frame = takebuf_array(io)

    @test bits(frame[1]) == (fin ? "1" : "0") * "000" * bits(op)[end-3:end]
    @test frame[2] == 126

    @test bits(frame[4])*bits(frame[3]) == bits(hton(@compat UInt16(len)))

    # add a mask
    maskkey = rand(UInt8, 4)
    data = vcat(
        frame[1],
        frame[2] | 0b1000_0000,
        frame[3],
        frame[4],
        maskkey,
        xor_payload(maskkey, frame[5:end])
    )
    frame_back = read_frame(IOBuffer(data))

    @test frame_back.is_last == fin
    @test frame_back.rsv1 == false
    @test frame_back.rsv2 == false
    @test frame_back.rsv3 == false
    @test frame_back.opcode == op
    @test frame_back.is_masked == true
    @test frame_back.payload_len == len
    @test all(map(==, frame_back.maskkey, maskkey))
    @test test_str == String(frame_back.data)
end

# TODO: test for length > typemax(Uint32)

#
close(io)

# Tests for is_websocket_handshake
chromeheaders = @compat Dict{AbstractString,AbstractString}(
        "Connection"=>"Upgrade",
        "Upgrade"=>"websocket"
    )
chromerequest = Request(
    "GET",
    "",
    chromeheaders,
    ""
    )

firefoxheaders = @compat Dict{AbstractString,AbstractString}(
        "Connection"=>"keep-alive, Upgrade",
        "Upgrade"=>"websocket"
    )

firefoxrequest= Request(
    "GET",
    "",
    firefoxheaders,
    ""
    )

handler = WebSocketHandler(x->x); #Dummy handler

for request in [chromerequest, firefoxrequest]
    @test is_websocket_handshake(handler,request) == true
end